New Insight into Parrots' Mitogenomes Indicates That Their Ancestor Contained a Duplicated Region.

Mitochondrial genomes of vertebrates are generally thought to evolve under strong selection for size reduction and gene order conservation. Therefore, a growing number of mitogenomes with duplicated regions changes our view on the genome evolution. Among Aves, order Psittaciformes (parrots) is especially noteworthy because of its large morphological, ecological, and taxonomical diversity, which offers an opportunity to study genome evolution in various aspects. Former analyses showed that tandem duplications comprising the control region with adjacent genes are restricted to several lineages in which the duplication occurred independently. However, using an appropriate polymerase chain reaction strategy, we demonstrate that early diverged parrot groups contain mitogenomes with the duplicated region. These findings together with mapping duplication data from other mitogenomes onto parrot phylogeny indicate that the duplication was an ancestral state for Psittaciformes. The state was inherited by main parrot groups and was lost several times in some lineages. The duplicated regions were subjected to concerted evolution with a frequency higher than the rate of speciation. The duplicated control regions may provide a selective advantage due to a more efficient initiation of replication or transcription and a larger number of replicating genomes per organelle, which may lead to a more effective energy production by mitochondria. The mitogenomic duplications were associated with phenotypic features and parrots with the duplicated region can live longer, show larger body mass as well as predispositions to a more active flight. The results have wider implications on the presence of duplications and their evolution in mitogenomes of other avian groups.

Introduction

Vertebrate mitochondrial genomes have usually been regarded as highly conserved because of the stable gene content and order (Clayton 1991; Boore 1999; Lavrov 2007) (fig. 1). The most frequently reported changes in the mitogenome organization are associated with tandem duplications comprising the control region (CR) with neighboring genes and subsequent degenerations of the duplicated copies. Such rearrangements were found in representatives of fishes (Lee and Kocher 1995; Lee et al. 2001; Inoue et al. 2003; Zhuang and Cheng 2010; Shi et al. 2014), amphibians (Mueller and Boore 2005; Kurabayashi et al. 2008; Kurabayashi and Sumida 2013), and reptiles (Kumazawa et al. 1996, 1998; Quinn and Mindell 1996; Kumazawa and Endo 2004; Fujita et al. 2007; Jiang et al. 2007).

. 1.

Comparison of mitochondrial gene orders. The region between gene ND5 and gene 12S rRNA is shown for: typical vertebrate gene order (A), typical avian gene order (B), hypothetical gene order ancestral for Psittaciformes (C) assuming tandem duplication of the cytb to CR segment, and gene orders with duplicated CRs found for parrot species (D). ND5, NADH dehydrogenase subunit 5 gene; cytb, cytochrome b gene; T, tRNA gene for threonine; P, tRNA gene for proline; ND6, NADH dehydrogenase subunit 6; E, tRNA gene for glutamic acid; F, tRNA gene for phenylalanine; 12S, 12S rRNA gene. Pseudogenes are marked by ψ and colored correspondingly to their functional gene copy. The solid blue line above the ND6 pseudogene in GO-1 gene order indicates a fragment identical with ND6 gene; the dashed line above this pseudogene and the tRNA-Glu gene indicates substantial degeneration of their sequences; the solid red line above this pseudogene with the letter “i” indicates an insertion within the pseudogene sequence.

Birds, however, are particularly noteworthy, because their mitochondrial genomes are characterized by a gene order (fig. 1) different from that in the majority of vertebrate mitogenomes (fig. 1). This deviation was initially found in the mitochondrial genome from Gallus gallus (Desjardins and Morais 1990), then it was recognized in many other birds. Therefore, this gene order is regarded as typical for Aves. As in other vertebrates, the fragment encompassing the CR with adjacent genes is also subjected to duplications and pseudogenization, which resulted in various genome rearrangements in many avian orders (Morgan-Richards et al. 2008; Sammler et al. 2011; Wang et al. 2015; Eberhard and Wright 2016; Huan et al. 2016; Akiyama et al. 2017; Hanna et al. 2017; Rodrigues et al. 2017).

Among bird orders, parrots (Psittaciformes) are especially interesting. This group comprises 362 species classified into 86 genera (Clements et al. 2017), which makes this group particularly rich in taxa. Parrot phylogeny is intensively studied, and many classification systems based on genetic as well as morphological data have been proposed (Joseph et al. 2012; Schodde et al. 2013). Therefore, the study of parrot mitochondrial genome evolution, including arrangements, seems to be notably attractive because it can be placed in a wider phylogenetic and taxonomic context.

Duplicated regions in parrot mitogenome were at first found in several Amazona species (Eberhard et al. 2001) (GO-4 in fig. 1). The same rearrangement was further identified in the closely related genus Pionus (Urantowka et al. 2017a). A new insight into the duplication in Psittaciformes was provided by the discovery of a new rearrangement in Prioniturus luconensis by Eberhard and Wright (2016). The presence of two CRs as well as the degenerated copies of genes for tRNA-Thr, tRNA-Pro, tRNA-Glu, and protein ND6 strongly suggested that the Prioniturus rearrangement could be the most similar to a hypothetically ancestral state after a duplication event. Moreover, the similarity between the duplicated region in P. luconensis and the fully duplicated region in other birds indicated that the gene order tRNA-Thr/tRNA-Pro/ND6/tRNA-Glu/CR was the original one that was subjected to duplication in parrot mitogenomes.

The analysis of 20 mitochondrial genomes from parrots showed that such duplications occurred independently at least six times across order Psittaciformes (Eberhard and Wright 2016). Moreover, the results suggested that the common ancestor of parrots had a single version of this region. It accords with the results obtained by Schirtzinger et al. (2012), who did not find any duplication in superfamilies Strigopoidea and Cacatuoidea, which are regarded as early diverged groups in parrot phylogeny (Wright et al. 2008; Joseph et al. 2012). However, it cannot be ruled out that the polymerase chain reaction (PCR) strategy missed identical repeats of cytb-CR or tRNA-Thr-CR regions because the methodology was designed based on the Amazona gene order, which was the only duplicated rearrangement known for parrots at that time. For this reason, we used a different PCR strategy to look for potentially omitted duplications within representatives of Strigopoidea and Cacatuoidea. The new data helped to elucidate the evolution of the parrot mitogenome in terms of duplications. As mitochondrial duplication can be associated with many phenomena, for example, greater longevity (Skujina et al. 2016) and stronger compositional asymmetry between DNA strands (Eberhard and Wright 2016), we also studied relationships of the mitogenomic duplications with the DNA asymmetry and life span as well as basic morphological characteristics of parrots. We also discussed our findings in the context of the mitogenomic duplications in other avian groups.

Results and Discussion

Gene Orders in Complete Mitochondrial Genomes of Cacatuoidea and Strigopoidea

