Complete Mitochondrial Genome of Turdus merula (Aves: Passeriformes: Turdidae) and Related Species: Genome Characteristics and Phylogenetic Relationships

The fast evolutionary rate, relatively conserved gene content and organization, maternal inheritance, and limited recombination make mitochondrial genomes an useful neutral molecular marker for studies related to species identification, verification taxonomic levels, and identification phylogenetic relationships (Boore and Brown, 1998; Anmarkrud and Lifjeld, 2017; Ingman et al., 2000; Wolstenholme, 1992; Bernt et al., 2013). The new sequencing and PCR technologies have made the utilization of the mitochondrial genome easier and more frequent in recent decades (Powell et al., 2013; Alam et al., 2014; Mu et al., 2015). Since the identification first complete mitochondrial genome of birds (Desjardins and Morais, 1992), nearly 2000 mtDNA sequence of birds have been submitted to NCBI (http://www.ncbi.nlm.nih.gov/). The vertebrate mtDNA is a compact double-stranded closed circular molecule that contains 2 ribosomal RNAs (rRNA), 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNA), and a non-coding control region (CR) (Simon and Frati, 1994). Some of these genes such as Cytb, COXI, and COXII have frequently been used for population genetic studies (Alam et al., 2014; Velez-Zuazo and Agnarsson, 2011).

Nevertheless, genetic information may lack sufficient resolution and may not be always rightly based on these short mitochondrial regions (Alam et al., 2014; Velez-Zuazo and Agnarsson, 2011). Compared with single or partial mitochondrial gene fragments, such as COX I and Cytb, mitochondrial genome sequence can provide more information and faster substitution rate and more insight and better resolution from higher-level groups to closely related species (Powell et al., 2013; Mu et al., 2015).

The monophyly of the order Passeriformes is strongly supported by the morphological characteristics and molecular data (Raikow, 1986; Edwards, 1991; Sibley and Ahlquist, 1992; Zhang et al., 2018). This order is divided into two major clades, the suboscines and oscines. The oscines can be further split into two suborders: Corvida and Passerida. Corvida is further divided into three superfamilies: Menuroidea, Meliphagoidea and Corvoidea. Passerida can be also divided into three superfamilies: Muscicapoidea, Sylvioidea, and Passeroidea (Sibley and Ahlquist, 1992; Zhang et al., 2018). However, the taxonomic status of passerines is rather confusing, since most of the evolutionary groups underwent very rapid radiation during the Paleogene (Feduccia, 1995). In addition, in some studies, the monophyly of Sylvioidea has not been supported (Sheldon and Gill, 1996; Barker et al., 2002; Alström et al., 2006), thus owing to the monophyletic taxonomic relationship of Passeridae and Sylviidae, Paridae should be removed from Sylvioidea. However, Alström et al. (2006) thought that the Alaudidae originally belonged to Passeroidea and should therefore be included in the Sylvioidea. These controversies are at a deep level of the taxonomic status between Muscicapoidea, Sylvioidea, Passeroidea, and Paridae (Barker et al., 2004; Johansson et al., 2008; Nabholz et al., 2010; Zhang et al., 2018).

Turdus merula, the national bird by Sweden, is a passerine bird of the family Turdidae (Passeriformes), a large family which is composed of 341 species in 24 genera (Dickinson, 2003; Sangster et al., 2010). They are widespread in temperate Eurasia, North Africa, the Canary Islands, and South Asia, and have been introduced to Australia and New Zealand (Long, 1981). The T. merula may be a resident, partially migratory, or fully migratory depending on latitude (Zheng, 2018). In China, the T. merula is a resident bird extending from the Yangtze River to the Tianshan Mountains (in Nanchong, Sichuan, and Zipeng mountain in Hefei province). However, in Hainan province, they are winter migratory birds (Luo et al., 2008; Ni, 2014; Zhou et al., 2001). In recent years, although T. merula population is increasing, due to their natural habitat loss and urbanization part of the T. merula has colonized from original forest to urban (Wang and Yin, 2016).Thus the T. merula was listed as less concern species by IUCN in 2021 (http://www.iucnredlist.org/search). To strengthen the protection of T. merula and marginal species of it, many scholers have studied the classification of T. merula and Turdidae, and some achievements have been made. Because of its easy amplification, lack of insertion and deletion, and sufficient variation, the COX1 gene in mtDNA plays an important role in bird identification and phylogeny (Seutin and Bermingham, 1997). Thus, Xu et al. (2010) studied 14 birds of the subfamily Turdidae and explored their phylogenetic relationships with COX1 genes in the year 2010, and the results showed that the family of Turdidae can be divided into 2 large branches, including Turdus, Zoothera; Phoeicurus, Tarsiger, Enicurus, and Myiophonus (Xu et al., 2010). Some phylogenetic relationships strongly suggest that Muscicapidae is not monophyletic or paraphyletic with mitochondrial genes (Zhang et al., 2018) and multilocus study (Sangster et al., 2010; Payevsky, 2014). Furthermore, the relationships between the family Turdidae species are suspected as a monophyletic group (Payevsky, 2014; Zhang et al., 2018). In addition, there is a parallel relationship between generic level in Turdidae or Muscicapidae (Sangster et al., 2010; Zhang et al., 2018).

In this study, we extracted, measured, and analyzed the characterization of the complete mitochondrial genome sequence of T. merula. Moreover we analyzed and compared its complete mitochondrial genome DNA sequence with 78 mitochondrial genome sequences of other Passeriformes to provide insight into their genome evolution as well as the phylogeny. We reconstructed the phylogenetic tree within this order based on 13 PCGs of the selected birds to understand the sequence characteristics and evolutionary status and to infer their higher phylogeny.

A single specimen of T. merula was obtained from Zitong County, Mianyang City, Sichuan Province in China (105°4’14.04”, 31°49’49.95”, 555m above sea level). Tissue samples (2 mL) were obtained and preserved with heparin anticoagulant tube prior to the analysis. The complete mitogenome was extracted and purified from the preserved muscle tissue by the Proteinase K/SDS digest extraction method followed by phenol–chloroform isolation and ethanol precipitation (Sambrock and Russel, 2001). After the quality of the obtained genomic DNA was checked by using 0.9 % agarose gel, it was stored at -20°C until needed for PCR.

The twelve overlapping fragments that we amplified with normal PCR and LA-PCR covered the entire mtgenome of T. merula, and each of them overlapped more than 120-230 bp. Complete mtDNA was amplified as concatenated sequences using selectively amplified mtDNA template and seven primer pairs derived from the literature (Zhao et al., 2012; Jiang et al., 2014). The remaining PCR primers were designed based on the alignments of the relatively conserved nucleotide sequences in 6 homologous Turdidae species in GenBank and designed by using Primer (Premier 5.0 software). Twelve pairs of primers (Table I) used for PCR amplification were used in the reaction volume of 25 μl which contained 2.5μl 10×loading buffer, 2.0μl of MgCl2 (2.5mol/L), 1.5μl dNTP mix (2.5mM/L each), 1.0 µL of each primer (10 μmol/L), 1.0 μl DNA template (20 ng/µL), 0.6μl (5U/µL) of LA Taq polymerase and sterilized distilled water. The thermal cycle comprised an initial denaturation at 94 °C for 2.5 min, 33 cycles each of denaturation at 94 °C for 45s, annealing at 50-61 °C for 40 s, extension at 72 °C for 80-180 s, and a final extension at 72 °C for 9 min. The products of PCR were separated by 1.0 % TAE agarose gel electrophoresis and recovered using a Gel Extract Purification System (Omega bio-tek, U.S.A) and sequenced by using an ABI 3730 sequencer, either directly or following sub-cloning into the pUC19 DNA vector (TaKaRa, Japan).

The sequence assembly and annotation were performed using DNA Baser Assembler and manual screening (http://www.DNABaser.com). The primary DNA sequence data were homologous alignment using BLAST searches at NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Protein-coding genes (PCGs) and rRNA genes were identified by multiple sequence alignments with previouslysubmitted sequences of three Turdus species as well as by secondary structure search (Barker, 2014; Yan et al., 2016; Gibb et al., 2015). Both transfer RNAs genes and their secondary structures of the stem-loop were identified by tRNAscan-SE Search Server v.1.21 (http://lowelab.ucsc.edu/cgi-bin/tRNAscan-SE2.cgi) using the default arguments (Lowe and Eddy, 1997). And the cloverleaf secondary structures of T. merula were painted by using Photoshop. The tandem repeats of noncoding regions were analyzed with the Tandem Repeats Finder program (http://tandem.bu.edu/trf/trf.advanced.submit.html) (Benson, 1999). The base composition of the complete mitochondrial genome, 13 PCGs, 22 tRNAs, 2 rRNAs, D-loop, codon usage of 13 PCGs were analyzed using MEGA 6.06 (Tamura et al., 2013). Composition skew was calculated according to the following formulas: AT-skew = [A–T]/[A + T] and GCskew = [G–C]/[G + C], respectively (Lobry, 1996; Perna and Kocher, 1995). The complete mtDNA sequence of T. merula reported in this article was deposited in GenBank under the accession number MN248536.

Note: The asterisk indicates the species in this study.

To explore the phylogenetic relationships between T. merula and other Aves (Passeriformes: Turdidae, Sturnidae, Muscicapidae and Paridae, and so on), all available complete mitochondrial genomes of them were downloaded at NCBI Genbank. Moreover, two species Menura novaehollandiae (NC_007883), and Climacteris picumnus (KY994598) were set as outgroup (Table II). All 13 protein coding genes were aligned using MEGA 6.06 with default settings (Tamura et al., 2013). Subsequently, all 13 PCG alignments were combined to create a concatenated data set. The optimal nucleotide substitution model was selected using jModeltest v.0.1.1 (Posada, 2008). Maximum Likelihood (with 1000 bootstrap replicates) and Akaike Information Criterion (AIC; Posada and Buckley, 2004) scores (Bayesian Information Criterion) were used for phylogenetic tree construction using PhyML v.3.0 (Guindon and Gascuel, 2003) and the Bayesian inference (BI) analysis was performed with MrBayes v.3.1.2 (Huelsenback and Ronquist, 2001). The best partitioning schemes for each partition were determined using PartitionFinder v2.1.1 (Lanfear et al., 2017). The mitogenome data matrix was divided into 39 partitions based on codon positions for each PCGs (13 genes × 3 codons = 39 partitions). The best fitting model of GTR+I+G (nst=6; rates=gamma) was selected for subsequent Bayesian phylogenetic analyses. Four Markov chains (one cold and three hot chains) were simultaneously run at five million generations, sampling every 1000 generations.

The complete mitogenomes of T. merula was a 16,734 bp in length. The circular DNA molecule (Table III, Fig. 1), which contains 13 protein-coding genes, 2 rRNA genes, 22 tRNA genes, and one non-coding region. Among the 37 fragment genes, 9 genes (including 1 protein-coding genes and 8 tRNA genes) were encoded by L chains, and the rest by H chains. Its overall base composition (H-strand) was: A, 29.34%; C, 32.50%; G, 14.82% and T, 23.34%, so the percentage of A and T (52.67%) was slightly higher than G and C (47.32%), within the range for avian mitogenomes (51.6-55.7%), similar to those of other thrushs (Haring et al., 2001; Rong et al., 2010; Yan et al., 2016; Li et al., 2016; Peng et al., 2016). The mitochondrial genome sequence of T. merula contained 42bp (0.25% of the whole sequence) overlapping nucleotides. Moreover, varied in length from 1 to 10 bp and the largest overlap (10bp) is located between 12S rRNA and tRNAVal, ATP6 and ATP8, respectively. A total of 74bp (0.44% of the whole sequence) intergenic nucleotides (IGN) were dispersed in 18 locations and ranged in size from 1 to 10 bp; the largest gene interval can be found in three places (tRNAAsp and COX2, ND5 and Cytb, tRNAPro and ND6). The details of gene location were given in Table III.

We identified two long gene overlaps in T. merula, namely, 10-bp long overlaps (ATGAACCTAA) and 7-bp long overlaps (ATGCTAA), which are also found in many other Passeriformes species (Fig. 2, Table III). The two overlaps were located between ATP8 and ATP6 and between ND4L and ND4 on the H-strand, respectively. The overlapped sequences between genes were thought to be translated as a bicstron (Stewart and Beckenbach, 2005).

* T represents incomplete stop codons. † Intergenic bp indicates gap nucleotides (positive value) or overlapped nucleotides (negative value) between two adjacent genes. ‡ H and L indicate genes transcribed on the heavy and light strands, respectively.

Another 9 bp long overlapping motif (CAAGAAAGG) was detected between COXI and tRNASer (UCN) in mitogenomes of T. merula, which was also present in other Passeriformes (Fig. 3). The sequenced motifs between ATP8 and ATP6, between ND4L and ND4, and between COXI and tRNASer (UCN) were relatively conserved in the Passeriformes species after mitogenomic comparisons. We found that gene overlap regions are ubiquitous in eukaryotic mitochondria. The existence of gene overlap regions enables limited genomic sequences to encode more genetic information, which is very economical and effective for the transmission of genetic information of species.

All the 13 protein-coding genes found in other animals were also presented in T. merula, including 7 NADH dehydrogenase subunits (ND1-6, ND4L), 3 cytochrome c oxidase subunits (COXI-III), 1 Cytb, 2 ATPase subunits (ATP6, ATP8), and one cytochrome b gene (Cytb). The 13 mitochondrial protein-coding genes were 11,394 bp in length, accounting for 68.09% of the entire mitogenome sequence. The base composition of 13 PCGs were shown in Table IV. The A+T content of the 13 PCGs was 51.67% and the AT skew (0.10) of PCGs was slightly positive, while the GC skew (-0.45) was strongly negative. Contents with A+T in the second position were slightly higher than G+C while those in the first and third positions were on the contrary (Table IV). Two start codons (ATG, GTG) (ATG account for 92.31% of all initial codons), four stop codons of three (TAA, AGG, TAG), and the single incomplete stop codon (T) were used for initiating and terminating the coding of mitochondrial 13 PCGs. Only one protein-coding gene (COXI) utilized GTG as the start codon and all the others used ATG as standard start codon. TAA was the most frequent termination codon and seven protein-coding genes (COXII, ATP8, ATP6, ND3, ND4L, ND5, Cytb) ended with it; whereas ND1 and COXI ended with AGG and ND6 with TAG, which were also found in other Turdidae birds (Li et al., 2016; Zhang et al., 2018). The uncanonical T termination codon was used in the other three PCGs (ND2, COXIII, and ND4), which may be completed by poly-adenylation of the 3’-end of the mRNA after transcription (Boore, 1999; Yoon et al., 2013). In addition, the codon usage was shown in Table V and Figure 4. Encoding 3,797 amino acids (excluding stop codons), the most frequent amino acids in the 13 PCGs of T. merula were Leucine (19.44%), then Isoleucine (11.43%), and the next Alanine (9.15%). The highly abundant like this were similar to mitochondrial proteins of other birds (Liu et al., 2015, 2016; Song et al., 2016; Yong et al., 2015). According to relative synonymous codon usage shown in Table IV and Figure 5, the RSCU values of UUC, CUC, CUA, AUC, and other 26 codons were great than or equal to 1, they were called preference codons of the T. merula mitochondrial genome.

The 12S rRNA was 984bp and the 16S rRNA was 1601bp in length, which were located between tRNAPhe and tRNALeu (UUR) genes, and separated by tRNAVal gene, as other avian rRNA genes (Yoon et al., 2015). The base composition of the two rRNA gene sequences are shown in Table V. The content of A+T (53.4%) was higher than that of C+G (46.6%).

Like most Turdus mtDNA, scattering throughout the mitogenome, the typical set of 22 tRNA genes was identified, ranging from 66 bp (tRNASer(AGY)) to 75 bp (tRNALeu and tRNAAsn ) in size (Table V). Among those tRNAs, fourteen tRNA were encoded on the H-strand and eight on the L-strand. Their secondary clover-leaf structures were predicted by tRNAscan-SE Search Server and presented in Figure 2. Only one (tRNASer(AGY)) of them can not fold into the canonical cloverleaf secondary structure, whose dihydroxyuridine arm is absent. These features are common in most metazoans mitogenomes (Ohtsuki et al., 2002; Yoon et al., 2015). Furthermore, except the tRNASer(AGY), the other tRNA genes each have a 7bp length on the amino acid acceptor arm and the anticodon loop; the length of the DHU arm is 3 to 4bp, the anticodon arm is 3 to 5bp, the TψC arm was 4 to 5bp; the DHU ring, variable ring, and T ring vary greatly. This may be one of the main reasons for the variation of the tRNA length. The sequence and structure of anticodon, amino acceptor, and TψC arms were highly conserved, while the structure of the variable loop was highly diverse with obvious indels polymorphisms (Chen et al., 2018a, b). In addition to the typical A-U and G-C pairing, there were 30 swinging mismatched base pairs (G-U) in the mitochondrial genome of tRNA secondary structure of T. merula, most of which were found in the amino acid acceptor arm (12 places) and the anticodon arm (9 places). In addition, the remaining contained 4 places in DHU arm and 5 places in the TψC arm, and there are also some other mismatched bases, such as U-U, C-C, A-A and A-C (Fig. 6). In most arthropod mtDNA, these mismatches might be corrected by RNA editing, so that they can not lead to obstruction in amino acid transportation (Yokobori and Pääbo, 1995; Dang et al., 2008; Liu et al., 2012). At the same time, Varani and Mcclain (2000) believed that the unmatched base pair G-U can play an important role in the stability of tRNA secondary structure.

The non-coding region of T. merula includ a 1,180 bp long control region (D-loop) and a few intergenic spacers. As is the case in most avian mitogenomes, the D-loop region is located between tRNAGlu and tRNAPhe on H strand. The nucleotide composition of the T. merula CR is 25.5% A, 29.4% T, 31.0% C, and 14.0% G, with a distinct bias against G. The AT-skew is slightly negative (-0.07) while the GC-skew is strongly positive (0.38). Comparing the conserved motifs with those of other avian CRs and according to the basis of the variability (Anderson et al., 1981; Brown et al., 1986; Doda et al., 1981), T. merula CR could be divided into three domains (Fig. 7): (i) the extended termination-associated sequence motif (ETAS) domain, which is associated with the termination of the newly synthesized H-strand during replication (Ryu and Hwang, 2012; Sbisa et al., 1997; Yoon et al., 2015); (ii) the central domain (CD); (iii) the conserved sequence block (CSB) domain. Locating between the 5’-end of the CR and the beginning of the F-box in the central domain, the ETAS domain contains conserved sequence blocks (CSB1-like) and two Extended termination associated motifs (ETAS1-2). Futhermroe closer to the 5’-end of the CR, a cytosine stretches sequence (CTCCACCCCCCCCCCTTCCCCCCC) was found, which also exists in many bird species (Randi and Lucchini, 1998; Ritchie and Lambert, 2000; Yoon et al., 2015). The central domain contains six highly conserved sequence blocks (F-box, E-box, D-box, C-box, b-box and B-box). The F-box is at the upstream of the CD region while the B-box at downstream. CSB region is known to contain three conserved regions (CSB1, CSB2, CSB3) in the mitochondrial of vertebrates (Walberg and Clayton, 1981), but there is only one CSB1 and a bi-directional transcriptional promoter (HSP/LSP box) in T. merula of this region. Gene duplications, rearrangements and missing may lead to gene order changes. Generally, there are two types of arrangements between Cytb and 12S rRNA genes of mitochondrial organisation: (i) tRNAThr/tRNAPro/NADH6/tRNAGlu/CR/tRNAPhe; (ii) tRNAThr/ CR/tRNAPro/NADH6/tRNAGlu/ NC/tRNAPhe and the T. merula follow the first pattern.

The bayesian inference (BI) and maximum likelihood (ML) phylogenetic trees, which are based on the concatenated 13 PCGs of mitochondrial gene sequences of 78 passeriformes species and two outgroup species (Menura novaehollandiae and Climacteris picumnus), were reconstructed. In our study, the concatenated PCG data set of the mitogenome sequences included 11,382 nucleotide positions, including 5,229 conserved sites, 6,153 variable sites, and 5,577 potentially parsimony-informative sites. The phylogenetic relationship among these species shows the tree topologies of ML and BI were similar, except for several Turdidae species (T. migratorius, T. rufiventris and T. merula) in relationships or placements (Fig. 8).

Our phylogenetic trees strongly support the monophyly of Passeriformes group, as suggested by a previous study (Sibley and Ahlquist, 1992). The topological structure of the phylogeny of passerines allowed Sibley and Ahlquist (1992) to verify the traditional division of passerines into two suborders, the Suboscines and Oscines forming two obvious monophyletic lineages. The Oscines can be further split into two suborders: Corvida and Passerida, and Passerida can be divided into three superfamilies: the Muscicapoidea, the Sylvioidea, and the Passeroidea (Sibley and Ahlquist, 1992; Zhang et al., 2018). Although Passeroidea is a sister taxon of Muscicapoidea and Sylvioidea, Muscicapoidea and Sylvioidea are more closely related to each other. This relationship is also inconsistent in previous studies (Barker et al., 2002, 2004; Marshall et al., 2013; Zhang et al., 2018). Furthermore, corvida and passerida cluster on a large branch, and form a sister classification relation. Terpsiphone atrocaudata (Monarchidae) is located at the base of the Corvidae and Laniidae clade, while Regulus calendula and R. regulus (Regulidae) combine into a clade at the base of the Passerida group (Fig. 8). The two species were previously placed in the family Sylviidae, which is inconsistent with the results of our phylogenetic tree, and it has been previously reported (Keith, 2014). The systematic evolutionary relationships of two species representing the Regulidae in Passerida belong to other taxonomic groups because of the high support value of nodes (1.00/975) (Fig. 8). However based on the phylogeny reasoning of our study, there may be a clear statement when the sampled taxa and molecular markers are added. Hyliota flavigaster (Hyliotidae) forms a sister group with the family Paridae at the base of this family. One obvious result of the Fuchs et al. (2006) and Johansson et al. (2008) studies was the identification of the presumptive sylvioid genus Hyliota as another deeply-diverging member of the Passerida with unclear affinities. In line with our phylogenetic tree, H. flavigaster has a distant relationship with the family Sylviidae, which is different to the previous study and considered it to family Sylviidae (Keith, 2014). Our data suggest that it should be classified into the superfamily Sylvioidea.

The phylogenetic analysis based on the 13 PCGs of the mitochondrial genome suggests that the species of genus Paradoxornis should not be placed in the superfamily Muscicapoidea or family Muscicapidae (Fig. 8). Based on our phylogeny, we propose that genus Paradoxornis should be placed in the superfamily Sylvioidea or the family Sylviidae. Moreover, the systematic taxonomic status of the P. humilis (Paridae) has long been controversial. The species is classified into the family Corvidae based on skeletal characteristics and taxonomic data, while Hope (1989) considered P. humilis not to be placed in Corvidae. Then, James et al. (2003) obtained similar results with different methods, namely, the skeletal characteristics, nuclear gene (c-myc), and mitochondrial gene (Cytb). Subsequently, some researchers (Ericson and Johansson, 2003; Johansson et al., 2008; Treplin et al., 2008) strongly supported that this family should be separated from all other Sylvioidea, and independently form a taxonomic group, which is consistent with previous studies (Alström et al. 2006). In our mitochondrial genome analysis, the P. humilis is nested within other species of the family Paridae (Fig. 8), and should belong to this Sylvioidea. However, this relationship is inconsistent with previous research and proposed that it should be placed within the Muscicapoidea (Ericson et al., 2003; Johansson et al., 2008; Treplin et al., 2008; Zhang et al., 2018). Furthermore Marshall et al. (2013) suggested that the family Paridae represented by P. humilis should be separated from Sylvioidea and moved to the superfamily Muscicapoidea. However our phylogenetic data, similar to Zhang et al. (2018), show that this family belongs to Sylvioidea and confirmed that P. humilis was included in the Paridae.

Within the taxa of the superfamily Muscicapoidea, Sangster et al. (2010) and Zhang et al. (2018) strongly suggest that Muscicapidae is not monophyletic, and exists extensive paraphyly at the family and genus levels within this group. If some species are not classified into other families, such as Paradoxornis, Luscinia cyanura, and Monticola gularis, our taxonomic relationships are similar to previous studies (Sangster et al., 2010; Zhang et al., 2018). Based on the topology tree we reconstructed, we suggest that the genus Paradoxornis needs to be classified as Sylviidae, while the two species (Luscinia cyanura and Monticola gularis) are placed in the family Muscicapidae (the specific classification is shown in Fig. 8). The formation of the sister taxa of these two families, Muscicapidae and Turdidae, indicates that they are closely related, and then they form the superfamily Muscicapoidea together with the Sturnidae families (Fig. 8).

In the family Turdidae, M. myadestinus was basal position of the other Turdidae species, indicating that its differentiation is relatively primitive. And T. merula form a big branch with T. obscurus, T. celaenops, T. eunomus, T. kesseri, T. mupinensis, T. naumanni, T. ruficollis, T. dissimilis, T. cardis and T. hortulorum, which indicates that they are closely related (Fig. 8). The species evolutionary relationship of genus Turdus we constructed is consistent with that of Xu et al. (2010), but different from that of Peng et al. (2016) (T. merula, T. migratorius; T. naumanni, T. hortulorum), Batista et al. (2020) (T. migratorius; T. rufiventris; T. eunomus; T. naumanni; T. obscurus; T. hortulorum; T. merula), Nylander et al. (2008) (T. naumanni; T. ruficollis; T. kesseri; T. obscurus; T. celaenops; T. cardis; T. hortulorum; T. dissimilis; T. migratorius; T. rufiventris; T. merula; T. abyssinicus; T. mupinensis) and Nagy et al. (2019) (T. naumanni; T. ruficollis; T. obscurus; T. kesseri; T. cardis; T. hortulorum; T. dissimilis; T. merula; T. rufiventris; T. migratorius). However in T. merula, T. migratorius, T. mupinensis and T. rufiventris there are some differences in the topological structure among different researchers. Hence, we suggest that the specific systematic taxonomic relationships of these species should be further analyzed and discussed. It may be possible that such a taxonomic status could be recovered using a more comprehensive taxa sampling of complete mtDNA genomes and supplementary nuclear gene molecular markers.

With regard to the New Zealand wattlebirds (Philesturnus carunculatus and Callaeas cinereus), our phylogenetic trees are strongly supported as a monophyletic group with both nuclear gene and mitochondrial DNA sequence data (Shepherd and Lambert, 2007; Gibb et al., 2015). Furthermore, the position of the New Zealand wattlebird taxa within the Oscines, as determined in our study (i.e. embodied in the oscines but excluded from the Passerida and Corvoidea), was consistent with the placement of P. carunculatus in previous studies (Barker et al., 2004; Shepherd and Lambert, 2007; Gibb et al., 2015). However, in our system topology, there are different classification arrangements (Fig. 8). For example, P. carunculatus and C. cinereus congregate in the same branch and are closely related to the Passerida, and they are located at the base of the order. Therefore, it is suggested that the two species should be classified into the order Passerida.

In this study, the complete mitochondrial genome of T. merula was determined and analyzed, and it is similar to other Passeriformes with many significant features. The motifs obtained by sequence comparison, “ATGAACCTAA” between ATP8 and ATP6, and “ATGCTAA” between ND4L and ND4, and “CAAGAAAGG” between COXI and tRNASer(UCN) were highly conservative in Passeriformes species. In the current research, the phylogenetic relationships based on nucleotide sequences of 13 PCGs showed that the phylogeny analyses of Passerida were conducted with the clear support of dividing the group into three superfamilies: the Muscicapoidea, the Sylvioidea, and the Passeroidea, and Passeroidea is a sister taxon for Muscicapoidea and Sylvioidea. We suggest that the genus Paradoxornis will be classified as Sylviidae, while these two species (Luscinia cyanura and Monticola gularis) are placed in the family Muscicapidae. Furthermore, Turdidae formed a sister group with Muscicapidae, and then they form the superfamily Muscicapoidea together with the Sturnidae families. The relationships between some species of the order Passeriformes may remain difficult to resolve despite an effort to collect additional characters and samples for phylogenetic analysis and our results could be useful in future research on population genetic structure, phylogeny, and conservation genetics.

We would like to express our gratitude to all those who helped us during the writing of this manuscript.

The animal handling procedures conformed with the China Animal Welfare guidelines and have been approved by the Ethic and Animal Welfare Committee and the Animal Protection and Use Commission of Mianyang Normal University (MNU: 19ZBJC-01).

This research was supported by the Research Project of Ecological Security and Protection Key Laboratory of Sichuan Province (No. ESP2003), the Research Project of Education Office Project of Sichuan Province (No. 18ZA0261), the Scientific Research Fund of Mianyang Teacher’s College (Nos. PY-2016-A03, MYSY2017JC02, Mnu-JY20035 and Mnu-19ZBJC-01).

This study was carried out in accordance with the guidelines of Animal Care and Use Committee at the Mianyang Normal University. Efforts were taken to minimize bird suffering by administering anesthesia. The study did not involve endangered or protected species.

The authors have declared no conflict of interest.