Multi-locus evaluation of gastrointestinal bacterial communities from Zalophus californianus pups in the Gulf of California, México

Sample collection

A total of 54 rectal cotton swabs (RCS) were collected from healthy and randomly selected sea lion pups (2–3 months old) located in six rookeries (n = 9 FCS per rookery) in the Gulf of California (Fig. 1). The RCS samples were collected in compliance with Mexican wildlife regulations (SGPVA/DGVS/003083/18). Before sampling, all pups were manually restrained and anesthetized with isoflurane (5%) to reduce mobility and prevent pain. Anesthetic agents were administered by veterinarians of the African Safari Zoo (Puebla, Mexico). All collected RCS were preserved in liquid nitrogen until DNA extraction.

DNA extraction and purification

Due to the complexity of fecal samples, DNA was extracted and purified using both the Wizard Genomic DNA purification kit (Promega Corporation, Madison, WI, USA) and the PureLink Invitrogen Genomic DNA kit (Thermo Fisher Scientific, Waltham, MA, USA) with protocol modifications. Briefly, 600 μL of nuclei lysis solution was added to 1.5-mL microcentrifuge tubes containing the frozen FCS heads. The tubes were incubated for 25 min at 80 °C. Contaminating RNA was degraded by adding 3 μL of RNAse solution and incubating the tubes for 30 min at 37 °C. Protein was precipitated by adding 600 μL of protein precipitation solution and incubating the tubes on ice for 5 min. Precipitated organic matter was removed by centrifugation for 3 min at 16,000 rpm, and the supernatants were transferred to new, sterile tubes. The DNA was precipitated by adding 600 μL of isopropanol to the supernatants, followed by inversion mixing and centrifugation for 2 min at 16,000 rpm. The supernatants were discarded, and the pellets were washed by adding 600 μL of molecular grade ethanol. The tubes were centrifuged at 16,000 rpm for 2 min before the ethanol was carefully removed and then drained at room temperature for 10–15 min to remove all residual ethanol. The DNA was resuspended in 50 μL of nuclease free water (NFW) and stored overnight at −20 °C. Final DNA purification was achieved following manufacturer protocols with 700 μL of wash buffer 1; wash buffer 2 was used to wash the DNA. A NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to measure the DNA concentration and quality before visualization in agarose gel (1.5%).

Sequencing and reference library construction

Prior to sequencing, the DNA extracts from RCS collected in the same rookery were quantified, normalized, and pooled in equal amounts, resulting in 18 pooled DNA samples. We amplified these samples using an Ion 16S™ Metagenomics kit (Thermo Fisher Scientific, Waltham, MA, USA) following manufacturer protocols. Amplification involved 25 PCR cycles with Ion Xpress Barcoded adapters. Emulsion PCR employed OneTouch™ 2,400-bp chemistry. Sequencing was performed on an Ion Torrent PGM using an Ion 510 chip Kit (Thermo Fisher Scientific, Waltham, MA, USA) with 400-bp read lengths. Torrent Suite v. 4.2.1. Ion Reporter Metagenomics 16S software (v. 5.2) was used for early read trimming and processing according to the following quality criteria: minimum base-calling error of Q20, no barcode mismatches, ≤3 primer mismatches, and length ≥150 bp. Primer detection was set to single end (minimum 90% alignment coverage between hit and query sequences). The raw DNA libraries were deposited in the National Center for Biotechnology Information (NCBI) database (BioProject PRJNA715917).

Query library preparation and bioinformatic analyses

Raw reads were initially trimmed with Trimmomatic (v. 0.38; Bolger, Lohse & Usadel, 2014) with the following parameters: Leading 3; Trailing 3; Slidingwindow 4:15; and Minlen 150. The trimmed reads were then imported into Quantitative Insights Into Microbial Ecology (QIIME 2; v. 2020.6; Bolyen et al., 2019). Dada2 (Callahan et al., 2016) was employed for denoising and detecting chimeras before the reads were clustered into amplicon sequence variants (ASVs; 100% similarity). The ASVs from the same HRs were filtered by de-novo clustering using VSEARCH (Rognes et al., 2016) and trained on each 16S rDNA reference library in the Ion Reporter Software. The ASVs were aligned with MAFFT (v. 7; Katoh & Standley, 2013), and the phylogenetic reconstruction was carried out with fastTree (v. 2.1; Price, Dehal & Arkin, 2010). The taxonomic position of each ASV was obtained using the BLAST+ consensus taxonomy classifier (Camacho et al., 2009) trained on the SILVA database (v. 132; Quast et al., 2013). To minimize the inclusion of spurious ASVs in the final data sets, singletons and ASVs composed of less than a specific number of HR-library reads (ranging from 15 for V2 to 72 for V8) were removed. Each variable region were then rarefied at the minimum sample read depth using the QIIME 2 core-metrics-phylogenetic plugin, and the resulting abundance tables were used in downstream analyses.