Using an appropriate PCR strategy (Gibb et al. 2007) shown in figure 2, we obtained full mitogenomic sequences for seven taxa: Nestor notabilis, Nymphicus hollandicus, Calyptorhynchus baudinii, Cacatua moluccensis, Cac. pastinator, Probosciger aterrimus goliath, and Eolophus roseicapilla. Interestingly, four of them (N. notabilis, Ny. hollandicus, C. baudinii, and Cac. pastinator) have genomes with the fully duplicated gene order (GO-FD in fig. 1). The presence of identical copies of tRNA-Thr, tRNA-Pro, ND6, and tRNA-Glu genes (supplementary fig. S1, Supplementary Material online) enabled us to precisely position the second CR. In the case of N. notabilis, Ny. hollandicus, and Cac. pastinator, the 5′ ends of two CRs directly follow the 3′ ends of corresponding tRNA-Glu copies. The number of nucleotides between tRNA-Glu and poly-C sequence located at the 5′ end of all known parrot CRs varies from 11 to 21 bp across the mitogenomes but is exactly the same for two CR copies in a given genome. The same rule was observed for 13 crane species (Akiyama et al. 2017). However, the C. baudinii mitogenome showed an additional spacer with the length of 37 bp between the tRNA-Glu2 and the second CR (supplementary fig. S1, Supplementary Material online). The alignment of sequences between CR1 and tRNA-Thr2 with cytb revealed in all four parrots a sequence identical or highly similar to the 3′ end of cytb (supplementary fig. S1, Supplementary Material online). As the cytb2 sequences are always deprived of the 5′ region compared with the full-length cytb1, they are most likely nonfunctioned pseudogenes.

. 2.

Location of primer pairs and relative differences in the length of diagnostic fragments used to determine CR state as single or duplicated. Segments 15, 16, and ND6 correspond to the strategy proposed by Schirtzinger et al. (2012) for mitogenomes comprising: typical avian gene order (A), duplicated gene order identified for Amazona species (B), tandem duplication of the tRNA-Thr to CR segment (C and D). A strategy proposed in this study for identification of all variants of gene orders with duplicated CRs using the diagnostic Fragment 5 (E). Sequences of diagnostic primers D-F and D-R are listed in table 1. ND5, NADH dehydrogenase subunit 5 gene; cytb, cytochrome b gene; T, tRNA gene for threonine; P, tRNA gene for proline; ND6, NADH dehydrogenase subunit 6; E, tRNA gene for glutamic acid; F, tRNA gene for phenylalanine; 12S, 12S rRNA gene. Pseudogenes are marked by ψ and colored correspondingly to their functional gene copy.

The other gene order that also includes a duplicated region (GO-1 in fig. 1) was found in Pro. aterrimus goliath, E. roseicapilla, and Cac. moluccensis. This rearrangement seems to be a degenerated form of the fully duplicated rearrangement GO-FD. It is characterized by shortening and pseudogenization of the second tRNA-Glu as well as the specific degeneration of the second ND6 copy. This copy contains a purine-rich insertion, which is variable in length and is preceded by a fragment (433–450 bp) almost identical with the first ND6 copy followed by a highly degenerated region (supplementary fig. S1, Supplementary Material online).

The tandem duplications found in the mitogenomes of N. notabilis, Ny. hollandicus, C. baudinii, Cac. pastinator, and Cac. moluccensis make them longer compared with their previous genomic versions assuming the typical avian gene order. The ND2 gene in the newly obtained Cac. moluccensis genome does not correspond to the previous one (JF414239.1), which occurred to be identical with the ND2 gene found in the previously sequenced Cac. pastinator mitogenome (JF414240.1). However, the regions flanking the ND2 gene are identical when two versions of Cac. moluccensis mitogenome are compared (supplementary fig. S2, Supplementary Material online).

Comparison of the Duplicated Regions within Cacatuoidea and Strigopoidea Mitogenomes

The fully duplicated region GO-FD found in four parrot species is generally characterized by high similarity between paralogous sequences, that is, copies found within the same mitogenome. The second copies of cytb, tRNA-Thr, tRNA-Pro, ND6, and tRNA-Glu are identical or differ from their corresponding first copies only in one nucleotide in the case of cytb from Ny. hollandicus and Cac. pastinator (table 1; supplementary fig. S1, Supplementary Material online). Although the cytb sequences show the high similarity, the second copies are substantially truncated from their originals, by 662 bp (Ny. hollandicus) to 1,108 bp (N. notabilis), leaving only the 3′ ends. The second copy of ND6 in Ny. hollandicus has a deletion at the end of the sequence, which is associated with the change of the original stop codon TAG to the other stop codon AGG (supplementary fig. S1, Supplementary Material online). Two CRs show a slightly greater variation in identity, from 96.8% (N. notabilis) to 98.3% (Ny. hollandicus). The difference is mainly located at their 3′ ends, except for Ny. hollandicus, whose CRs differ only at their 5′ ends.

Table 1.

Comparison of Genes and CRs in the First and Second copies of Duplicated Region in Mitogenomes from Seven-Parrot Taxa Newly Sequenced in This Study.

Species	Accession Number	Copy	Length (bp)						Number/Percent of Residues Identical between Two Copies of:
			cytb	tRNA-Thr	tRNA-Pro	ND6	tRNA-Glu	CR	cytb	tRNA-Thr	tRNA-Pro	ND6	tRNA-Glu	CR
Nestor notabilis	MH133967	1st	1,140	69	69	513	71	1,224	32/100	69/100	69/100	513/100	71/100	1,180/96.8
		2nd	32	69	69	513	71	1,293
Nymphicus hollandicus	MH133968	1st	1,140	68	71	516	73	1,131	477/99.8	68/100	71/100	515/100	73/100	1,111/98.3
		2nd	478	68	71	516	73	1,243
Calyptorhynchus baudinii	MH133969	1st	1,140	68	69	519	70	1,195	122/100	68/100	69/100	519/100	70/100	1,157/97.2
		2nd	122	68	69	519	70	1,251
Probosciger aterrimus goliath	MH133970	1st	1,140	68	70	519	75	1,209	130/97.7	68/100	70/100	488/94.9	31/63.3	1,157/95.7
		2nd	133	68	70	576 (ψ)	49 (ψ)	1,242
Eolophus roseicapilla roseicapilla	MH133971	1st	1,140	68	69	519	75	1,201	120/98.4	68/100	69/100	496/96.9	36/58.1	1,142/95.2
		2nd	122	68	69	698 (ψ)	69 (ψ)	1,236
Cacatua moluccensis	MH133972	1st	1,140	68	69	507	75	1,186	106/92.2	67/98.5	69/100	480/95.0	20/60.6	1,118/94.9
		2nd	115	68	69	644 (ψ)	33 (ψ)	1,243
Cacatua pastinator pastinator	MH133973	1st	1,140	68	69	519	75	1,322	114/99.1	68/100	69/100	519/100	75/100	1,278/97.6
		2nd	115	68	69	519	75	1,360

Note.—Pseudogenes were marked by ψ.

The newly defined GO-1 gene order in three parrots is characterized by a larger variation of the duplicated copies (table 1). Only paralogous tRNA-Thr and tRNA-Pro genes are identical or almost identical with one-nucleotide difference in the case of Cac. moluccensis tRNA-Thr copies. The 5′ ends of paralogous CRs in the mitogenomes with gene order GO-1 show a slightly lower (94.9–95.7%) identity than those in the mitogenomes with the GO-FD order (table 1). The second cytb copies are shortened by 1,007 bp (Pro. aterrimus) to 1,025 bp (Cac. moluccensis) and show identity to their first copies from 92.2% (Cac. moluccensis) to 98.4% (E. roseicapilla). However, the most distinctive feature of the GO-1 rearrangement is the advanced degeneration of the tRNA-Glu and ND6 second copies. The tRNA-Glu2 sequences are shortened and identical to their first copies in only 58.1% (E. roseicapilla) to 63.3% (Pro. aterrimus). The second copies of ND6 have a long insertion within their initial part, but other regions (433–450 bp) are almost identical with their first copies. Therefore, the global identity between the paralogous ND6 genes varies from 94.9% for Pro. aterrimus to 96.9% for E. roseicapilla. Interestingly, the comparison of paralogous CRs revealed that in the case of all seven examined taxa, the second CRs (CR2s) are longer.

