Mitogenomic Architecture
The mitochondrial genome of the Andean and California condor was 16.808 and 16.870bp in length, respectively (GenBank accession no. xxxx), containing all 37 genes typical of vertebrates, including 22 tRNAs, 2 rRNAs, 13 PCGs, and a control region (Fig. 1). As previously observed in partial mtDNA sequences [35], NW vultures exhibited the ancestral avian gene order, while OW species differed from the standard arrangement due to the presence of a pseudo control region [14;46–48]. The overall base composition of condors was similar (A = ~ 31%, T = ~ 23%, C = ~ 32%, G = ~ 13%), showing a slight A + T bias for most regions, except ND6 that showed a marked C + G bias, similar (r > 0.9; P < 0.001 in all cases) to other vultures (Figure S1). Our comparative identity analysis revealed a clear phylogenetic pattern for both PCGs and rRNAs, except for COX1 and COX3 that showed a relatively high similarity between OW and NW vultures, suggesting an important evolutionary constraint acting on these regions. Most tRNA sequences were highly conserved across species, except for trnR, S1, E, F, N and C. The latter two showed an unusually high divergence, even within OW vultures, suggesting a possible pattern of selective relaxation (Fig. 2a). The tRNAs of condors displayed the typical cloverleaf of secondary structures found in other vultures, exhibiting larger structural differences with OW vultures (Figure S2).
Our analysis of PCGs in condors showed that COX3, ND4 and CYTB, share the same transcriptional exceptions of Turkey vultures, where TAA stop codon is completed by the addition of two 3' A residues to the mRNA. This exception has been reported for COX3, ND2 and ND4 in the Griffon (Gyps fulvus) and Cape (Gyps coprotheres) vultures [47;48], and for ND2 in the Cinereous vulture and Himalayan griffon [14;46]. In addition, we found a frameshift in the ND3 gene of both condors, suggesting an alternative reading of mitochondrial mRNA. Previous studies have found that this feature is present in some birds, including some Cathartes species, the Cape and Griffon vultures, but seems to be absent in the Himalayan Griffon and Cinereous vultures, suggesting the existence of additional evolutionary constraints in the former taxa due to the maintenance of a secondary structure in the mRNA [49]. Overall, we found that the relative amino-acid and synonymous codon usage in condors was similar to other vultures, exhibiting a frequent use of Leu, Thr, Ala, Ile, Gly, Ser, and Pro, with a common codon usage of CUA, ACA, GCC, AUC, GGA, UCA, and CCA, respectively (Figure S3).
The length of the CR of the Andean condor (1,214 bp) was similar to the California condor (1,260), Turkey vulture (1,177 bp) and Gyps species (1,201–1,206 bp; A. monachus was excluded due to the poorly resolved CR). Our finding of high similarity of the Andean condor CR with that of the California condor (90%), Turkey vulture (88%) and Gyps species (~ 68%), similar to the observed in PCGs (Fig. 2a), is remarkable given that the CR of birds tends to evolve faster than PCGs [50]. We detected common conserved regions in all three domains, following the same order in all vultures (Fig. 2b; Table S2). In the first domain, we found a Poly-C region usually observed in raptors [14], but exhibiting a TA island in the Andean condor, California condor, Turkey vulture and the Griffon vulture, and we also identified the adjacent extended termination-associated sequence ETAS1. The central domain contained the short conserved blocks F, E, D, C and conserved sequence blocks CSBa, CSBb and the avian-specific Bird Box, while the third domain included CSB1 (Fig. 2b). However, we did not detect the presence of ETAS2, CSB2 or CSB3 in any species, suggesting that these sequences are not conserved among vultures. We detected a defined Short Tandem Repeat sequence (motif: AACAAAC) at the 3’ end of the third domain (Fig. 2b), in the Turkey vulture (14 repeats), the California (24 repeats) and Andean condor (20 repeats), previously reported as a source of length heteroplasmy in the latter [35]. Notably, this pattern is not present in OW vultures, where repeats are rich in T, and are located in the pseudo control region [14;46;47]. The STRs of the Turkey vulture was able to form a stable stem-loop secondary structure, similar to those reported in some Strigiforms [51], while the stem of both condors was found in the terminal region of the third domain adjacent to the STRs (Figure S4).
As expected, the first and third domains showed higher intergeneric variability between condors, than the second domain (Fig. 2b). Moreover, the intraspecific variation pattern reported in California condors [18](Fig. 2b) is consistent with the expectation that regions exhibiting the greatest interspecific variability tend to vary the most within species [50]. However, we found that most of the intraspecific variability reported in the Andean condor [2;35] was situated between the Bird Box and CSB1, while the third domain showed an unusually invariant pattern (Fig. 2b), suggesting the existence of functional constraints on this region. Our observation of high similarity between the CR of our NW vultures, with a divergent time as old as ~ 16 My (Fig. 3), along with a similar substitution pattern to those of PCGs, and the observation that STRs could be involved in secondary structures, suggest that the extremely low variability found in the Andean condor could be possibly explained by a slow evolutionary rate, especially in the third domain. Slow rates of evolution of the CR has been previously reported in some Galliforms, Gruiforms, Passeriforms and Charadriiforms [1;50], while STRs found in the third domain of owls and gulls suggested that these sequences might be functionally essential in terminating mitochondrial genome replication [1;51]. Thus, we suggest that as in the case of the California condor, targeting the first domain of the Andean condor (between ETAS1 and F-BOX), might provide higher intraspecific variability for population genetic studies, especially when using degraded samples for ancient DNA approaches.
Phylogenetic Position And Divergence Times
Our phylogenetic analyses were consistent across sampling schemes and between ML and BI methods, recovering the same phylogenetic positions for OW and NW vultures, except in the BI analysis of scheme 3 (all PCGs + nonPCGs), that failed to group NW vultures with Accipitriformes (Figure S5). Nodal supports obtained with PCGs (scheme 1) increased with the addition of rRNAs (scheme 2), but decreased with the inclusion of tRNAs (scheme 3), probably due to the conservative nature of these short sequences (Fig. 1a). Thus, we estimated the maximum clade credibility chronogram using 13 PCGs + 2 rRNAs. Overall, our results were similar to recent analysis using whole nuclear genome sequences of raptors [10;11], splitting our in-group into two basal branches: Afroaves, including Accipitridae (eagles, hawks, kites, harriers, buzzards and OW vultures), Pandionidae (ospreys), Sagitaridae (secretary birds), Cathartidae (NW vultures), Strigidae (owls); and Australaves, including Falconidae (falcons and caracaras; Fig. 3). Hence, the inclusion of mitogenomic data of condors supports the hypothesis that NW and OW vultures belong to Accipitriformes [10;52] and that vultures are not closely related to Falconiforms as previously believed.
Among the Accipitriformes, Cathartidae was the first family to split off, suggesting that the striking similarities between NW and OW vultures could be regarded as an exemplar case of convergent evolution driven by the selective pressures of the obligate scavenger lifestyle. Some evident distinctions between these groups, such as the functional hind toe, lack of syrinx, internal separation of the nostrils, or the lack of squirting behavior on the legs in OW vultures, among others [7;9], suggest that evolutionary constraints (parallel evolution sensu neo-Gouldian [53]) are not driving their overall phenotypic similarities. Moreover, NW vultures seem to be more sensitive to lead contamination than OW species [20], while the striking tolerance to Diclofenac in the Turkey vulture (> 100 times than in OW vultures) suggest that NW vultures are less vulnerable to the toxic effects of non-steroidal anti-inflammatory drugs [54]. If true, components of genetic variation of the detoxification metabolism have diverged between NW and OW vultures. Future comparative analysis of the nuclear genome between these groups would help to understand to what extent NW and OW vultures exhibit divergent genetic pathways. This information will be critical not only to improve diagnosis and treatments, but also to assist in management plans aimed at maximizing adaptive genetic variability of threatened populations.
Our molecular dating of the origin of NW vultures in the early Paleogene (58.8 Mya; CI: 67.4–50.4; Fig. 3) differs from the estimates of recent analysis using whole nuclear genome data of raptor species [11], but is in agreement with previous studies including wider taxon sampling [10;52]. A likely explanation for this discrepancy is that Zhou et al.,[11] had a sole representative taxon of NW vultures (Turkey vulture) and lacked genomic resources for the extant representatives of Pandionidae (osprey) and Sagitaridae (secretary bird), resulting in a delayed estimate of the origin of Cathartidae (36 Mya). In fact, the fossil record indicates that Sagitarid, Pandionid and Accipitrid lineages split before the Oligocene, with the oldest putative representatives Pelargopappus schlosseri, Palaeocircus cuvieri, and Milvoides kempi, occurring in the early Oligocene, and mid-late Eocene, respectively [55;56]. Moreover, our estimated divergence time of NW vultures was robust to the exclusion of the calibration node using the controversial fossil P. cuvieri (results not shown), and agrees with the age of the oldest known fossils of the stem group of Cathartidae, Diatropornis ellioti and Paracathartes howardae, around the early-mid Eocene, almost 40 to 55 Mya, respectively [56;57]. In addition, our determination of the diversification time of extant Cathartids between the mid and late Miocene (10–16 Mya; Fig. 3), agrees with idea that small vultures (Cathartes) are more primitive than condors and is consistent with the earliest record of the condor-like vulture Hadrogyps aigialeus from the Barstovian period (13–16 Mya) of California [58]. This result supports the early arising of this group, probably initiated by the paleo-climatic conditions and diversification of megafauna during the Miocene [59;60]. The historical record shows evidence that NW vultures were highly diverse during that period and harbored large populations until the early Holocene when the climate became warmer and humid, resulting in weaker thermal updrafts, limiting soaring flight, which together with the loss of megafauna has likely contributed to their demise [58;59;60]. The oldest undisputed records of extant Andean condors from Pleistocene deposits found in Peru and Chile [61;62] suggest the recent origin of the species in South America during the rapid rise of the Andes. Notwithstanding, previous studies proposed that condors evolved from a smaller species in North America and expanded to the south crossing the ancient coastline (before the Isthmus of Panama), during the late Miocene - early Pliocene [61;63] (but see [62]). If true, the characteristic large body size of California and Andean condors may represent a case of rapid convergent evolution driven by the selection pressures of the obligate scavenger lifestyle in mountainous environments. Although more efforts are still needed to elucidate the natural history of vultures, the use of mitogenome data in conjunction with an increased taxon sampling seems a promising approach to untangle the long-standing controversial relationships within this group.
Selective Relaxation
Our selection analyses suggest that vultures experienced similar levels of genetic constraints to the rest of accipitrimorphs, except for the selective relaxation detected in ND1 and ND4 codons in the former group (Fig. 4 and S6; Table S3). A possible explanation for these differences might be related to the contrasting flight abilities of these raptors. Our finding of selective relaxation in those genes is in line with previous studies that found larger substitution rates in ND genes in weak versus strong locomotive birds due to the stronger selective pressure to eliminate deleterious mutations in the latter [64]. Thus, it is possible that the stronger purifying selection observed in non-vulture accipitrids is related to the active metabolism of flapping flight, while the specialization of soaring flight in vultures resulted in the relaxation of selective constraints on the energy metabolism [23;24]. In fact, recent reports demonstrated that vultures are able to fly with a strikingly minimal amount of flapping [65], which in the Andean condor can be as low as 1% of the total flight time, a record for any living flying bird [25]. In addition to the locomotive strategy, our results also suggested the implication of phylogenetic constraints influencing mitochondrial evolution. NW vultures showed unique and radical non-synonymous substitutions in several genes compared to the rest of Accipitrimorphs (Figure S7). Furthermore, our ω estimates for NW species were smaller, suggesting stronger purifying selection (Figure S6), while our CmC analysis showed that OW vultures exhibited significant relaxation in ND3 and ND6 codons (Fig. 4; Table S3). This is consistent with our observation of large numbers of non-synonymous substitutions fixed in these genes of OW vultures (Figure S7), and lower ω values in the terminal branches of NW species compared to OW vultures (0.04 and 0.28, respectively; P < 0.005 in both cases). Notwithstanding, a larger number of species is needed to confirm the phylogenetic inertia of substitution patterns.
Our CmC analysis of selective signals according to the body size of vultures, revealed that larger species exhibited significant relaxation in ND6 codons, while smaller species showed selective relaxation in ND1 codons (Fig. 4). However, our regression analyses showed no relation between the body mass of vultures and their ω values in any gene (P > 0.22 in all cases). Furthermore, after dis-aggregating the CmC analysis within OW and NW vultures, our results suggested that the observed relaxation pattern wight be explained by the strong relaxation signal of large OW species, especially in ND6 (Fig. 4, Table S3). This was further confirmed by our BM and BSM analysis (P < 0.001), suggesting the involvement of positive selection acting on ND6 sites (Table S3). Thus, our results propose that the evolution of body size in this guild is not associated with a common or simple mitochondrial genetic basis. Given that large OW species inhabit higher mountains and fly at higher elevations than condors, our detection of positively selected ND6 sites suggests a possible signal of mitochondrial adaptation to hypoxia. Previous studies have shown that positive selection of the NADH dehydrogenase complex plays an important role in adaptation to high-altitude in mammals species and in galliform and passeriform birds [6;66]. If true, our observed selection footprint should be also present in smaller high-altitude vultures. Thus, the inclusion of the mitogenome of the Ruppell’s griffon (Gyps rueppelli) will be instrumental to test this hypothesis, as this medium-sized vulture (~ 7 kg) can fly as high as 11.300 m [31]. Finally, our CMC and BM analyses focused on NW vultures revealed a selective relaxation signal in condor´s CYTB gene, respect the Turkey vulture´s ortholog (P = 0.028; Table S3). This result could be related to the migratory behaviour of Turkey vultures [67], which might exert higher selective pressures on mt-genes related to the energy metabolism of this small vulture.
Overall, our results suggest that the exceptional specialization of soaring flight of vultures has compensated the evolution of body mass without pressing directional changes on mtDNA. The energetic cost of flapping flight is expected to be extremely high compared to soaring or gliding flight, and lighter species tend to flap more than large ones [25]. Given that vulture species are among the largest flying birds and evolved the most extreme use of soaring flight [23–25;67], the relaxation of selective constraints on mt-genes may indicate the degeneration of flapping flight ability. Our results are consistent with previous studies that found small differences in the energy consumption between flight and resting in large vultures [25;65], suggesting that low metabolic rates might have buffered the selection pressure on mitochondrial efficiency. Notwithstanding, nuclear genomic sequences and an increased number of species are needed to better understand how the evolution of the obligate scavenger lifestyle affected metabolic selection.