Ecological and functional prediction analyses

Bacterial communities were initially evaluated via rarefaction curves with the median read frequency used as the sequencing depth. Both alpha- and beta-diversity were determined using the rarefied number of reads.

A permutational multivariate analysis of variance (PERMANOVA) was used to evaluate structural differences between rookeries with Monte Carlo tests and 4,999 permutations in PRIMER+P (v. 6; Clarke & Warwick, 2001) using Bray-Curtis, Jaccard, and weighted and unweighted UniFrac distances. However, the principal coordinate analysis (PCoA) results were visualized with EMPeror (Vázquez-Baeza et al., 2013) using a Jaccard distance of the most abundant ASV (>1% of the read frequency). Finally, to evaluate the principal bacterial taxa driving the dissimilarities, a similarity percentage (SIMPER) analysis of the ASVs detected within V6–V7 was performed in PRIMER+P (v. 6; Clarke & Warwick, 2001).

PICRUSt (v. 2.2.0) beta pipelines (Douglas et al., 2019) were used to predict the potential functional capabilities of the ASVs as well as their potential contributions to enriching the predicted ecological functions among all HRs. The abundances of the detected ASVs were normalized by the predicted 16S rRNA gene copy number present in the closest reference genome. Subsequently, functional predictions were obtained with the predicted, orthologous genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database of each ASV. For quality control purposes, the nearest sequenced taxon index (NSTI) was determined. This index reflects the accuracy of the functional predictions for each ASV by conveying the average genetic distance (measured as the number of substitutions per site) between an ASV and the 16S rRNA of the most similar reference genome (Langille et al., 2013; de Voogd et al., 2015; Douglas et al., 2019). Following the suggested guidelines (Douglas et al., 2019), ASVs with NSTI values >2 were removed. To visualize the functional dispersion present among bacterial communities, PCoA plots were created with STAMP (v. 2.1.3; Parks et al., 2014) based on fourth-root normalized metabolic KEGG pathway counts. Significant functional differences between rookeries were evaluated by PERMANOVA with Monte Carlo tests based on 4,999 permutations in PRIMER+P (v. 6; Clarke & Warwick, 2001).

Sequencing results

The yield of each HR was evaluated by read counts, and microbial composition and diversity were determined. A total of 281,630 raw reads (254-bp average length) were sequenced, with 266,177 (93.4%) reads meeting quality control test standards and being assigned to 1,852 ASVs. The number of detected ASVs and reads for each HR (read number reported in parentheses) were 107 (55,861), 87 (45,981), 53 (71,467), 114 (87,249), 7 (89,181), and 17 (33,398) for V2, V3, V4, V6–V7, V8, and V9, respectively. The rarefaction curves were asymptotic in shape, suggesting that the major fraction of the bacterial communities were sampled for V2, V3, V4, and V6–V7; however, the ASVs detected for V8 and V9 were mostly assigned to Proteobacteria, suggesting that the bacterial community was not entirely characterized with those regions (Fig. S1).

Bacterial structural compositions among HRs

The capacity of each HR to determine the composition of the fecal microbiome at different taxonomic levels was analyzed. At the phylum level (average read frequency in parentheses) and among all HRs, Proteobacteria (up to 100%) predominated followed by Bacteroidetes (up to 30.26%) and Firmicutes (up to 14.73%). Nevertheless, distinctive arrangements were observed among HRs (Fig. 2). For example, Fusobacteria were mostly detected with V6–V7, while Deferribacteres and Gemmatimonadetes were exclusively detected with V2 and V4, respectively (Fig. S2). Similar results were also observed at lower taxonomic levels. Distinctive abundance rearrangements of the 29 predominant bacterial orders (>1% abundance cutoff; Fig. 3) were observed for each tested HR. Given the average values of all ecological indices considered in this study, the highest bacterial community diversity was generally detected with V6–V7, whereas the least complex bacterial community was detected with V8 (Table 1).