Mitogenomic Duplications in Other Parrots

The PCR strategy (fig. 2) was also conducted for other 12 Cacatuidae species (supplementary table S1, Supplementary Material online). In all cases, we obtained diagnostic fragments with estimated sizes from 1,900 bp (Lophochroa leadbeateri) to 2,250 bp (Cac. sulphurea sulphurea). The length of these fragments together with the PCR strategy implies that other cockatoo species also contain the tandem duplication of at least one CR because the primers used (supplementary table S1, Supplementary Material online) were designed to anneal the central regions of CRs in a way that prevents amplification of products from a single CR (fig. 2). Moreover, the similarity in the length of these fragments and analogous fragments 5 obtained for fully sequenced mitogenomes (supplementary table S1, Supplementary Material online) indicate that gene orders GO-FD and GO-1 or their degenerative variants should be also present in all Cacatuidae mitogenomes.

The complete mitochondrial genome of N. notabilis belonging to Strigopoidea also contained a duplicated region. However, the mitogenome of another member of this superfamily, Strigops habroptila, showed the typical avian gene order, that is, without any duplication (Harrison et al. 2004). This order was deduced only from the direct adjacency of tRNA-Thr and tRNA-Pro genes, which is not present in the Amazona mitogenome with the duplication. Nevertheless, such gene order does not rule out the possibility of duplication because the duplicated region of tRNA-Thr/tRNA-Pro/ND6/tRNA-Glu in N. notabilis contains also adjacent tRNA-Thr and tRNA-Pro genes. Therefore, we can assume that the mitogenome of S. habroptila also shares the same gene order and contains a duplication. The region in question was sequenced only partially from gene tRNA-Phe to CR, for which only the first 69 nucleotides were obtained. The misidentification of duplication in the Strigops mitogenome seems to be likely because the sequencing of this region proved difficult without the cloning of DNA fragments on account of the presence of repeats and heteroplasmy (Harrison et al. 2004). Interestingly, the 69-bp sequence annotated as the beginning of the CR differs substantially from all other CRs which were sequenced so far for Psittaciformes species and it rather resembles the GA-rich insertions identified in the present study within the ND6 copy of Pro. aterrimus, E. roseicapilla, and Cac. moluccensis. The mitochondrial gene order with this insertion (GO-1 in fig. 1) looked like a degenerated variant of the GO-FD gene order (fig. 1), in which pseudo ND6 and pseudo tRNA-Glu genes are present. The sequence identity between the pseudogenes and their corresponding functional copies was much higher for the ND6 gene (94–95.8%) than for tRNA-Glu gene (50.7–63.3%). However, similarity between tRNA-Glu copies was high enough to detect marks of the duplication. The GO-FD gene order, shown for the first time in our study for Psittaciformes, was previously found in Bucerotiformes, Gruiformes, Pelecaniformes, Procellariiformes, and Suliformes avian orders, for example (Abbott et al. 2005; Gibb et al. 2007, 2013; Cho et al. 2009; Eda et al. 2010; Morris-Pocock et al. 2010; Sammler et al. 2011; Zhou et al. 2014; Lounsberry et al. 2015; Akiyama et al. 2017; Rodrigues et al. 2017).

Extraordinarily long first CRs and specific gene orders in other parrot mitogenomes prompted us to search for potentially hidden pseudogenes in these regions (see section Reannotation of Gene Orders in text file in Supplementary Material online). The homology analyses revealed the presence of degenerated copies of cytb in P. luconensis, Psittacus erithacus, Poicephalus gulielmi and Neophema chrysogaster as well as ND6 in Forpus passerinus and Ps. erithacus (supplementary figs. S3 and S4, Supplementary Material online). Moreover, the Neo. chrysogaster mitogenome showed the first case of a remnant variant of CR1 (supplementary fig. S4, Supplementary Material online). So far, only vestigial CR2s were reported. As a result of these reannotations, we determined four additional gene orders in mitogenomes: GO-2, GO-3, GO-5, and GO-7 (fig. 1).

The newly sequenced and reannotated mitogenomes increased the number of parrot genomes with duplicated CRs from 16 to 24 and enabled their comparison to draw some general conclusions about these regions (see section Comparison of Length and General Structure of Duplicated CRs in text file in Supplementary Material online). We found that the second CRs are longer than the first ones in 16 out of 21 mitogenomes (supplementary table S2, Supplementary Material online). This tendency seems common because it is also present in other groups of birds (supplementary table S3, Supplementary Material online). Out of 121 pairs of the avian CRs, the CR2 is longer in 109 cases. The median lengths of CR1 and CR2 for all bird mitogenomes are 1,136 and 1,265 bp, respectively, and the difference is statistically significant (supplementary table S4, Supplementary Material online).

Phylogeny and Evolution of Duplicated Regions in Parrots

To infer evolutionary relationships between CRs, especially those that were duplicated, we performed phylogenetic analysis using all known sequences of parrot CRs. Three methods produced almost identical tree topologies (fig. 3). Although not all relationships between groups at higher taxonomic levels were well recovered because of a high divergence rate of these sequences, phylogenies at the lower taxonomic level were reliably resolved. In the tree including all available CR sequences, 18 out of 20 paralogous CRs were significantly clustered together (fig. 3). Only the relationships within Amazon parrots were not well resolved. Therefore, we carried out additional analyses including longer sequences only from Amazona and Pionus selected as an outgroup. The inferred phylogenetic trees confirmed closer relationships between two CRs present in the mitochondrial genome from the same species, that is, between paralogous regions (CR1 and CR2), rather than corresponding regions from different species, that is, orthologous regions (CR1s or CR2s) (fig. 3). Our results based on sequences from Amazon parrots correspond to those obtained by Eberhard et al. (2001).

. 3.

Phylogenetic trees obtained in PhyloBayes based on all available parrot sequences of CRs (A) and those from Amazona and Pionus (B). The blue and red colors correspond to the first and the second copy of CR, respectively. The values at nodes, in the order shown, indicate posterior probabilities found in PhyloBayes (PB) and MrBayes (MB) as well as SH-aLRT (SH) and nonparametric bootstrap (BP) percentages calculated in IQ-TREE. The posterior probabilities <0.5 and the percentages <50% were omitted or indicated by a dash “-.”

All the phylogenetic analyses and the high similarity between two CRs in the same mitogenome molecule indicate that the duplication of these regions or at least a large portion of them are specific to species lineages in all parrots studied so far and are subjected to concerted evolution. Several mechanisms of this phenomenon were proposed (Kumazawa et al. 1996, 1998; Arndt and Smith 1998; Eberhard et al. 2001; Abbott et al. 2005; Shao et al. 2005; Gibb et al. 2007; Kurabayashi et al. 2008; Cadahia et al. 2009; Eda et al. 2010; Morris-Pocock et al. 2010; Sammler et al. 2011; Schirtzinger et al. 2012) and the most probable involve: 1) gene conversion of only selected conserved parts and independent evolution of intervening sections, or 2) gene conversion of the whole region followed by stronger selection on the functional part and quick divergence or degeneration of other parts. An alternative possibility assuming independent recent duplications in each lineage is much less likely because it would require a large number of such events followed by the same pattern of replication errors, duplication and degradation (Shao et al. 2005; Morris-Pocock et al. 2010; Zhou et al. 2014) as well as the entire duplicated region should have shown evenly distributed high sequence similarity (Eberhard and Wright 2016).

