Time-calibrated phylogeny and full mitogenome sequence of the Galapagos sea lion (Zalophus wollebaeki) from scat DNA

Introduction

The family Otariidae is a main clade of pinnipeds (Berta & Churchill, 2012; Berta, 2018), including seven genera and 14 species. It is morphologically characterized by an auditory pinna, and posterior limbs directed forward (Wursig, Perrin & Thewissen, 2009; Berta & Churchill, 2012). The distribution of otariids comprises the coasts of North and South America, the coasts of North Asia, Australia, New Zealand, and in addition several island systems including the Galapagos archipelago (Wursig, Perrin & Thewissen, 2009). The Galapagos sea lion, Zalophus wollebaeki (Sivertsen, 1953), is an endemic and endangered otariid species (Trillmich, 2015), which concentrates its largest populations in the southeastern part of the archipelago (Riofrío-Lazo, Arreguín-Sánchez & Páez-Rosas, 2017; Páez-Rosas et al., 2021). Given its position as top predator and its high nutritional efficiency, it is considered as a sentinel species of the ecosystem (Páez-Rosas & Guevara, 2017), and its permanence in the region is assumed to regulate the health and functioning of the Galapagos marine ecosystem at all trophic levels (Fariña et al., 2003; Riofrío-Lazo et al., 2021). However, it remains exposed to periodic oceanographic-atmospheric events in the region, such as warm phases of the El Niño–Southern Oscillation, which has caused drastic population decreases in the last four decades (Páez-Rosas et al., 2021).

Despite the advent of high-throughput genomic approaches, mitochondrial DNA remains an important tool in phylogenetic and population genetic inference (DeSalle & Hadrys, 2017), and has also been used for population genetics of Z. wollebaeki (Arnason et al., 2006). In animals, due to its high copy numbers compared to nuclear DNA, it is a particularly suited marker to be obtained from low-quality DNA samples (DeSalle & Hadrys, 2017). For the genus Zalophus, complete mitogenome sequences have been published both for the California sea lion, Z. californianus, and the extinct Japanese sea lion, Z. japonicus (Arnason et al., 2006; Kim et al., 2021). Here, we complement the mitogenomic resources for Zalophus genus, by a complete mitogenome, sequenced from DNA extracted from a scat sample of Z. wollebaeki. In the same way, we constructed a molecular clock using the complete mitogenomes of otariid species.

Materials and Methods

A fresh scat sample was obtained from an Zalophus wollebaeki individual on the El Malecón rookery at San Cristóbal island, Galapagos archipelago (Lat: −0.90072, Long: −89.610117) in July 2021 as part of a trophic ecology project. The sample was stored in the Molecular Biology Laboratory of the Galapagos Science Center (GSC) in a −20 °C freezer. Total DNA was extracted with the QIAamp Fast DNA Stool Mini Kit (Qiagen, Hilden, Germany) from 220 mg of sample following manufacturer instructions with some modifications such as decreased incubation temperature from 70 to 50 °C, and the final elution of extracted DNA in 100 µl of TAE buffer. Agarose gel electrophoresis (1%) was performed to check purity and integrity of DNA, while a NanoDrop 2000 spectrophotometer was used to check quantity, yielding 56.1 ng/µl of DNA.

A library of the total extracted DNA was prepared using an Illumina TruSeq stranded protocol and sequenced together with 12 unrelated (RNAseq) samples on an Illumina NextSeq system using a 2 × 150 cycles using a High Output kit. The obtained raw reads (3,511,116) were processed using the GetOrganelle software (Jin et al., 2020) to filter the mitochondrial Zalophus DNA reads (remaining ∼16% of raw reads) and assemble them. GetOrganelle was chosen because it deals appropriately with contamination, an evident issue with the scat sample we work with (Freudenthal et al., 2019). Stool contains DNA from other sources, including mainly bacteria but also some prey, parasites, and environmental contaminants. GetOrganelle firstly filters read mapping to a mitochondrion reference, discarding most of the contaminant reads in this step (Jin et al., 2020). Then, it employs the SPAdes assembler that is designed for single-cell sequencing, so it allows to discriminate sequences from different sources with different coverages (Bankevich et al., 2012). GetOrganelle then recognizes the target mitogenome based on connections, coverages, and BLAST hits, putting aside any remaining contamination sequences (Jin et al., 2020).