Bacterial compositions between rookeries

The predominant orders among all rookeries were Enterobacteriales and Pasteurellales (Fig. 3); however, structural differences were also observed between sampling areas. Specifically, the SIMPER analysis outcomes (Table S1 from A to O) suggested that the major fraction of geographic variability was explained by shifts in the relative abundance of common ASVs mainly assigned to Proteobacteria, Bacteroidetes, and Fusobacteria phyla rather than the presence of exclusive bacterial ASVs among rookeries. Significant differences in bacterial communities were observed between all distance matrices (PERMANOVA; pseudo-f: >2.33; p < 0.0158) for V3, V4, and V6–V7. Moreover, within the V2 region, no significant differences were observed using Bray-Curtis distances between Coloradito and other regions (PERMANOVA; pseudo-f: <0.19; p > 0.072) nor were significant differences observed with the V8 and V9 regions (PERMANOVA; pseudo-f: <1.69; p > 0.11). The results of the PCoA based on Jaccard distance are reported in Fig. 4.

Despite not observing unambiguous geographic patterns, higher bacterial α-diversity was generally detected in Cantiles and Partido compared to the values observed in Coloradito, Granito, and Machos (Table 1).

Functional microbiome capabilities

Functional predictions were estimated by using 7 (V8) to 110 (V6–V7) ASVs. During these analyses, two, three, and four ASVs were removed from the V2, V4, and V6–V7 libraries, respectively, as these ASVs presented NSTI values >2. For the remaining 16S rRNA regions, all detected ASVs presented NSTI values <2. A higher number of functional pathways was observed for the V3, V6–V7, and V4 regions for which 147, 142, and 139 pathways were predicted, respectively (Tables S2–S7). Overall, functional capabilities (gene count percentages in parentheses) were mainly assigned to the metabolic pathways (up to 73.44%) of carbohydrate metabolism (14.92%), cofactor and vitamin metabolism (12.44%), and amino acid metabolism (13.02%). The major bacterial families responsible for assigning these functions were Enterobacteriaceae (up to 94.67%) and Pasteurellaceae (up to 18.7%). Additional bacterial families with relatively high functional contributions (threshold > 3%) in at least one HR are reported in Table 2.

Similar abundances in functional pathway profiles were observed among rookeries (KEGG hierarchical level 1; Fig. 5). Notably, seven bacterial genera were associated with the assigned metabolic activities and included Escherichia-Shigella (~22%), Fusobacterium (~16.22%), Otariodibacter (~14%), Bacteroides (11.88%), Ornithobacterium (6.77%), Porphyromonas (5.91%), and Pseudomonas (5.5%). Moreover, the metabolic functions in Coloradito and Granito were mostly enriched by Otariodibacter and Ornithobacterium, whereas Escherichia-Shigella and Bacteroides mainly contributed to the same functions in Macho, Partido, Rasito, and Cantiles, as did Fusobacterium in Macho and Rasito. The contributions of Pseudomonas and Porphyromonas were consistent among all rookeries. The PCoA results based on the normalized count numbers of the pathways at KEGG hierarchical level 3 suggested functional splits between most of the sampled rookeries for the V2, V3, V4, and V6–V7 regions although not for V8 and V9 (Fig. 6). Nevertheless, the PERMANOVA results revealed consistent differences only between Granito and Partido and Granito and Rasido for the V4 region (pseudo-f: >2.26; p < 0.0482).

Bacterial taxa abundance in pinnipeds

Technological advances in high throughput sequencing have provided unprecedented opportunities to characterize bacterial communities from almost any type of environmental DNA sample, including samples from the GI tracts of host species. However, accurately characterizing these communities can be challenging due to the unequal amplification of bacterial taxa during PCR (Brooks et al., 2015). Moreover, the main limitation of metabarcoding surveys can be the incompleteness of reference databases (Martinez-Porchas et al., 2017; Gold et al., 2021). The latter is particularly evident for external HRs, which are known to have incomplete reference sequences (Martinez-Porchas et al., 2017). Moreover, metabarcoding surveys measure relative rather than absolute biological abundance (Barlow, Bogatyrev & Ismagilov, 2020). To circumvent or at least minimize these difficulties in this study, we integrated the information from six HRs of 16S rRNA to obtain high-resolution microbial community profiles.