Our findings show that the homogenization of duplicated regions occurs within individual species lineages at a frequency higher than the rate of speciation. Eberhard et al. (2001) estimated the rate of this process is once in 34,670 ± 18,400 years in the case of Amazona mitogenomes. Interestingly, this process can take only 850 ± 479 years in mangrove killifishes (Kryptolebias marmoratus) (Tatarenkov and Avise 2007) or even can occur in every generation as it was suggested for Philippine hornbills (Sammler et al. 2011).

The concerted evolution seems to be a universal feature of avian mitogenomes and was also suggested for other avian groups: albatrosses (Abbott et al. 2005; Eda et al. 2010), spoonbills (Cho et al. 2009), shorebirds (Verkuil et al. 2010), boobies (Morris-Pocock et al. 2010), Philippine hornbills (Sammler et al. 2011), ardeids (Zhou et al. 2014), cranes (Akiyama et al. 2017), passerines (Gibb et al. 2015; Caparroz et al. 2018) as well as gannets, darters, and cormorants (Gibb et al. 2013). The analyses showed that not only CRs are subjected to concerted evolution but also the genes adjacent to the CR, such as tRNA-Thr, tRNA-Pro, ND6 and tRNA-Glu, and a portion of cytb are evolving in concert. Our analyses also revealed a high similarity between paralogous gene copies in the same mitogenome, which also suggests their concerted evolution. However, in contrast to CRs, many duplicated genes are truncated and have deletions, indicating that they are more often or faster subjected to pseudogenization and final loss from the genomes. Exceptions are mitogenomes of Nestor, Calyptorhynchus, Nymphicus, and Cac. pastinator, where only cytb was pseudogenized.

Relationships between Mitogenomic Duplications and DNA Strand Asymmetry

Asymmetry between DNA strands of mitochondrial genomes was measured by normalized differences in complementary nucleotides, that is, AT-skew and GC-skew (Perna and Kocher 1995; Reyes et al. 1998). These parameters assess differences in directional mutational pressure associated with replication between the heavy (H) and light (L) strands in the mitochondrial genome. During replication, the H strand stays for a longer time in a single-stranded state, in that it is more prone to accumulation of substitutions C→U and A→hX (hypoxanthine), as a result of spontaneous deamination of C and A. This leads to compositional asymmetry between these strands. It is expected that the bias is stronger in genomes with slower replication when their H strand remains exposed to mutations for a longer time (Reyes et al. 1998).

Our analyses of all sites in 47 parrot mitogenomes showed no significant difference in AT-skew and GC-skew between mitogenomes with and without the duplicated region (supplementary table S5, Supplementary Material online and fig. 4). However, the difference became significant when all sites in protein-coding sequences (PCGs) and only 4-fold degenerate sites (4FD) were analyzed (supplementary table S5, Supplementary Material online and fig. 4). The mitogenomes with the duplicated region are characterized on average by a greater bias in the two asymmetry parameters than the mitogenomes with the single region. The difference is statistically significant for GC-skew in all PCG and 4FD sites, marginally significant for AT-skew in all PCG sites but not significant for AT-skew in 4FD sites. Eberhard and Wright (2016) did not find any significant difference in the asymmetry for the entire genome either but only for AT-skew when all sites from the H-strand-transcribed genes were considered.

. 4.

Box plots of AT-skew (ATs) and GC-skew (GCs) calculated for the entire genome as well as all sites in PCGs and 4-fold degenerate (4FD) sites in mitogenomes with duplicated (D) or single (S) regions. The thick line indicates median, the boxes show quartile range, and the whiskers denote the range without outliers. P-values of the null hypothesis assuming no differences in DNA asymmetry between the mitogenomes are shown.

The weaker compositional bias in the entire genomic sequence and the lack of a significant difference between mitogenomes with and without duplicated regions result likely from the stronger selective constraints on genes encoding rRNA and tRNA, whose sequences are present in the analysis of the whole genomic sequence. Mitochondrial PCGs are characterized by a greater substitution rate and can accumulate more mutations associated with replication than the RNA genes (Pacheco et al. 2011; Urantowka et al. 2017b). On the other hand, DNA asymmetry is most pronounced in 4FD sites because they are subjected to weaker selection pressure although they can be also under the influence of selection for translational efficiency (Jia and Higgs 2008; Sun et al. 2009; Lin et al. 2012; Brule and Grayhack 2017) and selection acting at the amino acid level (Morton 2001; Błażej et al. 2017).

In agreement with findings by Eberhard and Wright (2016), the presented results indicate that the duplicated region can extend time in which the H strand stays in a single state. It may be associated with slowing the replication of mitochondrial DNA, extending the replication time by the additional duplicated region or increasing the number of nascent H strands due to the larger number or extension of the replication initiation region. It is not inconceivable that the second CR can initiate a next replication round before the first one is finished, similar to bacteria that divide quicker than its chromosome replicates (Cooper and Helmstetter 1968; Skarstad et al. 1986). It would be interesting to experimentally verify these hypotheses.

Relationships of Mitogenomic Duplications with Morphological Features and Life Span

We tested several hypotheses that asked whether possessing mitogenomes with a single or duplicated region might be associated with basic morphological features (fig. 5 and supplementary table S6, Supplementary Material online). The analyses showed that parrots with a duplicated region in their mitogenome are characterized by significantly greater (more than two times) body mass (M) and wingspan (W) but shorter tails (T) than parrots without such region. These two groups do not differ in the total body length (L) but if we compare their body length without the tail contribution (L – T), the parrots with duplication are longer than those without a duplicated region. After the normalization of T and W by M and the inclusion of the body length without the tail (L – T), the parrots with a duplicated region also have a significantly shorter tail relative to M and L. However, their wingspan, longer before the normalization, becomes significantly shorter in relation to the body mass (W/M) but not to the body length (W/[L – T]).

. 5.

Box plots of morphological features for parrots with mitogenomes containing duplicated (D) or single (S) regions. M, body mass; W, wingspan; T, tail length; L, total body length. The thick line indicates median, the boxes show quartile range, and the whiskers denote the range without outliers. P-values of the null hypothesis assuming no differences in a given feature between the parrots are shown.

As the animals’ respiration metabolic rate shows negative allometric relationship with body mass (Kleiber 1932; West et al. 1997; Dodds et al. 2001; Glazier 2005; White and Kearney 2014), larger-bodied parrots should have lower body mass-specific metabolic rates compared with smaller-bodied parrots. If we assume that mitogenomes with the duplicated CR can increase the overall number of genome copies per mitochondrion and efficiency of energy production, the parrots with larger body mass might benefit from the duplication to increase their metabolic rate compared with the expectation. The more efficient energy production in the parrots with the duplicated CR might also compensate their smaller wingspan in relation to body mass and the greater effort that should be put into flight of the heavier birds. Generally, as avian body size increases, flight performance declines and total flight power does not increase as fast as body mass does (Tucker 1973; Greenewalt 1975; Rayner 1995; Lovette and Fitzpatrick 2016).

