Characterization of the Complete Mitochondrial Genome of Sea Duck Mergus serrator and Comparison with other Anseriformes Species

Mitochondrial genome (mitogenome) harbors useful molecular information that can be used to infer phylogenetic relationships among Anseriformes different taxa (Livezey, 1997; Donne-Goussé et al., 2002; Jiang et al., 2009; Sun et al., 2017), one of the best-studied groups of birds (Livezey, 1986). Anseriformes includes 3 families, 8 tribes, and more than 150 species (Johnsgard, 2010; Ringelman et al., 2015). Due to gene rearrangement, genetic introgression, gene flow, and hybridization, the taxonomic statuses and phylogenetic placements within the order remain in dispute (Scribner et al., 2001; Gonzalez et al., 2009; Lavretsky et al., 2021). The Mergini tribe includes 20 species of which two are extinct (Johnsgard, 1966). As a key group in the Anatidae, the Mergini tribe (sea ducks) is strongly structured, with clear genetic assignments and phylogenetic relationships (Johnsgard, 2010). Whereas, the relationships among the Mergini tribe species are not certain based on multiple approaches and different genes so far. For example, the phylogenetic status of Mergus serrator is inconsistent based on different genes (Liu et al., 2012, 2014; Lavretsky et al., 2021). There is a remarkable difference in the life history and behavioral traits among different tribes in the Anseriformes (Peters et al., 2012). Most tribes tend to live in inland streams, rivers, or ponds, such as the Anatini, Aythyini, and Anserini tribes (Debela et al., 2021; Ponomarenko et al., 2021). In contrast, the Mergini tribe roosts and forages in sea habitats (Johnsgard, 2010). The tribes also differ in terms of breeding strategy, longevity, and bionomic strategy (Savard et al., 1998; Johnsgard, 2010). The Mergini tribe as an excellent tribe to investigate molecular adaptions to elevate metabolic burden associated with other life-history characteristics. Phylogenetic analyses can help investigate the relationships between mitochondrial DNA (mtDNA) molecular evolution and metabolic performance (Sruoga et al., 2008; Sonsthagen al., 2011).

MtDNA is typically double-stranded and maternally (Wolstenholme, 1992; Roques et al., 2004), it is valuable for studies of avian taxonomy and phylogeny (Shen et al., 2009; Jetz et al., 2012; Ren et al., 2014), genetic structure (Ruan et al., 2018), biological identification (Lecroy and Barker, 2006; Zhu et al., 2021), and conservation genetics (Allendorf et al., 2010), as it is a neutral molecular marker (Torroni et al., 2001; Nabholz et al., 2008). As the main energy suppliers of the cell, mitochondria produce approximately 90% of cellular ATP reserves through the highly efficient chemical osmotic coupling of electron and proton transfer to ATP synthesis (Saraste, 1999; Starkov, 2008). Evaluating the selective pressure of environmental temperature and oxygen availability on changes in mtDNA molecular could provide key insights for the adaptive evolution of mitogenomes (Luo et al., 2013). The key to successful migration to a new habitat is whether it can adapt to different energetic demands. Energetic demand is essentially related to the ability of mitochondrial to generate energy through the process of oxidative phosphorylation, such as in certain positions encoded by the NADH dehydrogenase genes (ND4 and ND5) (Janssen et al., 2004; Liu et al., 2018). However, single or combined genes have obvious shortcomings in phylogenetic analysis, such as insufficient information (Liu et al., 2012). These are overcome when using complete sequences of mitochondria, and it is becoming the first-choice marker for resolving controversial species relationships (García et al., 2014; Peng et al., 2017).

In this study, we sequenced the mitogenome of a typical sea duck, M. serrator, to analyze the genetic structure of and phylogenetic relationships among the Mergini tribe. In recent years, an increasing number of mitogenomes of the Mergini tribe have been sequenced. This study aims to shed increased light on the phylogenetic status of M. serrator in the Mergini tribe and to elucidate the evolutionary rates and molecular signatures of natural selection for all 13 mtDNA protein coding genes (PCGs) in sea ducks.

The sample collection was strictly conducted under national ethical guidelines (Regulations for Administration of Affairs Concerning Experimental Animals, China, 1988) for animal husbandry and humane treatment.

We collected a dead specimen of M. serrator from Dalian (38°51′52.22″N 121°39′39.11″E; Liaoning, China). A muscle sample was transferred to a 5 mL centrifuge tube, 1 mL 75% ethyl alcohol was added, and the sample was stored in a refrigerator at –20°C. The sample (specimen voucher: NIES20210320HXQSYS01) was deposited in the Museum of Laboratory of Biodiversity Investigation of Nanjing Institute of Environment Science.