The combined information from the tested HRs suggests that the bacterial communities in the rectal samples of California sea lion pups are mainly composed of five phyla: Proteobacteria, Bacteroidetes, Epsilonbacteraeota, Firmicutes, and Fusobacteria. Despite most of these taxa being common and important for pinniped health, variations in their abundance remains a topic in need of discussion. For example, the findings from this study suggest that Firmicutes may cumulatively constitute up to 30% of GI bacterial communities in sea lions; however, Delport et al. (2016) reported that Firmicutes may constitute approximately 76% of the GI bacterial communities in Australian sea lions (Neophoca cinerea). Comparisons between studies should be made cautiously due to differences in methodology (e.g., variations in DNA extraction methods, PCR-targeted 16S rRNA regions, and downstream bioinformatic methods) and other factors (e.g., environment, health, inheritance, or age). Nevertheless, several studies that have used the universal primers 27F (5′-AGAGTTTGATCCTGGCTCAG-3′; e.g., Ludwig, Mittenhuber & Friedrich, 1993) and 519R (5′-GWATTACCGCGGCKGCTG-3′; e.g., Ruff-Roberts, Kuenen & Ward, 1994) that span the V1–V3 regions have usually identify Firmicutes as the predominant bacterial species in both pinniped (Nelson et al., 2013; Acquarone, Salgado-Flores & Sundset, 2020) and non-pinniped host species (Cicala et al., 2018; Brown, Wiens & Salinas, 2019). Conversely, when alternative primer sets targeting V4 have been implemented, the reported Firmicutes abundances are more consistent with the findings of this study, as evidenced by Phoca vitulina richardii (~26%), Leptonychotes weddellii (~23%), and Arctocephalus pusillus (~25%) abundance (Banks, Cary & Hogg, 2014; Pacheco-Sandoval et al., 2019; Bai et al., 2021).

Discrepancies in bacterial taxa abundance can be extended to include Proteobacteria and Epsilonbacteraeota. To the best of our knowledge, Epsilonbacteraeota is one of the predominant bacterial taxa in the GI microbiota of Cape fur seals within V4 (Bai et al., 2021); however, this bacterial group has not been previously identified in the GI bacterial community of Zalophus californianus nor those of other sea lion species. Moreover, in the present study, Epsilonbacteraeota was largely detected within the V3 region. Together, these findings suggest that Epsilonbacteraeota may be considered a predominant bacterial group in sea lions, even though we speculate that this bacterial taxon has not been previously detected in sea lions given that V3 remains a poorly utilized region in metabarcoding surveys. Nevertheless, more evidence is needed to support or disprove this hypothesis. Our results also identified Proteobacteria as one of the main taxa in bacterial communities, with abundances as high as 60% in the fecal samples of sea lion pups. Conversely, several authors have reported abundances of this phylum as low as 1–2% based on the genetic information from V1–V2 or V3–V5 (Smith et al., 2013; Bik et al., 2016). Thus, further research is required to determine the exact proportion of Proteobacteria in pinniped GI tracts.

Studies examining GI bacterial diversity using six HRs in sea lions or other pinniped species are lacking; thus, our results may encourage a reevaluation of previously reported bacterial abundance profiles in pinnipeds. Moreover, the proposed observations may provide starting points for further investigations, as the most comprehensive description of the bacterial community was obtained by integrating information from the V3 and V6–V7 regions. These HRs may reflect GI changes better than the other 16S rRNA regions (Vargas-Albores et al., 2017); hence, future bacterial surveys of sea lions, and possibly those of other pinniped species, may benefit from sequencing approaches employing the V3 and V6–V7 regions.

Bacterial community variation between rookeries

Given that distinctive HR-related bacterial consortia were detected, we evaluated which among the tested HRs might be most suitable for identifying structurally related microbiome features, such as geographic changes or functional capabilities (reported in the next paragraph). Even after considering the low number of samples processed for this analysis (n = 3 per rockery), consistent structural variation among rookeries was observed among V3, V4, and V6–V7, whereas marginal geographic changes were observed within V8 and V9.

Notably, structural geographic differences were mainly related with variations in the relative abundances of Proteobacteria, Fusobacteria, and Bacteroidetes rather than through the acquisition of new bacterial species. This observation may further support the existence of a bacterial core in pinnipeds (Glad et al., 2010; Nelson et al., 2013; Numberger et al., 2016).