The other distinguishing feature of parrots with the duplication is that the short tail is significantly associated with not only body mass and body length, but also possibly high speed of flying because the elongation of the tail increases the drag that resists a bird’s forward motion decreasing their flight speed and increasing metabolic costs (Balmford et al. 1993; Thomas 1993; Norberg 1995; Carey 1996; Hedrick and Biewener 2007; Clark and Dudley 2009; Lovette and Fitzpatrick 2016). On the other hand, long tails with stiff feathers can generate some lift as well as help in maneuvering and steering (Carey 1996; Evans 1998; Rowe et al. 2001; Matyjasiak et al. 2004; Hedrick and Biewener 2007). These findings could imply that parrots with the duplicated CR in their mitochondrial genomes are more active aviators flying with higher speed, which may be associated with more efficient production of energy by mitochondria containing a larger number of replicating genomes per organelle. However, to confirm the above-mentioned relationships, a better sampling of parrot species is necessary, especially in groups where mobility, dispersal capacity, and migratory status vary a lot among closely related species.

The parrots possessing mitogenomes with the duplication are also characterized by significantly greater longevity, measured by the maximum life span, than parrots without the duplication (fig. 6, supplementary table S6, Supplementary Material online). However, it may be a consequence of larger body size because body mass is positively correlated with longevity in homeothermic animals, including birds (Lindstedt and Calder 1976, 1981; Hulbert et al. 2007; Skujina et al. 2016). In our data set, the relationship between maximum life span Smax and body mass M was described by the equation ln Smax=0.392⋅ln M + 1.091 with adjusted R-squared = 0.453 and all coefficients significant with P < 2.2E-16. Separate linear modeling for parrots with and without duplicated regions in their mitogenomes demonstrated that the linear function for the first group is characterized by a larger slope and intercept, that is, 0.388 ± 0.048 (SE) and 1.224 ± 0.264, than those in the second group, that is, 0.379 ± 0.073 and 1.079 ± 0.367, respectively. The coefficients were significantly greater than zero with P < 0.0048. The results might suggest that parrots with the duplication generally live longer in relation to their body mass than parrots without the duplication. Similarly, a study of 92 avian families including parrots also demonstrated that longer‐lived avian species tend to have mitogenomes with the duplicated CR (Skujina et al. 2016). Such organisms with the extra copies might live longer because it may protect their mitochondrial DNA from loss of the region due to deletion or increasing flexibility of mitochondrial response to environmental changes associated with elevated metabolism. However, additional direct biochemical and physiological analyses are necessary to prove these hypotheses.

. 6.

Box plots of the maximum life span in years for parrots with mitogenomes containing duplicated (D) or single (S) regions. The thick line indicates median, the boxes show quartile range, and the whiskers denote the range without outliers. P-value of the null hypothesis assuming no differences in the maximum life span between the parrots is shown.

It was also shown for birds that body size and metabolic rate may additively contribute to an increase in overall substitution rate under selection for small body and/or population size (Berv and Field 2018). This increase can result from bigger production of mutagenic oxygen radicals in smaller-bodied organism with higher body mass-specific metabolic rates (Gillooly et al. 2005, 2007). Assuming that the duplicated regions in mitogenomes could also cause a more intense production of energy by mitochondria and free radicals, we should expect an elevated mutation and substitution rate in such genomes. However, the relationship between the substitution rate and the presence or absence of the duplicated region was not found (Eberhard and Wright 2016). An additional study based on larger sampling can verify the relationship but this analysis should take into account that the observed substitution rate is also under a strong influence of body and population size, which changed during the evolution of Aves (Berv and Field 2018).

Mapping Mitogenomic Duplications onto the Phylogenetic Tree of Parrots

To analyze mitogenomic duplications from the evolutionary point of view, we inferred phylogenetic relationships between as many members of Psittaciformes order as possible using 245 sequences of NADH dehydrogenase subunit 2 (ND2) gene, which is one of the most commonly used molecular markers for parrots (supplementary figs. S5–S8, Supplementary Material online). However, the single-locus phylogenies can suffer from a poor phylogenetic signal. Therefore, we also constructed phylogenetic trees based on the concatenated alignment of 37 mitochondrial genes to maximize the number of informative sites (supplementary figs. S9–S12, Supplementary Material online). Three approaches, two Bayesian and one maximum likelihood (ML), produced trees with very similar and quite well-resolved relationships. The same main taxonomic groups from superfamily to tribe levels can be well recognized in these trees. The relationships between the main groups correspond well to recent large-scale phylogenomic studies including many parrot species (Provost et al. 2018) and representatives of major parrot clades (Prum et al. 2015). However, some groups are differently clustered. Two Bayesian approaches, MrBayes and PhyloBayes, grouped the clade of Amoropsittacini and Forpini with Androglossini, while IQ-TREE with Arini in the ND2 trees (supplementary figs. S5–S8, Supplementary Material online). In all three trees based on 37 mtDNA genes, representatives of Amoropsittacini were absent and Forpini clustered directly with Androglossini (supplementary figs. S9–S12, Supplementary Material online). It contrasts with the sisterhood between Forpini and Arini obtained by other authors (Schweizer et al. 2014; Provost et al. 2018); however, the relationships between the above-mentioned tribes are not well resolved in any phylogenetic trees. It may result from a rapid divergence of these lineages, which is demonstrated by short internal branches in the trees (supplementary figs. S6–S8 and S10–S12, Supplementary Material online). Depending on the method used, different relationships were also obtained within Psittaculidae. Other approaches also provided inconsistent or unresolved tree topologies (Schirtzinger et al. 2012; Schweizer et al. 2013; Provost et al. 2018). In our study of ND2 gene, members of Agapornithinae created a monophyletic clade in the MrBayes tree (supplementary fig. S7, Supplementary Material online), while in the IQ-TREE (supplementary fig. S6, Supplementary Material online) and PhyloBayes (supplementary fig. S8, Supplementary Material online) trees this group was separated into two lineages: Agapornis and Loriculus with Bolbopsittacus. All methods and data sets showed the subfamily Platycercinae polyphyletic because its two tribes Pezoporini and Platycercini were separated in the phylogenetic trees. Such separation was highly supported, especially in the trees based on 37 mtDNA genes (supplementary figs. S9–S12, Supplementary Material online). These tribes did not form a monophyletic clade, as was the case in trees obtained by Schirtzinger et al. (2012) and Provost et al. (2018) in contrast to Schweizer et al. (2013). Relationships between other groups and within tribes were consistent with the results of other researches.

We mapped collected data on the presence and absence of the mitogenome duplication onto the tree topology obtained for ND2 in one of the approaches (IQ-TREE) and inferred the ancestral states by ML and maximum parsimony (MP) methods for two data sets, one including 238 parrot species with both known and unknown duplication states (supplementary figs. S13 and S14, Supplementary Material online) and the other containing 141 species with a verified duplication (figs. 7 and 8). These methods provided the same results and all approaches indicated that the mitogenome of the common ancestor of all known parrot taxa contained a duplication. The ML method provided the probability of 0.997 for this state. The duplication state was inherited by the ancestor and descendants of early diverged lineages, Strigopoidea and Cacatuoidea as well as Cacatuoidea+Psittacoidea, with the probability >0.997. The mitogenome of the Psittacoidea ancestor could also contain a duplication with the probability of 0.963, which was lost in the ancestor of Psittaculidae with the probability of 0.981. In the inference by the ML method, however, the duplication occurred again in the Loriinae with the probability of 0.966 and some other lineages, that is, Neophema chrysogaster, Loriculus vernalis, and Prioniturus with the probability of 0.967. If we assume alternative relationships within Psittaculidae obtained in our Bayesian trees and by other authors (Schirtzinger et al. 2012; Schweizer et al. 2013), the general result stays the same. However, the probability of single region in the ancestor of Psittaculidae decreases to 0.67 supposing that Lorinae is sister to Agapornithinae, while Platycercinae and Psittaculinae cluster with Psittacellinae (Schweizer et al. 2013).