In order to adjust the reading frame of the obtained assembly and to facilitate its alignment with published mitogenomes, its reverse complement was obtained, and the circularization cut-off was set in the position 373; in this way we succeeded obtaining a final quasi-mitogenome assembly for Z. wollebaeki. Quality control of the final mitogenome assembly was assessed by getting a complete annotation of it using the MITOS Web server 2.0 (Donath et al., 2019), and by comparing the nucleotide similarity (via NCBI Megablast) of this mitogenome with the references of Z. japonicus, Z. californianus and Eumetopias jubatus (Steller’s sea lion). Syntheny across the three Zalophus species was assessed by comparing their mitogenomes’ annotations; a codon usage analysis of the three species was also performed employing the Sequence Manipulation Suite (Stothard, 2000). Finally, the Z. wollebaeki mitogenome map was constructed using OGDRAW 1.3.1 (Greiner, Lehwark & Bock, 2019).

A time-calibrated phylogenetic tree was constructed including the newly assembled Z. wollebaeki mitogenome together with other complete otariid mitogenomes published at GenBank (Table 1). The harbor seal (Phoca vitulina), grey seal (Halichoerus grypus), and walrus (Odobenus rosmarus) mitogenomes were set as outgroups. All the mitogenomes were aligned in MAFFT online (Katoh, Rozewicki & Yamada, 2019). Then, JModelTest 2.1.10 (Darriba et al., 2012) was employed to choose the most appropriate evolutionary model, using all the best model selection tests (i.e. AIC, BIC and DT). All tests indicated that the best model for the data set was GTR+I+G. With the information from JModelTest it was possible to adjust the initial parameters for Gamma, portion of invariant sites, nucleotide frequencies and nucleotide exchange rates that were included in the subsequent Bayesian analysis to determine divergence ages. MrBayes 3.2.7a (Ronquist et al., 2012) was used for Bayesian and coalescent analysis, fitting the parameters retrieved from JModelTest, and including the calibrations (priors) from Table 2. A MCMC of 3,000,000 generations were run (with a burn-in of 25%). Sampling was taken every 100 records and the molecular clock was set to strict. The resulting parameters of interest showed outstanding convergence. The time scale was adjusted in the tree using 2% divergence. The resulting tree was visualized in FigTree version 1.4.4 (Rambaut, 2010).

Results

We obtained a final Zalophus wollebaeki mitogenome assembly of 16,676 bp in length (GenBank: OP115880), with a GC content of 40.5%. The average coverage of the mitogenome was 55×. The newly assembled Z. wollebaeki mitogenome consists of 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNA), and two ribosomal RNAs (rRNA) (Fig. 1), agreeing that seen in the other Zalophus mitogenomes. Some gaps (labelled as Ns in our mitogenome) were observed in the non-coding control region (D-loop) when aligning the Z. wollebaeki mitogenome with the published ones, these being potential segments that couldn’t be successfully assembled. However, the total length of these gaps did not exceed the 100 bp total.

Sequence comparison revealed that Z. californianus was 99.42% similar to our Z. wollebaeki mitogenome (16,346/16,439 bp were identical), with only one gap (1 bp indel) found. In the same way, Z. japonicus showed 98.65% of similarity (16,215/16,440 bp were identical), with four gaps found. On the other hand, the Eumetopias jubatus mitogenome displayed a similarity of 94.11% with respect to Z. wollebaeki.

The codifying elements at the mitogenomes of the three Zalophus species were perfectly collinear and no major indels, inversions or translocations were found when comparing them in syntheny analysis. Codon usage analysis revealed that Z. wollebaeki and Z. californianus are biased towards the same codon at 16/20 aminoacids plus the stop codon. When comparing Z. wollebaeki and Z. japonicus, they are biased for the same codon at only 13/20 aminoacids plus the stop codon (Table 3).

The time-calibrated phylogenetic tree (Fig. 2) showed a monophyletic group containing the three species of Zalophus, supporting the sister group relationship of Z. californianus and Z. wollebaeki with a divergence time estimation of 0.65 mya. It also shows the genus Zalophus diverged from E. jubatus 7.98 mya. Our phylogeny supports Eumetopias being sister to the genus Zalophus, while the clade containing these two genera would be sister to the remaining otariids (Southern hemisphere otariids). The divergence times confidence intervals are mostly quite tight, demonstrating the high degree of convergence of the substitution rates and thus of the node ages.

Discussion

The general features of Zalophus wollebaeki mitogenome—such as length, GC content, and PCGs, rRNAs, and tRNAs found—were consistent to those of the other two Zalophus mitogenomes: Z. japonicus (16,698 bp; Kim et al., 2021) and Z. californianus (16,677 bp; Arnason et al., 2006). It must be noted that Z. wollebaeki was sequenced with the Illumina method, while Z. japonicus and Z. californianus were sequenced through Sanger (Arnason et al., 2006; Kim et al., 2021). All three sequences presented a high similitude regardless of the sequencing method used. Confirming that any type of sequencing can recovered a complete mitogenome.

