Characterization of the complete mitochondrial genome of Otus lettia: exploring the mitochondrial evolution and phylogeny of owls (Strigiformes)

Abstract

Large-scale molecular phylogenetic studies of the avian order Strigiformes have been performed, and numerous mitochondrial genomes have been determined. However, their intergeneric relationships are still controversial, and few comprehensive comparative analyses of mitochondrial genomes have been conducted on Strigiformes. In this study, the mitochondrial genome of Otus lettia was determined and compared with other Strigiformes. The O. lettia mitochondrial genome was 16,951 bp in size. For Strigiformes, atp8 can be used as a suitable molecular marker for population genetic diversity, while cox1 is a candidate barcoding marker for species identification. All protein-coding genes may be under strong purifying selection pressure, and one extra cytosine insertion located in nad3 is common to all owls except Tyto longimembris, T. alba, and Athene noctua. Four different mitochondrial gene arrangement types were found among the Strigiformes mitogenomes, and their evolutionary relationship between each other can be perfectly explained by the tandem duplication and random loss model. The phylogenetic topologies using the mitochondrial genomes showed that target species O. lettia had a closer relationship with O. scops + O. sunia than O. bakkamoena, the genus Glaucidium was paraphyletic, and the Ninox clade was located at the basal position of Strigidae lineage. Our phylogenetic trees also supported the previous recommendations that Sceloglaux albifacies, Ciccaba nigrolineata, and Ketupa flavipes should be transferred to Ninox, Strix, and Bubo, respectively. These findings will be helpful in further unraveling the mitochondrial evolution and phylogeny of Strigiformes.

Keywords:

• Strigiformes
• mitochondrial gene rearrangement
• extra cytosine insertion
• molecular phylogeny

1. Introduction

Birds are primarily diurnal, but one notable exception is the nocturnal owls (Aves: Strigiformes) (Zhou et al. 2019b). Owls display a cosmopolitan distribution (except the Antarctic, most of Greenland, and some remote islands) (Wink et al. 2009; Zhou et al. 2019b). According to the latest data, the order Strigiformes comprises more than 240 recognized species assigned to two families: Tytonidae (barn owls) and Strigidae (typical owls) (Handbook of the Birds of the World and BirdLife International 2019; Gill et al. 2021). Although recent large-scale molecular phylogenetic studies have been performed on the Strigiformes, most of them relied on a single or few markers (Wink and Heidrich 2000; Wink et al. 2009; Sun et al. 2016; Zhang et al. 2016), moreover, the intergeneric relationships inferred by these studies remain controversial (Wink et al. 2009; Kang et al. 2018; Lee et al. 2018; Spiridonova and Surmach 2018; Park et al. 2019a, 2019b; Zhou et al. 2019a; Sun et al. 2020; Zhong et al. 2020), and extensive paraphyly was revealed in Strigidae (Salter et al. 2020; Sun et al. 2020). Therefore, the phylogenetic relationships of Strigiformes species still need further studies.

Owing to several features, including multiple copies in a cell, extremely compact gene packing, maternal inheritance, low sequence recombination, and generally rapid evolutionary rates, mitochondrial genome (mitogenome) has been used extensively in studies of phylogenetic relationship and comparative genome (Liu et al. 2014; Yan et al. 2017; Sun et al. 2020). Dozens of Strigiformes mitogenomes have been deposited in the GenBank now, and many studies have successively reported their characteristics (Zhang et al. 2016; Zhong et al. 2020). Meanwhile, several mitochondrial gene arrangement orders have been found, such as duplicated control regions (Hanna et al. 2017; Kang et al. 2018). Comparative mitogenomic analyses among multiple related species can provide sufficient sequence data (e.g. gene organization, arrangement, and codon usage) and it is of great significance for understanding biodiversity and evolution (Bondarenko et al. 2019). Yet, few comprehensive comparative mitogenome analyses have been performed on Strigiformes. Hanna et al. (2017) confirmed the presence of a duplicated control region in Strigiformes mitogenome for the first time, Kang et al. (2018) conducted the mitochondrial gene rearrangements analyses of Strigidae for the first time and determined four gene orders. As more Strigiformes mitogenomes are sequenced, will we find more novel rearrangement types?