. 7.

Maximum likelihood reconstruction of ancestral states and mapping of mitogenomic duplications onto the ND2 phylogenetic tree of parrots with verified presence or absence of the duplication. The area of colors at nodes corresponds to the probability of the given state, single or duplicated region. The probability value for a more probable state was given.

. 8.

Maximum parsimony reconstruction of ancestral states and mapping of mitogenomic duplications onto the ND2 phylogenetic tree of parrots with verified presence or absence of the duplication. Two-colored branches correspond to the equal probability of two states, single or duplicated region.

The duplication was also maintained up to the ancestor of Psittacidae with the probability of 0.980 and Arinae with the probability of 0.947 but was lost independently with the probability of 0.967 in the lineage of Brotogeris and Myiopsitta as well as Amoropsittacini and core Arini after divergence of Deroptyus and Pionites (figs. 7 and 8). The adoption of different phylogenetic relationships within Arinae obtained by other authors (Schweizer et al. 2014; Provost et al. 2018) did not change the conclusion about the independent loss of the duplication in these three lineages and the presence of the duplicated region in the ancestral mitogenome of Arinae.

Analyses of complete mitochondrial genomes and mapping their duplicated features onto 37-gene phylogeny also demonstrated that the parrot ancestor contained a duplication in its mtDNA genome with the probability of 0.992 (fig. 9). There was no difference in the result of the ML and MP methods. This state was maintained in the ancestors of subsequently diverged from main parrot lineages with the probability >0.978 to 0.997. According to the result of the ML inference with the complete mitochondrial genome data, the duplication was lost independently at least five times in representatives of Psittrichasidae, Psittaculidae, and Psittacidae. The biggest number of genomes with the single region are classified to Arini, whose ancestor could lose the duplication with the probability of 0.926. However, there is a difference between the mapping duplications on the ND2 and 37-gene phylogenies. In the first case, the ancestor of Psittaculidae likely lost a duplicated region, which occurred secondarily in the ancestor of Lorinae as well as some members of Platycercinae, Psittaculinae, and Agapornithinae (figs. 7 and 8). In the second case, the ancestor of Psittaculidae (and also Psittrichasidae) still contained a duplicated region, which was later lost in some genomes (fig. 9). Additional analyses are required to confirm one of these scenarios.

. 9.

Maximum likelihood (A) and maximum parsimony (B) reconstruction of ancestral states and mapping of mitogenomic duplications onto the parrot phylogenetic tree based on 37 mitochondrial genes. Other explanation as in figures 7 and 8.

The structure of the sequenced duplicated region and the distribution of individual gene orders in the phylogenetic tree indicated that the mitogenome of the parrot ancestor could comprise copies of all genes (GO-A, fig. 1). The cytb gene was probably the first to degenerate (GO-FD in fig. 1) as it was observed in the earliest evolved lineages of parrots, Nestor, Calyptorhynchus and Nymphicus as well as in the later derived Cac. pastinator. The more degenerated level of the duplicated region with degraded second copies of ND6 and tRNA-Glu (GO-1 in fig. 1) is present in three Cacatuinae members. The different way it has degraded involving long insertions in gene ND6 suggests an independent pseudogenization process. The degeneration of the duplicated region in the P. luconensis (Psittaculidae) mitogenome was also independent because it involved the first copy of ND6 and tRNA-Glu instead of the second one (GO-2 in fig. 1). A slightly more advanced degeneration involving three tRNA genes is represented by closely related Ps. erithacus and Po. gulielmi, the earliest branch of Psittacidae (GO-3 in fig. 1). Pionus and Amazona belong to much later diverged Psittacidae group and create a monophyletic clade sharing the same type of gene order (GO-4 in fig. 1) with the loss of one cytb gene, one tRNA-Thr gene, and one tRNA-Pro gene, which suggests that this process could have begun in their common ancestor. Forpus passerinus, representing an independent lineage of Psittacidae, lost furthermore a copy of the tRNA-Glu gene (GO-5 in fig. 1). The most degenerated duplicated region is present in Melopsittacus undulatus which in addition to the above-mentioned genes also lost one ND6 gene (GO-6 in fig. 1). A unique degeneration way of the duplicated region may be exhibited by Neophema chrysogaster (Psittaculidae), which maintained a pseudogene of cytb but lost the internal part of the first CR besides other genes (GO-7 in fig. 1).

Conclusions

The presence of mitogenomes with the duplicated region in the earliest diverged parrot lineages, Strigopoidea and Cacatuoidea, as well as superimposing duplication data from all parrot mitogenomes onto their phylogeny indicate that the duplication could be a plesiomorphic state for the common ancestor of known parrots. The state was inherited by subsequent ancestors of main parrot groups and was then lost at least four times in some lineages. In three cases, the duplication in the mitogenome reoccurred secondarily. Our findings imply that once duplications appeared in the ancestral mitochondrial genomes, they tended to evolve in concert rather than independently, most likely by gene conversion. The CRs were subjected to the homogenization, while adjacent duplicated genes were usually prone to pseudogenization and loss. The concerted evolution occurred quite often within individual species’ lineages at a frequency higher than the rate of speciation. This process made the CRs identical or very similar.

It must be noted that the proposed scenario is most parsimonious for the current data but may change when additional mitogenomes of parrots are checked for a duplicated region. We can rather expect that more mitochondrial genomes of parrots will reveal the presence of duplication undetected so far, which will support the view that the ancestor of Psittaciformes and the main parrot lineages contained a mitogenomic duplication. However, it is difficult to distinguish individually between a duplication inherited from an ancestor and a newly obtained one because the concerted evolution erases the differences accumulated between the copies since their duplication in the past.

The statistical models assumed the equal probabilities of gain and loss of duplications in mitogenomes. However, the real gain and loss probabilities can be different. The potential bias can be modified on two levels, mutational and selective. On one hand, the bias can result from different probabilities of duplications and deletions of DNA fragments in the genomes due to molecular mechanisms. On the other hand, selection can favor mitogenomes with duplicated regions, for example, because of more efficient energy production, or eliminate them, for example, because of a longer replication time or a greater energy consumption in the synthesis of the longer genomes. More mitogenomes studied in the context of duplication are necessary to estimate the real duplication gain/loss ratio and more detailed knowledge about these processes can give deeper insight into the mitogenome evolution.