Total genomic DNA was extracted from muscle tissue using the Easy Pure genomic DNA kit (Trans Gen Biotech Co, Beijing, China). The original sequence data (15.65 G) was deposited in NCBI’s Sequence Read Archive (SRA; accession: SRR13516389), and complete mitogenome sequencing (15.62 G) was performed on an Illumina NovaSeq6000 platform (Novogene Bioinformatics Technology Co. Ltd., Tianjin, China) (Tang et al., 2021). The mitogenome was assembled using the program NOVO Plasty 3.7 (Dierckxsens et al., 2017) and then adjusted manually. The annotated using the MITOS Web Server (Bernt et al., 2013). The assembled mitogenome sequences have been deposited in GenBank under accession numbers MZ365040.

All sequences were aligned on Clustal X using the default options (Thompson et al., 1997). PCGs were identified by comparison with the corresponding known complete mtDNA sequence of M. merganser and M. squamatus using Sequin 11.0. The 22 tRNA genes of cloverleaf secondary structures and anticodon sequences were identified from MITOS WebServer (http://mitos.bioinf.uni-leipzig.de/index.py). The control region (CR) was identified by sequence homology analyses (Cadahia et al., 2009; Barker et al., 2012). The start and stop codons of 13 PCGs were identified in the M. serrator mitogenome by ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/) and use the CGView Server software to generate a graphical map of the mitogenome and modified it manually (Table I). The codon number and relative synonymous codon usage (RSCU) of three species of Mergus mitochondrial PCGs were obtained using MEGA X and PhyloSuite.

Negative values indicate overlapping; H, heavy strand; L, light strand.

Retrieved and downloaded the complete mitogenomes of 53 Anseriformes species from GenBank, and including 2 families and 46 genera (Table II). We used MEGA5.0 (Tamura et al., 2011) and DNASTAR (Burland, 2000) to compare the sequences of these birds and to analyze their base content. The calculation for the skewing of nucleotide composition was AT skew = (A-T) / (A + T), GC skew = (G-C) / (G + C) (Perna and Kocher, 1995).

We collected the complete mitogenome of 53 Anseriformes species for phylogenetic analysis (Table II). Corresponding Gallus gallus, Hydrophasianus chirurgus, and Porzana fusca sequence were used as outgroups. Phylogenetic trees were estimated using ML and BI methods. The optimal parameter model selected by MrModel Test 3.06 (Nylander et al., 2004) and PAUP 4.0b10 software (Swofford, 2003). The GTR + I + G model that was selected as the best fit model for nucleotide phylogenetic analysis. The BI trees of Anseriformes birds were constructed using MrBayes 3.1.2 software (Ronquist et al., 2012). The ML tree was constructed using the MEGA5.0 software (Tamura et al., 2011), as follows: run the files in Phy format, select the best model using the Mr Model Test 3.06 software filter, set the ML in bootstrap mode, repeat bootstraps 50 times, and run the whole process 10,000 times to obtain the ML tree. Finally, Treeview32 (Saldanha, 2004) software was used to create and annotate the evolutionary tree diagram.

In order to test the possible influence of topology on the inference of selected locations, estimated ω (dN/dS) rate using the topology inferred from genomic data (Nadeau et al., 2007; Stoletzki and Eyre-Walker, 2011). dN is the non-synonymous substitution rate, or the rate at which changes in nucleotide sites lead to changes in new amino acid chains. Things are the opposite for dS (synonymous substitution). The dN/dS ratio provides evidence for selective restraint and pressures acting on PCGs (Shen et al., 2009; Botero-Castro et al., 2017). If natural selection is favoring amino acid changes, ω > 1 indicates positive selection. A ratio below 1 indicates purification selection, and a ratio at 1 indicates neutral evolution.

The phylogenetic tree has many branches, and test the role of selection on the Mergini branch. We employed CodeML branch models that estimate the selective pressure between Mergini tribes and other tribes, which allow the estimation of dN/dS of specific branches and clades of interest related to other clades and other parts of the tree (Zhang et al., 2005). In the dual scale and all branch site models described below, we label these species in the Mergini tribes as the foreground branches in a separate analysis. Use branch-site analysis to determine whether a portion of sites was subject to positive selection. To further evaluate these sites, we performed a positive selection analysis using the site model in PAML4. Selection Test was carried out using likelihood-ratio tests (LRT). The model was compared with a null model, in which the ω of the foreground partition class of the selected location was restricted to 1 (Zhang et al., 2005). LRT was used to identify the significance of selective pressure between different pedigrees.

M. serrator is a typical sea duck in the Mergini tribe and its complete mitogenome was found to be 16,627 bp (Table I). So far, the length of the mitogenome of birds that has been sequenced is from 16,300 bp to 23,500 bp (Gao et al., 2009; Lei, 2015). Among Anseriformes species, the shortest mitogenome belongs to Oxyura jamaicensis (15,914 bp) and the longest is that of Anseranas semipalmata (16,868 bp). In the Anseriformes taxa, the length and gene arrangement of the mitogenome are quite conserved. The genomes of M. serrator contained a typical set of 37 genes (Fig. 1), including 13 PCGs, 2 ribosomal RNAs, 22 transfer RNAs, and a non-coding putative CR (Table I, Fig. 2). The structure and gene arrangement were basically the same as most birds (Gao et al., 2009; Lei, 2015). An ND6 gene and eight tRNA genes (tRNAGln, tRNAAla, tRNAAsn, tRNACys, tRNATyr, tRNASer (UCN), tRNAPro, and tRNAGlu) were encoded by the L-strand, whereas other genes were encoded by the H-strand (Table I).

The overall information about the mitogenome of M. serrator was as follows: T: 21.72%, C: 33.38%, A: 28.77%, and G: 16.13% and the content of GC is 49.51% with A+T-rich feature, same as in other species of Anseriformes. The mean base composition of Anseriformes is T (22.31 ± 0.51%), C (32.63 ± 0.64%), A (29.36 ± 0.64%), and G (15.71 ± 0.52%). Additionally, in our study, the sequence of base richness was mainly C > A > T > G, and only Chauna torquata had base contents of A > C >T > G. The AT content range is 50.27% ~ 55.31%, with an average value is 51.67 ± 1.10%, which is higher than the GC content and similar to that of birds in general (50.5% to 57.7%) (Gao et al., 2009).

Vertebrate mtDNA usually shows significant chain bias in nucleotide composition, and this chain bias can be measured as AT and GC skew. We measured the AT skew (0.14) and GC skew (-0.35) of the H-strand in the M. serrator mitogenome. The mean AT skew was 0.14 ± 0.01, and the GC skew was -0.35 ± 0.01. It shows that the nucleotide composition of the complete mitogenome of Anseriformes has a slight specific bias against A and C.

The geneme of M. serrator had 13 PCGs, and 12 PCGs were encoded on the H-strand except ND6. The overall length of 13 PCGs was 11,229 bp, which accounted for 67.53% of the whole mitogenome (Table I). The longest was the ND5 gene (1,824 bp) and the shortest was the ATP8 gene (168 bp). Among the 13 PCGs, the start codon of 9 PCGs was ATG, and the nonstandard start codon GTG was found in COI, COII, and ND5 genes. In addition, the start codon of ND1 gene was ATA (Table I). The most common stop codon was TAA, which has been used six times. TAG and AGG were used twice as stop codons, and the terminal T probably served as the stop signal in the ND2, ND4, and COIII genes. In the ND3 gene, the base A of the stop codon was not translated.

The RSCU values of 13 PCGs are shown in Figure 3. In addition, the start and stop codons, the mitogenome consists of 3662 codons. The most commonly used amino acid was leucine (13.63%) in the use of 13 PCGs of M. serrator, followed by proline (11.47%), serine (9.59%), and threonine (8.06%). In three Mergini species, eight amino acids (Pro, Thr, Leu1, Arg, Ala, Ser1, Val, and Gly) were encoded by four different codons, and the others were encoded by two different codons. The overall AT-skew and GC-skew of the 13 PCGs was 0.06 and -0.36, respectively. The base composition of PCGs has a higher A + C contents, and a preferred to use A or C nucleotides in the first codon position (Fig. 3). This phenomenon also occurs in other species of Mergini, especially M. squamatus and M. merganser.

In the M. serrator mitogenome, ribosomal RNA genes include 12S rRNA and 16S rRNA, which are 995 bp and 1600 bp in length, respectively. They are located between the tRNAPhe and tRNALeu genes, and separated by the tRNAVal, a feature often found in the mitogenomes of birds (Cheng et al., 2011; Jin et al., 2012). Two ribosomal RNAs were identified on the H-strand. The A + T content of both rRNA genes were collectively 52.5% and they had a positive AT-skew (0.28).

A total of 22 tRNAs were found in the mitogenome, and they were scattered between rRNAs and PCGs. Their sizes is between 64 bp and 76 bp (Table I, Supplementary Fig. 1). Three tRNA clusters, namely, IQM (tRNAIle tRNAGln tRNAMet), WANCY (tRNATrp tRNAAla tRNAAsn tRNACys tRNATyr), and HSL (tRNAHis tRNASer (AGY) tRNALeu (CUN)), were determined. Sequences of the tRNA genes could be folded into a canonical cloverleaf secondary structure.

The CR of mtDNAs of M. serrator was 1,067 bp, which was located between the tRNAGlu and tRNAPhe genes. The CR could be divided into three domains according to the base composition (5′-peripheral domain, Domain I; central conserved domain, Domain II; and 3′-peripheral domains, Domain III). In Domain I, extended termination-associated sequences 1 (ETAS1) was observed at position nt 77 - 98 (Fig. 2). A sequence block was also found in domain I (nt 298–323), which is similar to the conserved sequence block (CSB1-like). Four conserved sequence boxes (C, D, E, and F) were identified in Domain II. CBS1 (nt 407–589) and the H- and L-strand transcription promoter sites (nt 806–822) occurred in Domain III.

The complete sequences should be suitable for the phylogenetic analysis and research of Anseriformes. Reconstruction of BI and ML trees using complete mitogenome sequences of 56 species. Use the best-fit model of GTR+G+I (-lnL = 258460.5951, AIC= 516937.2376) and calculated with the Modeltest and Mrmodeltest programs for analysis. Since the generated BI and ML trees show the same topological structure, only the BI tree is displayed (the numbers on the branches represent the bootstrap support for BI and ML trees).

Traditional morphological analysis divides Anatinae into Tadornini, Tachyerini, Cairinini, Anatini, Aythyini, Mergini, Oxyurini, and others. In the literature, Anseriformes species are relatively clearly classified at the level of the family but not the subfamily and genus (Liu et al., 2012; García et al., 2014; Hu et al., 2017). The replacement rate and mutation rate of the mitochondrial sequence of Anseriformes species are high (Fain et al., 2007; Yang et al., 2010). The BI and ML trees showed the complete mitogenome has the ability to distinguish 53 Anseriformes species. The phylogenetic analysis based on the complete mitogenome strongly supports the monophyleticity of the order Anseriformes. The analyses indicated that the Anseriformes include the families Anseranatidae, Anhimidae, and Anatidae. Anseranas semipalmata was included in the branch of Anseranatidae, and Chauna torquata was included in the branch of Anhimidae, and other species belong to the branch of Anatidae (Fig. 4). The branch of Anatidae included Dendrocygninae, Oxyurinae, Stictonettinae, Anserinae, Tadorninae, and Anatinae. The Anatinae had the largest number (32) of species. The Anatinae included Aythyini, Anatini, Somaterini, and Mergini. The results suggest that M. serrator is sister to M. merganser and M. squamatus, which is slightly different from a previous study (Liu et al., 2014). However, both studies suggest that this group belongs to the Mergini.

Most mitochondrial genes besides the ND6 gene are conserved and evolved under purifying selection (Fig. 5). The branch model detected significant signals of positive selection for the ND6 gene tested in all tribes of Anseriformes. It also fitted better for most PCGs when the Mergini tribe was labeled as a foreground branch (Fig. 5) (LRT p value < 0.05). In addition, the ND6 gene showed the highest ω values (ω > 1) both for the foreground (ω = 4.59) and background (ω = 3.28) branches, compared to other genes. The ω value of the ND3 gene in the Mergini tribe is lower than those for other tribes. The ω values of other genes in the Mergini tribe are higher than those for other tribes. The free-ratio model fitting the data of 13 PCGs is significantly better than the null hypothesis (one-ratio model), indicating that there are different evolution rates among the tribes included in our data set.

The mitogenome (16,627 bp) was larger than those of most Anseriformes species. A positive AT skew consistent with other Anseriformes species was found, which implies a slightly specific bias towards A and C in the Anseriformes. PCGs accounted for 67.53% of the whole mitogenome and most of these genes were encoded by the L-strand. ATG and TAA were the most frequent start codon and stop codon, respectively. We summarized the RSCU values for the 13 PCGs and 3788 codons and high A+C contents were found, which also occurred in other species of Mergini. Two ribosomal RNA genes, 12S rRNA and 16S rRNA, were encoded on the H-strand and had a positive A-T skew. Sequences of 22 tRNA genes ranged from 64 bp to 76 bp and could be folded into a canonical cloverleaf secondary structure. The phylogenetic relationships were analyzed and M. serrator was sister to M. merganser and M. squamatus, and formed a closely evolved Mergini clade. Compared to other species of Anseriformes, most PCGs evolved under purifying selection except for the ND6 gene. Different evolutionary rates between the Mergini tribe and other tribes were found.

We thank researcher Yu Dandan for her help in writing this article. This study was supported by the Basic Research Business Project of National Scientific Institute (GYZX210405) and the National Science Foundation of China under Grant (32001111).

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/20220314030351

The authors have declared no conflict of interest.