Our phylogenetic reconstruction supports the genus Zalophus being monophyletic, with Z. californianus being the sister species and hence a possible antecessor of Z. wollebaeki (Wolf, Tautz & Trillmich, 2007; Wolf et al., 2008). The perfect collinearity and high sequence similarity (>98%) across these three species further support monophyly in Zalophus. Although a recent phylogenomics study failed in supporting the monophyly of the Zalophus genus (Lopes et al., 2021), most of the body of evidence in otariid phylogeny supports consistently this pattern (e.g., Yonezawa, Kohno & Hasegawa, 2009; Churchill, Boessenecker & Clementz, 2014). Our results support that Eumetopias jubatus is sister to the genus Zalophus, while Eumetopias and Zalophus together would be sister to the remaining otariids. This also coincides with other phylogenetic studies in Otariids, from those based on mitochondrial sequences only (Yonezawa, Kohno & Hasegawa, 2009; Churchill, Boessenecker & Clementz, 2014; Austin, Eriksen & Bachmann, 2023) to those based on phylogenomics (Lopes et al., 2021).

Our phylogeny, high nucleotide similarity and close codon usage patterns between Z. wollebaeki and Z. californianus supports your relationship of sister species (Athey et al., 2017). However, enough nucleotide differences (including 16 substitutions at the D-loop region) would exist for supporting Z. wollebaeki being a separated species from Z. californianus, as also suggested from the Z. wollebaeki paratype mitogenome analysis (Austin, Eriksen & Bachmann, 2023). According to our estimations, the Z. wollebaeki and Z. californianus split would be dated 0.65 mya, a more recent date than the 2.5 mya estimated by Wolf, Tautz & Trillmich (2007). Such a recent divergence between both species agrees with the evident morphological similarity they share, being size the only distinguishable difference between them (Z. wollebaeki is smaller) (Sivertsen, 1953; Orr & Scheffer, 1958); reason why until recently, Z. wollebaeki was classified as a subspecies of Z. californianus (Wolf, Kauermann & Trillmich, 2005; Wolf, Tautz & Trillmich, 2007). Our estimation of 0.65 mya also matches with the Kansas glaciation, an especially strong glacial period occurred 0.7−0.6 mya (Aber, 1991). Global cooling events as this one increased ocean productivity even in equatorial waters, facilitating the spread of pinnipeds across hemispheres and tropical environments (Churchill, Boessenecker & Clementz, 2014; Lopes et al., 2021). Therefore, the arrival of the ancestors of Z. wollebaeki into the Galapagos Islands and its divergence from Z. californianus around this time, is a perfectly plausible event.

We successfully annotated in Z. wollebaeki mitogenome all the genes expected in every mammal mitogenome. Thus, our study demonstrates the prospect of recovering complete and accurate mitogenomes of the host from scat samples in the field. A similar method has previously been used for mitogenome sequencing of mammals, including artiodactyls (i.e., bharal Pseudois schaeferi) (Dong et al., 2019), and primates (i.e., toque macaque Macaca sinica) (Roos, 2018). In the latter case, DNA extracted from feces yielded complete mitogenomes by Sanger sequencing, with sufficient coverage to assemble it without inconvenience, leaving aside the rest of reads present in the feces (Roos, 2018). Note we achieved highly supported results, with accurate time estimates and no phylogenetic ambiguities, by analyzing entire mitogenomes without partitioning them for modeling. Thus, our work joins other phylogenetic studies among closely related species (e.g., Zhang et al., 2021; Zawal et al., 2022; Sun et al., 2023) where the use of whole mitogenomes instead of partitioning yields better results by preventing the addition of unnecessary complexity and bias in our analyses, while providing a greater amount of information for inference (Duchêne et al., 2011).

Despite the unspecific approach used herein (shotgun sequencing of total DNA, which mostly yielded DNA reads of bacteria), assembly was straightforward and did not require any filtering approach. This suggests that host mitochondrial DNA was common in the DNA extracted from the scat sample, and contaminants from other metazoans did not influence the assembly in a meaningful way. The assembly process in our study certainly was also facilitated by the existence of mitogenome sequences of related species of Zalophus. We obtained, with a fecal DNA extraction kit, a DNA concentration of 56.1 ng/µl which was sufficient for library preparation and assembly of a sequence of 16,676 bp, but we assume that with less fresh samples, DNA yield and proportion of predator mtDNA may be less favorable. So, we expect the novel method presented in this work to be applicable for fresh scat samples of other pinniped species.

Conclusions

We stress the novelty of our study, obtaining a complete mitogenome (excepting by a <100 bp stretch) for Zalophus wollebaeki from fresh scats collected in the field, in a non-invasive sampling. With this information we were able to compare distinct sequences and estimate the divergence time of Z. wollebaeki, which matched with a global cooling event. Additionally, we supported the sister-species relationship of Z. wollebaeki and Z. californianus, as well as the monophyly of the Zalophus genus.