It is interesting to consider whether other bird groups that are sister to parrots also had an ancestral mitogenomic duplication and whether this feature may have arisen earlier in avian evolutionary history. The most comprehensive phylogenetic analyses indicate that Passeriformes is the most closely related order to Psittaciformes (Jarvis et al. 2014; Prum et al. 2015). In passerines, various gene rearrangements were also detected, including duplications (Mindell et al. 1998; Bensch and Harlid 2000; Gibb et al. 2015; Caparroz et al. 2018). Out of 298 passerine mitogenomes annotated in GenBank, 72 contain the duplication, 202 do not have it, and 24 are incomplete. Two species have a unique duplication cytb/tRNA-Thr1/tRNA-Pro1/ND6-1/tRNA-Glu1/CR1/tRNA-Pro2/ND6-2/tRNA-Glu2/CR2 (Gibb et al. 2015). Forty-six species share the same gene order with the parrot M. undulatus (GO-6 in fig. 1), while 24 species differ from the GO-6 in the presence of a remnant CR2. However, the current study of the distribution of mitogenomic characters on the phylogeny of passerines suggests that their ancestor did not contain the mitogenomic duplication, which occurred independently later at least six times in some lineages (Caparroz et al. 2018). It would imply that the duplication was plesiomorphic only for parrots. On the other hand, Psittaciformes and Passeriformes are related to Falconiformes, whose all known mitogenomes comprise exclusively a duplication (Gibb et al. 2007; Ryu et al. 2012; Dou et al. 2016; Wang, Zhang, et al. 2016; Sveinsdottir et al. 2017). Besides parrots and falcons, the single origin of the duplicated region was also proposed for cranes (Akiyama et al. 2017), sulids (Morris-Pocock et al. 2010), and ardeids (Zhou et al. 2014). Interestingly, all Gruiformes and Suliformes mitogenomes containing the duplication maintain also the full set of all duplicated gene copies, that is, the gene order GO-FD, which was most likely present in a potential parrot ancestor. This gene rearrangement was also identified in some taxa belonging to Bucerotiformes and Procellariiformes (Abbott et al. 2005; Gibb et al. 2007, 2013; Sammler et al. 2011; Lounsberry et al. 2015) as well as the majority of Pelecaniformes (Zhou et al. 2014). The duplication at various levels of degeneration was also found in representatives of Accipitriformes (Gibb et al. 2007), Coraciiformes (Huang et al. 2016), Cuculiformes (Pratt et al. 2009; Pacheco et al. 2011; Wang, Liang, et al. 2016), Phoenicopteriformes (Morgan-Richards et al. 2008; Luo et al. 2016), Piciformes (Gibb et al. 2007), and Strigiformes (Hanna et al. 2017). Therefore, it is not inconceivable that the mitogenomic duplication can be a plesiomorphic feature not only for parrots but also for other avian groups or even the whole Neoaves. On the other hand, there are many Neoaves orders whose members contain exclusively mitogenomes without the duplication: Caprimulgiformes, Cathartiformes, Columbiformes, Eurypygiformes, Otidiformes, Phaethontiformes, Podicipediformes, and Trogoniformes. These groups are spread over the whole avian phylogenetic tree and mix with those containing the mitogenomic duplication. It cannot be ruled out that they can contain a hidden duplication. Then, an extensive survey of mitogenomes from many bird representatives is necessary to reconstruct the evolution of mitochondrial genome duplication in Aves at the global scale.

The preservation of two CRs can give selective advantage due to more efficient initiation of replication or transcription resulting from more replication and transcription events starting per unit time, which can increase the overall number of genomic and transcript copies per mitochondrion (Kumazawa et al. 1996; Arndt and Smith 1998; Tang et al. 2000; Umeda et al. 2001; Jiang et al. 2007). Actually, it was demonstrated that snake and human mitogenomes with two CRs can replicate more efficiently than variants with only one CR (Kumazawa et al. 1996; Umeda et al. 2001). Moreover, it was observed that the population of mitogenomes in human cells being initially a mixture of genomes with one and two CRs, shifted over time toward genomes with two CRs (Tang et al. 2000). Two identical or similar copies of CR generated by gene conversion can make the regulation of replication and transcription also more effective and consistent because these processes can be then controlled by the same enzymes and factors. In the case of parrots, the duplications can also compensate for a slow rate of genome replication (Eberhard and Wright 2016). Duplication variants including ND6 may be especially favored because the majority of initiation events can be located in this gene as was shown for the chicken mitogenome (Reyes et al. 2005).

As a result of duplications, energy production by mitochondria can be more efficient and higher. Our study demonstrated that keeping two copies of CR can be associated with a longer life span and morphological features related to more active flight consuming a lot of energy. Molecular dating and phylogeography analyses demonstrated that long-distance and even transoceanic dispersals well explain the present distribution and diversification of parrots (Wright et al. 2008; Schweizer et al. 2010; White et al. 2011). It suggests that ancestors of modern parrots must have been much better fliers than the present species, which are not long-distance migration birds, on the continental and intercontinental scales. It corresponds well to our findings showing that the ancestors contained mitogenomes with an additional copy of the CR, which can be associated with more efficient functioning of mitochondria and production of the energy desirable in active flight. In the light of this argument, it is interesting to consider reasons to lose the duplication in some parrot lineages. This process could occur without any selective benefits via neutral evolution or be associated with selection for genome compactness and increase the speed of replication. Alternatively, it could be related to the evolution of parrots with less active and short-distant flight type or favoring a more sedentary lifestyle, spending more time perching or climbing in trees. Nevertheless, these hypotheses wait for verification in biochemical and physiological investigations.

Materials and Methods

Samples and DNA Extraction

Samples from parrots were obtained from Rosewood Bird Gardens & Breeding Farms, in Rosewood, Australia; Loro Parque Foundation, Avda Puerto de la Cruz, Tenerife, Spain; and Zoological Garden in Wrocław, Poland (supplementary table S1, Supplementary Material online). Blood samples were taken as dry blood spots on a fiber filter for laboratory analysis and were preserved in parafilm-sealed Eppendorf tubes at −20 °C until use to avoid dampness. Shafts of freshly collected breast feathers were cut off and stored in the same conditions until use. Total DNA was extracted from both tissue types with Sherlock AX Kit (A&A Biotechnology) according to the manufacturer’s protocol.

PCR Strategy for Tandem Duplication Survey

The amplification of diagnostic fragments (Seg 15, Seg 16, and Seg ND6 in fig. 2) varying in length depending on the presence or absence of a duplicated CR was used to screen 117 parrot species from 79 genera. The chosen strategy allows to distinguish the typical avian gene order (fig. 1) from that identified in Amazona sp. (GO-3 in fig. 1). However, this strategy may miss the duplicated regions when they are nearly identical and comprise tRNA-Thr/tRNA-Pro/ND6/tRNA-Glu genes followed by a CR (GO-FD in fig. 1). Such gene order would be expected considering the results for P. luconensis mitogenome (Eberhard and Wright 2016), in which a degenerated rearrangement was found. In the case of such duplication, amplified diagnostic PCR fragments should have the same length as those obtained for the typical avian gene order. Moreover, the identical length of two copies of segment 16 (fig. 2) makes it impossible to distinguish them. As it was shown in figure 2, only the increased elongation time can allow detection of the presence of two copies of ND6 gene and CR. However, the presence of two copies of some genes within the long Seg 15 and long Seg 16 amplicons makes it difficult or even impossible to establish the order of these genes. Moreover, this PCR strategy will produce not only long variants but also many short variants perturbing the sequencing reaction. Therefore, to verify the presence of tandem duplication within the mitochondrial genome of selected parrots, we used a strategy, which was proposed by Gibb et al. (2007) and is presented in figure 2. We designed primers (supplementary table S1, Supplementary Material online) to amplify the appropriate fragment 5 between two potential CRs. As the paralogous CRs are usually nearly identical as a result of concerted evolution (Eberhard et al. 2001; Akiyama et al. 2017), suitable primers were designed to anneal the central regions of CRs. This strategy is diagnostic and the amplification occurs only when two CRs are present in the genome. In the case of a single CR copy no product is amplified (fig. 2).

PCR Strategy for Amplification of the Whole Mitochondrial Genomes