In this study, we sequenced the mitogenome of the Collared Scops Owl (Otus lettia), which is listed as National Class II Key Protected Wild Animals in China. Based on our new data together with the other 32 available Strigiformes mitogenomes from GenBank, we aimed to: (i) determine the comparative mitogenomic features of owls, including genome size, nucleotide composition, nucleotide diversity, gene evolutionary rate, features of a protein-coding gene, and gene arrangement; (ii) reconstruct the phylogenetic relationships for owls using mitogenomes.

2. Materials and methods

2.1. Specimen collection and DNA extraction

The sample used in this experiment was an illegally hunted owl (Location: Zhaotong City, Yunnan Province, China), and it was confiscated by the forest police and identified as Otus lettia erythrocampe by morphological characters and molecular identification. The forest police provided us with a few individual feathers for this study, which we then preserved in 95% ethanol (Specimen ID PB2019-211-1, presently stored at the Chengdu Research Base of Giant Panda Breeding, contact person Jiaojiao Yu, email [email protected]). Genomic DNA was extracted from the feather samples using the QIAGEN^® DNeasy Blood and Tissue Kit, following the manufacturer’s instructions.

2.2. DNA amplification and sequencing

Overlapping fragments were amplified by the ordinary PCR method. PCR primer sets were chosen from a previous study (Sorenson et al. 1999). Reaction volumes and cycling conditions of PCR were shown in Supplementary Table 1 and Supplementary Figure 1, respectively. The DNA purification and sequencing were carried out as previously described by Liu et al. (2019). The sequences were checked using Chromas (http://technelysium.com.au/) and assembled using SeqMan from the DNASTAR Lasergene package (https://www.dnastar.com/), then checked manually.

2.3. Genome annotation

The preliminary predictions of mitochondrial ribosomal RNA genes (rRNAs), transfer RNA genes (tRNAs), protein-coding genes (PCGs) and control region (D-loop), and the definitions of gene boundaries were carried out using MITOS (http://mitos2.bioinf.uni-leipzig.de/index.py) (Bernt et al. 2013) and MitoAnnotator (http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html) (Iwasaki et al. 2013). The PCGs were rechecked using NCBI ORFfinder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). The tRNAs were further confirmed with ARWEN (http://mbio-serv2.mbioekol.lu.se/ARWEN/) (Canbäck and Laslett 2008) and tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan 2016) via their proposed cloverleaf secondary structures and anticodon sequences. Additionally, all results were verified via alignment with their closely related species’ homologs to further improve the genome annotation. The mitogenome map was drawn by the online tool CGView Server (http://cgview.ca/) (Grant and Stothard 2008).

2.4. Data acquisition and processing

Over 50 Strigiformes mitogenomes have been deposited in GenBank now, but many of them belong to the same species, thus we retrieved 32 representative mitogenomes belonging to 27 species (Table 1). The taxonomic classification of all species was following the current version of HBW and BirdLife Taxonomic Checklist (Handbook of the Birds of the World and BirdLife International 2019). Our preliminary investigation found that the annotations of nad3 genes in the sequences KP889214 and MG916810 were inconsistent with the descriptions of their corresponding studies (Zhang et al. 2016; Lee et al. 2018). Thus, we re-annotated all the nad3 genes from all retrieved Strigiformes mitogenomes. Finally, we found nad3 sequences of 12 mitogenomes showing differences between original annotation and re-annotation results. The comparison of the original annotation and re-annotation results of these 12 sequences were shown in Supplementary Table 2.

Table 1. Mitogenome sequences were used in this study.

Species (scientific name)	Sequence (accession ID)	Size (bp)	Usage		Notes
			Comparative analyses^a	Phylogenetic analyses
Strigidae
Otus lettia	MW364567	16,951	√	√	Complete, this study
Otus sunia	MN276053	17,604	√	√	Complete
Otus sunia	MF346692	17,835	√	√	Complete
Otus scops	KY471456	17,317	√	√	Complete, re-annotated nad3
Otus scops	KT340630	17,413	√	√	Complete
Otus bakkamoena	KT340631	17,389	√	√	Complete
Asio flammeus	KP889214	18,966	√	√	Complete, re-annotated nad3
Asio otus	MG916810	17,555	√	√	Complete, re-annotated nad3
Strix aluco	MN122823	16,490		√	Partial, trnP-nad6-trnE is unconfirmed, D-loop is partial, re-annotated nad3
Strix uralensis	MG681081	18,708	√	√	Complete
Strix uralensis	MG681082	18,709	√	√	Complete
Strix occidentalis	MF431746	19,889	√	√	Complete, re-annotated nad3
Strix varia	MF431745	18,975	√	√	Complete, re-annotated nad3
Strix leptogrammica	KC953095	16,307		√	Complete, nad5-trnF may be incorrect
Bubo scandiacus	MG681084	18,734	√	√	Complete
Bubo bubo	MG681083	18,952	√	√	Complete
Bubo bubo	MK656281	18,954	√	√	Complete
Bubo blakistoni	LC099104	19,484	√	√	Complete, re-annotated nad3
Ketupa flavipes	LC099100	19,447	√	√	Complete, re-annotated nad3
Glaucidium brodiei	KP684122	17,318	√	√	Complete
Glaucidium cuculoides	KY092431	17,392	√	√	Complete
Glaucidium brasilianum	MN356303	17,717	√	√	Complete
Athene noctua	MN651450	17,772		√	Complete, atp8 is unverified
Athene noctua	MN122903	15,776	√	√	Partial, D-loop is unconfirmed, re-annotated nad3
Athene brama	KF961185	16,194	√	√	Partial, D-loop and trnF are partial, re-annotated nad3
Ninox scutulata	KT943750	16,208	√	√	Complete
Ninox strenua	KX529654	16,206	√	√	Complete
Ninox novaeseelandiae	AY309457	16,223	√	√	Partial, D-loop is partial
Ciccaba nigrolineata	MN356178	14,875		√	Partial, trnP-nad6-trnE-D-loop is unconfirmed
Sceloglaux albifacies	KX098448	15,565	√	√	Partial, D-loop is unconfirmed
Tytonidae
Tyto alba	EU410491	16,148	√	√	Partial, D-loop and trnF are partial, re-annotated nad3
Tyto longimembris	KP893332	18,466	√	√	Partial, D-loop is partial
Phodilus badius	KF961183	17,086	√	√	Partial, D-loop is partial, re-annotated nad3
Outgroups (Coliiformes, Coliidae)
Colius striatus	MN356125	17,283		√	Complete
Urocolius indicus	MN356373	17,267		√	Complete

^aSince atp8 of MN651450 has not been verified, cytb of KC953095 may be incorrect, and nad6 of MN356178 and MN122823 have not been identified due to incomplete sequencing, these four mitogenome sequences were excluded from the comparative genome analyses (nucleotide composition, A + T content, AT skew, GC skew, statistics of variation sites, nucleotide diversity, and Ka/Ks value).

2.5. Comparative mitogenomic analyses

Twenty-nine mitogenome sequences were used in the comparative genome analyses (Table 1). The nucleotide composition and substitution were calculated using MEGA X (Kumar et al. 2018). The nucleotide skew was measured using the formulas ‘AT-skew = (A − T)/(A + T)’ and ‘GC-skew = (G − C)/(G + C)’ (Perna and Kocher 1995). The nucleotide polymorphism and the ratios of non-synonymous substitutions rate (Ka) to synonymous substitutions rate (Ks) for each PCG among Strigiformes mitogenomes were estimated by DnaSP (Rozas et al. 2017).

2.6. Molecular phylogenetic analyses

The phylogeny of 33 Strigiformes species was constructed using two rRNAs and 10 PCGs excluding nad6, atp8, and cytb due to unverified nad6 in Ciccaba nigrolineata MN356178 and Strix aluco MN122823, atp8 in Athene noctua MN651450, cytb in Strix leptogrammica KC953095 (Table 1). According to genome-scale phylogeny results (Prum et al. 2015), the Coliiformes species Urocolius indicus and Colius striatus were chosen as the outgroups.

The rRNAs sequences were aligned in batches with MAFFT (Katoh and Standley 2013) integrated into PhyloSuite (Zhang et al. 2020), using ‘–auto’ strategy and ‘normal’ alignment mode. The PCGs sequences were aligned with the codon-aware program MACSE (Ranwez et al. 2018) integrated into PhyloSuite (Zhang et al. 2020), based on ‘Vertebrate mitochondrial genetic code.’ Subsequently, gap sites were recognized and removed with trimAl (Capella-Gutiérrez et al. 2009) integrated into PhyloSuite (Zhang et al. 2020), using the ‘automated1’ command. Finally, the trimmed alignments of 12 chosen genes (10 PCGs: atp6, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5; two rRNAs: 12S, 16S) were concatenated into two datasets using PhyloSuite (Zhang et al. 2020), 10PCG dataset (9531 bp) and 10PCG + 2rRNA dataset (12,017 bp), which were used for the Maximum Likelihood (ML) and Bayesian inference (BI) analyses.

The ML analyses were carried out using IQ-TREE (Minh et al. 2020) under Edge-linked partition style for 10,000 replicates of ultrafast bootstrap (Hoang et al. 2018), as well as the Shimodaira–Hasegawa–like approximate likelihood-ratio test (SH-aLRT) (Guindon et al. 2010). The BI phylogenies were inferred using MrBayes (Ronquist et al. 2012) under partition model (2 parallel runs, 10,000,000 generations, the minimum estimated sample size values were more than 100, the average standard deviation of split frequencies values were below 0.01, and the potential scale reduction factor values were close to 1.0000), in which the initial 25% of sampled data were discarded as burn-in. The partition-specific substitution model (Supplementary Table 3) was chosen by ModelFinder (Kalyaanamoorthy et al. 2017) integrated into PhyloSuite (Zhang et al. 2020) according to BIC criterion.

3. Results

3.1. Characterization of the Otus lettia mitogenome

The mitogenome of O. lettia was a circular genome with 16,951 bp in size. The overall nucleotide composition was biased toward C and A (32.86% C, 30.59% A, 22.72% T, and 13.83% G), and G + C content was lower than A + T content. It included the canonical two rRNAs, 13 PCGs, 22 tRNAs, and one D-loop region (Figure 1, Supplementary Table 4). Among these, only nad6 and eight tRNAs (trnQ, trnA, trnN, trnC, trnY, trnS2, trnP, and trnE) were encoded on the light strand. A total of 86 nucleotide intervals were noted at 20 locations, and a total of 29 overlapping nucleotides were dispersed at five gene junctions (Supplementary Table 4). The D-loop region of 1378 bp was located between trnE and trnF genes, and possessed an A + T content of 57.70%. A repeat pattern of 22 bp (ACACACACATACCATCCATCCA) with 2.9 tandem repeats was found at 5′-end of the D-loop region. We deposited the O. lettia mitogenome in GenBank with the accession ID MW364567.

Figure 1. The graphical visualization of Otus lettia mitogenome.

3.2. Comparative mitochondrial genome

3.2.1. Nucleotide composition

We calculated that the average genome size of 29 complete Strigiformes mitogenomes was 18,003 bp (Table 1) and their overall average base composition was 23.10% T, 32.80% C, 30.32% A, and 13.78% G. For the gene region (Figure 2(a)), nad6 had more G base than the other three bases (G > T > 9.99% C ≈ 9.78% A), whereas the other 12 PCGs did the exact opposite (C > A > T > G), which was closely related to the encode strand; in addition, there was not much difference between 12S (31.46% A ≈ 30.94% C > G > T) and 16S (A > C > 19.01% T ≈ 18.87% G). More interestingly, more than 50% A + T content was observed in PCGs and rRNAs except for nad6 and 12S; atp8 had the highest A + T content (Figure 2(b)), which was also found in Charadriiformes mitogenomes (Hu et al. 2017).

Figure 2. The base composition (a) and A + T content (b) of 29 Strigiformes mitogenomes. Each circle represents a mitogenome.

For the nucleotide skew of Strigiformes mitochondrial genes (Supplementary Figure 2), nad6 showed marked negative AT-skew and positive GC-skew values, whereas rRNAs and other PCGs showed negative GC-skew and positive AT-skew values (except for AT-skew of nad1 in Athene brama KF961185, Athene noctua MN122903 and O. bakkamoena KT340631, nad4L in Asio flammeus KP889214, Strix uralensis MG681081 and S. uralensis MG681082, cox1 in O. bakkamoena KT340631 and Glaucidium brodiei KP684122, and nad3 in O. bakkamoena KT340631 was slightly negative).

3.2.2. Gene evolutionary patterns

For the rRNAs and PCGs, the percentage of variable sites was from 41.38% in cox1 to 71.43% in atp8, the nucleotide diversity was from 0.14288 in cox1 to 0.24713 in atp8, there was a similar trend between the two observations (Figures 3(a,b)). For the PCGs, the Ka, Ks, and Ka/Ks of atp8 were the highest, while the Ka and Ka/Ks of cox1 and the Ks of nad4 were the lowest (Figures 3(c–e)). Obviously, differences of Ka basically reflected the differences of Ka/Ks among the 13 PCGs, which was congruent with previous studies (Hu et al. 2017; Bi et al. 2019).

Figure 3. Statistics of variation sites (a), nucleotide diversity (b), Ka (c), Ks (d), and Ka/Ks (e) of 29 Strigiformes mitogenomes.

3.2.3. Feature of nad3

As previously reported (Mindell et al. 1998), we also discovered an extra cytosine insertion causing a frameshift at position 174 of nad3 gene in O. lettia mitogenome (Figure 4), but we found that the nad3 cytosine insertion annotations of some other Strigiformes mitogenomes deposited in GenBank were not carried out as recommended by Mindell et al. (1998). Because of this, we re-annotated them (Supplementary Table 2).

Figure 4. Alignment of the mitochondrial nad3 gene region spanning the extra cytosine insertion for Strigiformes species. Additional three nad3 sequences EU601009, MN356151, and EU601014 were used to show the intraspecific differences of inserted cytosine.

3.2.4. Mitochondrial gene rearrangements

Our comparison results show that the gene reorganization was located in the fragment from cytb to trnF, and duplication of D-loop was the most prominent rearrangement feature (Figure 5(a)), indicating that the region spanning D-loop was a hotspot of Strigiformes mitogenomes rearrangement. We summarized four gene arrangement types among the Strigiformes mitogenomes (Figure 5(a)), although a few sequences (e.g. Sceloglaux albifacies KX098448, Athene noctua MN122903, Ciccaba nigrolineata MN356178, Strix leptogrammica KC953095, and Strix aluco MN122823) were not determined due to the incomplete sequencing of important fragments (e.g. the part adjacent to the D-loop). Type I, which was a typical avian gene order and considered ancestral for all Aves, so we named it the presumed ancestral avian gene order (Figure 5(a)). Type II, which was characterized by two copies of trnT-trnP-nad6-trnE-D-loop, so we called it the duplicate trnT∼D-loop gene order (Figure 5(a)). Type III was featured with a pseudo trnE (p-trnE) and duplicate D-loop, and type IV was featured with only a duplicate D-loop (Figure 5(a)). Type I was the most prevailing arrangement out of the four ones, and it was found in all Tytonidae and most Strigidae species (Figure 5(b)). Type II and III emerged in Bubo + Ketupa clade, and type IV appeared in Strix + Ciccaba clade (Figure 5(b)).

Figure 5. Strigiformes mitochondrial gene arrangement types (a) and BI and ML phylograms inferred from 10 PCGs and two rRNAs (b). The orders of the genes outside the region depicted in this figure are the same. The phylograms were visualized and annotated using FigTree (http://tree.bio.ed.ac.uk/software/figtree/). Values at nodes are Bayesian posterior probabilities, SH-aLRT bootstrap percentages, and ultrafast bootstrap percentages, respectively. I–IV denote five different gene arrangement types. ‘?’: the arrangement type is uncertain due to incomplete sequencing.

3.3. Mitochondrial phylogenomics

In this study, both ML and BI methods of both datasets (10PCG dataset and 10PCG + 2rRNA dataset) provided identical and well-supported tree topologies and indicated that all owl species grouped into two large lineages, representing two families: Tytonidae and Strigidae (Figure 5(b), Supplementary Figure 3). The Ninox clade was located at the basal position of the Strigidae lineage. The genus Glaucidium was paraphyletic. Otus lettia had a close relationship with O. scops + O. sunia. Sceloglaux albifacies, Ciccaba nigrolineata, and Ketupa flavipes were respectively nested within other genus clades, implying their classification status needs to be concerned.

4. Discussion

4.1. Comparative mitochondrial genome

The O. lettia mitogenome has many features similar to other owls, such as genome size and content, nucleotide composition, and skew. AT and GC skews are used to measure compositional asymmetry (Hu et al. 2017). We found nad6 of Strigiformes have negative AT-skew and positive GC-skew values, which has been observed also in Accipitriformes (Jiang et al. 2015, 2019). However, Hu et al. (2017) found that there were positive AT-skew and negative GC-skew values in nad6 of 20 Charadriiformes mitogenomes. This completely opposite phenomenon caught our attention, and eventually, we detected that Hu et al. calculated the skew values of heavy strand corresponding to nad6 location.

To gain a better understanding of the evolutionary patterns of the mitochondrial gene of Strigiformes, we compared the evolutionary rate (Figure 3). Based on the percentage of variable sites and nucleotide diversity, we can infer that atp8 was the most variable mitochondrial gene in Strigiformes, and cox1 was the most conserved. Without considering gene size, atp8 is suitable as a molecular marker for population genetic diversity analyses, while cox1 is a candidate barcoding marker for species identification. Moreover, in Strigiformes mitogenomes, the highest Ka/Ks value was detected for atp8 and the lowest for cox1, similarly in the Accipitriformes, Charadriiformes, and Piciformes mitogenomes (Jiang et al. 2015; Hu et al. 2017; Bi et al. 2019). In addition, the Ka/Ks ratios of all 13 PCGs were lower than 1, indicating they may be under strong purifying selection pressure.

Yan et al. (2017) once found that the nad3 of all Strigiformes, Galliformes, and Piciformes birds had an untranslated extra cytosine insertion. Whereas some exceptions were found when more sequences were available for comparative analysis (Li et al. 2014; Bi et al. 2019), including Athene noctua (MN651450 and EU601009), Tyto alba (EU601014), and T. longimembris (KP893332) analyzed in this study (Figure 4). Moreover, the presence or absence of inserted cytosine appeared to randomly happen in closely related species (Jiang et al. 2015; Hu et al. 2017), and even in the same species, such as in Chroicocephalus saundersi (Charadriiformes; KJ631624 and JQ071443), A. noctua (Strigiformes; MN122903, MN651450, and EU601009), and T. alba (Strigiformes; EU410491, MN356151, and EU601014) (Figure 4). However, we cannot rule out whether the interspecies/intraspecific differences are due to mis-editing of the original sequencing files during sequence inspection, so more data are needed in order to verify this. The phenomenon of extra-base insertion was not only common in birds (Li et al. 2014; Jiang et al. 2015; Hu et al. 2017; Yan et al. 2017; Bi et al. 2019), but also in turtles (Mindell et al. 1998). Mindell et al. (1998) speculated the extra-base appears to be not processed during translation. Therefore, we also speculated that the extra cytosine insertion may be the remnant of avian mitogenomes during the rapid explosive evolution (Li et al. 2014).

Kang et al. (2018) described the Strigidae mitochondrial gene orders for the first time. Compared with them, we added more species and mitogenome sequences, aiming to explore gene rearrangement on a larger scale. Among the Strigiformes mitogenomes, we also summarized four gene arrangement types, which were identical to Kang et al. (2018) (Figure 5(a)). The tandem duplication and random loss model can be used to explain most vertebrate mitogenome reorganization (Yang et al. 2018). Based on this model, we hypothesized the evolutionary relationship of four Strigiformes mitochondrial gene arrangement types (Figure 5(a)). Type I generated type II by one tandem duplication, then type II generated type III and type IV by random loss or degeneration. Type I was a typical avian gene order and considered ancestral for all Aves. In our study, Type I was the most prevailing gene order out of the four ones, and it was found in all Tytonidae and most Strigidae species. Type II and III emerged in Bubo + Ketupa clade, and type IV appeared in Strix + Ciccaba clade (Figure 5(b)). We confirmed that Type II has arisen in the common ancestor of Bubo + Ketupa clade and Strix + Ciccaba clade (Hanna et al. 2017). In type II and type III, the pseudogene trnE (p-trnE) and gap degenerated from their functional genes trnE and nad6, respectively. In contrast to trnE, p-trnE lacked anticodon triplet TTC in the anticodon loop. Kang et al. (2018) reported S. leptogrammica KC953095 possessed type I and thought it was a single reversal. However, as previously reported by Hanna et al. (2017), our re-annotation results noted that this sequence potentially contained significant errors in the fragment from nad5 through nad6 to trnF, so we couldn’t determine its correct rearrangement type. Moreover, for the mitochondrial gene order of S. leptogrammica, we thought the single reversal interpretation was debatable. Therefore, it is necessary to collect the samples of S. leptogrammica and re-sequence its mitogenomes in future studies.

4.2. Mitochondrial phylogenomics

Inferring phylogenetic relationships among taxa is a prerequisite for evolutionary and conservation biology studies (Delsuc et al. 2005). The introduction of the use of molecular data in phylogenetics had created splendid progress (Delsuc et al. 2005). However, information from a single gene is often insufficient to obtain robust statistical support for particular nodes of a phylogeny (Delsuc et al. 2005). With the development of sequencing technology, this predicament was solved by genomic data, which could provide a more credible result than a single gene (Song et al. 2015; Liu et al. 2019). In the present study, two datasets (10PCG dataset and 10PCG + 2rRNA dataset) connected from 12 mitochondrial genes were used.

Target species O. lettia had a closer relationship with O. scops + O. sunia than O. bakkamoena, which was inconsistent with previous studies based on cytb dataset (Wink et al. 2009; Zhang et al. 2016). More importantly, the taxonomic status of the three species was a matter of concern. Firstly, Sceloglaux albifacies were nested within the Ninox clade. S. albifacies was an extinct species and was once found in New Zealand. Its taxonomic issue had been widely debated (Wink et al. 2009). Previous molecular phylogenetic study based on mitochondrial cytb also demonstrated this species should be referred to as Ninox albifacies (Wood et al. 2016). Secondly, our results also showed the nested position of Ciccaba nigrolineata within genus Strix, which supports the taxonomic view (Strix nigrolineata) of the Avibase online database (Lepage et al. 2014) and the IOC World Bird List (Gill et al. 2021). Thirdly, our phylogenetic topology indicated the nested position of Ketupa flavipes within the genus Bubo clade. The scientific discussion about the genus Ketupa lasted for a long time. Earlier, the morphological study had recommended that Ketupa was incorporated into Bubo genus (Amadon and Bull 1988), although there were different opinions. Moreover, previous molecular phylogenetic studies based on a single gene (cytb) and the combined dataset (cytb + rag-1) had also proposed merging Ketupa into Bubo (Wink and Heidrich 2000; Wink et al. 2009). However, Spiridonova and Surmach (2018) recently kept the taxonomic status of K. flavipes, K. zeylonensis, and K. ketupu, and suggested that it was more reasonable to assign B. blakistoni along with these three species to a separate Ketupa genus based on mitogenomic phylogeny. Now different taxonomic results are still used in these primary world bird lists including Avibase on-line database (Lepage et al. 2014), Clements Checklist (Clements et al. 2019), HBW and BirdLife Taxonomic Checklist (Handbook of the Birds of the World and BirdLife International 2019), and IOC World Bird List (Gill et al. 2021); for example, B. flavipes was classified into Ketupa genus by HBW and BirdLife Taxonomic Checklist referred by this study (Handbook of the Birds of the World and BirdLife International 2019). According to our phylogenetic topology, we supported the recommendation that merges Ketupa into Bubo.

There was no consistent systematic relationship at the genus level in previous studies. In this study, both ML and BI methods provided identical and well-supported tree topologies (Figure 5(b), Supplementary Figure 3), which are in congruence with Sun et al. (2020). The Ninox clade was located at the basal position of Strigidae lineage, supporting the conclusion of Sun et al. (2016) and Zhou et al. (2019a), but refuting those of other studies (Wood et al. 2016; Spiridonova and Surmach 2018; Park et al. 2019a, 2019b; Zhong et al. 2020). In our present study and several others (Park et al. 2019b; Salter et al. 2020; Sun et al. 2020), the genus Glaucidium was paraphyletic. Given the considerable morphological differences and deep mitogenome divergence between G. brodiei and other Strigidae species, Spiridonova and Surmach had pointed out the need for the revision of the Glaucidium genus and its redescription as a separate family (Spiridonova and Surmach 2018). We recommend against such a taxonomy because this division would result in a paraphyly of the Strigidae lineage.

Author contributions

Jiaojiao Yu and Jiabin Liu: conceptualization (equal); methodology (equal); investigation (equal); data curation (equal); formal analysis (equal); visualization (equal); writing-original draft (equal); writing-review and editing (equal). Cheng Li: investigation (equal); data curation (equal); formal analysis (equal); writing-review and editing (equal). Wei Wu: resources (equal); investigation (equal); data curation (equal); writing-review and editing (equal). Feifei Feng: resources (equal); investigation (equal); writing-review and editing (equal). Qizhi Wang and Xiaofeng Ying: resources (equal); writing-review and editing (equal). Dunwu Qi: conceptualization (equal); writing-review and editing (equal); funding acquisition (supporting). Guilan Qi: conceptualization (equal); formal analysis (equal); validation (lead); project administration (lead); supervision (lead); writing-review and editing (equal). All authors contributed to and approved the final manuscript.

Acknowledgments

We would like to thank Dr. Jacob R. Owens at Los Angeles Zoo and Botanical Gardens for revising the manuscript.

Disclosure statement

The authors declare that they have no competing interests.

Additional information

Funding

This study was funded by the National Natural Science Foundation of China [31772484], the Science and Technology Department of Sichuan Province [2018SZDZX0037], and the Chengdu Giant Panda Breeding Research Foundation [CPF2017-20]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.