Geographic GI microbiota patterns have been reported for sea lions as well as other pinnipeds, and these are thought to be associated with several factors, including diet, site fidelity, and colony density (Nelson et al., 2013; Delport et al., 2016; Pacheco-Sandoval et al., 2019). Although we did not directly address the contributions of these factors on GI bacterial diversity, our experimental design may be useful for explaining additional causes of local variations in GI microbiota. Indeed, we analyzed GI samples from pups 2–3 months old. At this age, sea lion pups are not able to move between colonies and are still dependent on maternal feeding. Therefore, the observed GI geographic variation in pups may reflect the local preferences of their mothers. If accepted, this finding may also suggest that GI microbiomes in pups are vertically inherited from mother to offspring at birth and as a result of nursing, as has recently been proposed for humans (Blekhman et al., 2015; Duranti et al., 2017; Meehan et al., 2018) and other mammalian species (Antwis et al., 2018), although additional research is needed to support this hypothesis. Further, our findings support the hypothesis that once an initial microbiota core is established, it may persist to adulthood (Nelson et al., 2013), even though external factors may modulate the abundance of certain bacteria that promote the most beneficial configuration for the host. For example, Smith et al. (2013) reported that milk-based diets in seal pups promote the growth of Lactobacilli, which were present in lower abundances in adults with solid, marine-based diets.

Finally, our results do not support the observation that higher microbial diversity is usually detected in high-density colonies (Nelson et al., 2013; Delport et al., 2016). Relatively high values of bacterial community diversity were observed in Partido, one of the colonies with the lowest number of individuals (Szteren, Aurioles & Gerber, 2006; Masper et al., 2019). This observation suggests that colony density per se is not a force that shapes GI bacterial communities and that individuals with close ecological ties (e.g., similar diets) are more likely to converge and exhibit similar microbiome compositions, regardless of colony density.

Microbiome functional predictions

Given that we detected distinctive HR-related bacterial consortia, we evaluated if these alternative arrangements could lead to different functional capabilities of the microbiota. Despite the different number of molecular pathways identified among HRs, no significant HR-driven changes were apparent. In addition, similar functional profiles were also observed between rookeries. The simplest explanation of such consistent results may lie in the proposed methodological approach. The results of comparative analyses have recently suggested that the sensibility of PICRUSt and other bioinformatics tools for in silico functional predictions of microbial community is largely limited in non-human samples (Sun, Jones & Fodor, 2020). Nevertheless, the high degree of genetic similarity between the detected ASVs and reference genomes (highlighted by the low NSTI values) may indicate reliable prediction accuracy, and thus additional explanations may be considered.

In this context, we observed that the majority of the predicted functional activities were enriched by seven bacterial genera: Escherichia, Fusobacterium, Otariodibacter, Bacteroides, Ornithobacterium, Porphyromonas, and Pseudomonas. Hence the GI microbiota of sea lion pups may be characterized by a high degree of functional redundancy (i.e., different bacterial species bearing similar ecological functions). As has been reported in other species, functional redundancy appears to be an intrinsic property of mammalian GI ecosystems that sustains homeostatic conditions, which ultimately guarantees host health (Lavery et al., 2012; Pérez-Cobas et al., 2013; Moya & Ferrer, 2016). This condition may be particularly evident in sea lion pups given that their GI microbiomes are primarily involved in metabolic activities. According to recent genomic annotations, the genomes of Bacteroidetes and Fusobacteria are known to present a high number of genes involved in the degradation of glycans, proteins, and complex oligosaccharides (Wang et al., 2016; Birhanu et al., 2019). Accordingly, these genera may confer ecological advantages in sea lions by increasing fat deposition, which is necessary for thermal isolation in the aquatic environment, where heat transfer is several times faster than the terrestrial environments and supplies energy (Lavery et al., 2012; Acquarone, Salgado-Flores & Sundset, 2020).

Notably, some of these genera are known to be pathogenic members of human GI communities; however, their high abundance in predicted KEGG pathways suggests mutualistic relationships with sea lion pups. We speculate that sea lions and other pinnipeds species may act as natural reservoirs for these human pathogenic bacteria (Stewart et al., 2008; Zheludkov & Tsirelson, 2010; Sato et al., 2019), which switch from being commensals to pathogen species after their introduction into “non-native” host species, although additional research is needed to support this hypothesis.