We amplified the whole mitochondrial genomes of N. notabilis, Ny. hollandicus, C. baudinii, Cac. moluccensis, Cac. pastinator, Pro. aterrimus goliath, and E. roseicapilla (table 1). They were amplified in six overlapping fragments. As fragments 5, designed for tandem duplication survey, comprise potentially the second part of CR1 and the first part of CR2, appropriate PCR reactions were performed to complete the missing parts of CR1 and CR2 (fig. 2). The first reaction amplified fragment 4 including the whole ND5 gene, while the second one produced fragment 6 with the complete 12S rRNA gene. The appropriate elongation times were used to avoid amplification of long variants (about 2 kb longer) of fragments 4 and 6, which would contain two copies of some genes and/or CRs (fig. 2). Fragment 5 is the only one whose length depends on the presence of genes located between two potential CRs. The other five fragments used to amplify the whole genome were designed to be from 4 to 5 kb in length in order to minimize the potential amplification of nuclear mitochondrial DNA inserts (NUMTs) whose average size is usually below 1 kb (Richly and Leister 2004).

DNA Amplification and Sequencing

The PCR amplifications were performed in 25 μl reaction mixture containing 50 ng of the DNA template, 1U DreamTaq Green DNA Polymerase (Thermo Fisher Scientific), 2.5 μl of 10× buffer, 0.6 μl of 10 mM dNTPs, and 0.6 μl of each primer (10 μM). The reaction conditions were as follows: 94 °C for 5 min; 94 °C for 30 s, 58 °C for 30 s, 72 °C for 180 s repeated 35 times; and 72 °C for 5 min in the case of fragments 1, 2, 3, and 4. For fragments 5 and 6, we used the following program: 94 °C for 5 min; 94 °C for 30 s, 58 °C for 30 s, 72 °C for 120 s repeated 35 times; and 72 °C for 5 min. For each fragment, the appropriate amount of the PCR reaction mixtures was cleaned with the use of Clean-up Kit (A&A Biotechnology) to obtain the final volume of 100 μl with the concentration of at least 50 ng/μl. Two DNA strands of the cleaned PCR products were sequenced using Primer Walking method (Wyzer Biosciences Inc., Cambridge, MA). Overlaps between the amplified fragments from 1 to 6 were sufficient to assemble the whole circular genomic sequence with the use of appropriate reference mitogenomes of Cacatuidae taxa or N. notabilis. The annotation of genes was performed in MITOS (Bernt et al. 2013). Intergenic and long CRs were searched for remnant CRs or pseudogenes using the optimal global: local algorithm (glsearch) form FASTA package version 36.3.8g (Pearson et al. 1997).

Phylogenetic Analyses

Phylogenetic relationships between parrot species were inferred using two sets: the alignment of 245 nucleotide sequences of NADH dehydrogenase subunit 2 (ND2) gene and the concatenated alignment of 13 protein-coding genes, 12S and 16S rRNA as well as 22 tRNA sequences from 48 parrot species and five passerines used as an outgroup (see section Phylogenetic Analyses in text file in Supplementary Material and supplementary tables S7 and S8, Supplementary Material online). Phylogenetic analyses of CRs were also based on two sets: all available 67 CR sequences and 18 CRs from Amazona and Pionus. The alignments were obtained in MAFFT (Katoh and Standley 2013), and sites suitable for phylogenetic study were selected in GBlocks (Talavera and Castresana 2007). Three phylogenetic approaches were applied: the ML method in IQ-TREE (Nguyen et al. 2015), as well as two Bayesian analyses in MrBayes (Ronquist et al. 2012) and PhyloBayes (Lartillot and Philippe 2004). The best-fitted substitution models considering all possible codon and gene partitions in the ND2 and 37-gene alignments were found in ModelFinder (Chernomor et al. 2016; Kalyaanamoorthy et al. 2017) and PartitionFinder (Lanfear et al. 2012). The best model for CR alignments was used according to jModelTest 2.1 (Darriba et al. 2012). We also applied mixed models in MrBayes and the CAT model in PhyloBayes. The full list of the best substitution models found and used for the phylogeny construction is included in section Phylogenetic Analyses in text file in Supplementary Material and supplementary table S8, Supplementary Material online. SH-aLRT test assuming 10,000 replicates and nonparametric bootstrap with 1,000 replicates were used in IQ-TREE.

Data about the presence or absence of duplication in the parrot mitogenomes were collected based on PCR screening obtained by Schirtzinger et al. (2012) and ourselves as well as genomes deposited in GenBank and obtained in this study (supplementary table S7, Supplementary Material online). The data were mapped on the IQ-TREE tree obtained for ND2 gene and mtDNA markers using Mesquite (Maddison WP and Maddison DR 2017). Two ND2 data sets were considered. The first data set included 238 parrot species with known and unknown information about the mitogenomic duplications. The second set included 141 taxa with confirmed the presence or absence of duplication. We applied MP and ML reconstruction methods. In the latter case, we used Mk1 model (Markov k-state 1-parameter model) because it fit the data better according to Akaike information criterion than the alternative AsymmMk model (asymmetrical Markov k-state 2-parameter model).

Analyses of DNA Asymmetry, Morphological Features, and Life Span

The DNA asymmetry was measured by normalized differences in complementary nucleotides, that is, AT-skew=(A − T)/(A + T) and GC-skew=(G − C)/(G + C), which results from different mutational patterns in the heavy (H) strand than the complementary light (L) strand of mitochondrial DNA during replication (Reyes et al. 1998). These parameters were calculated for the L strand sequence of 47 complete parrot genomes, 23 with the duplicated region and 24 with the single version. We also analyzed all PCG and 4FD sites in this way. The composition of 4FD sites in the ND6 gene, whose sense strand is located in the genomic H strand, was determined from the complementary L strand.

The data on body mass (M), total body length (L), wingspan (W) and tail length (T), as well as the maximum life span (S) of parrot species were collected from AnAge (http://genomics.senescence.info/species), Global Species (https://www.globalspecies.org) and Species360 (https://www.species360.org) as well as references (Juniper and Parr 1998; Brouwer et al. 2000) (supplementary table S9, Supplementary Material online). Body mass was expressed as an average of values reported in these resources. To eliminate the tail length contribution to the total body length, we calculated L − T as well as ratios of W/M, W/(L − T), T/M, and T/(L − T) to normalize the wingspan and tail length by measures associated with body size. As body mass is positively correlated with the maximum life span in birds (Skujina et al. 2016), we modeled it using the linear model ln Smax=a⋅ln M + b. The model was fitted using data for 178 parrot species with known body mass (M) and the maximum life span (Smax). We also carried out separate modeling for parrots with and without the duplicated region.

To check whether the analyzed variables are normally distributed, we applied the Shapiro–Wilk test. As the variables were not shown to be normally distributed, nonparametric tests were applied. The unpaired Wilcoxon–Mann–Whitney test was used in the comparison of two groups. The paired version of this test was applied in the comparison of the length between two duplicated CRs within a given avian order, while the Dunn’s test following the Kruskal–Wallis test was used in the multiple comparisons of the CR lengths between the avian orders. In this multiple testing, P-values were corrected using the Benjamini and Hochberg method. The likelihood ratio (G-test) goodness-of-fit test with the conserved Yates correction was applied to verify whether the bias in the number of second CRs longer and shorter than the first ones deviates from the equal distribution. The null hypotheses were rejected at the 0.05 level. The statistical analyses were carried out in R package 3.3.3 (R Core Team 2017).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.