Skin microbiome of beluga whales: spatial, temporal, and health-related dynamics

Abstract

Background

Host-specific microbiomes play an important role in individual health and ecology; in marine mammals, epidermal microbiomes may be a protective barrier between the host and its aqueous environment. Understanding these epidermal-associated microbial communities, and their ecological- or health-driven variability, is the first step toward developing health indices for rapid assessment of individual or population health. In Cook Inlet, Alaska, an endangered population of beluga whales (Delphinapterus leucas) numbers fewer than 300 animals and continues to decline, despite more than a decade of conservation effort. Characterizing the epidermal microbiome of this species could provide insight into the ecology and health of this endangered population and allow the development of minimally invasive health indicators based on tissue samples.

Results

We sequenced the hypervariable IV region of bacterial and archaeal SSU rRNA genes from epidermal tissue samples collected from endangered Cook Inlet beluga whales (n = 33) and the nearest neighboring population in Bristol Bay (n = 39) between 2012 and 2018. We examined the sequences using amplicon sequence variant (ASV)-based analyses, and no ASVs were associated with all individuals, indicating a greater degree of epidermal microbiome variability among beluga whales than in previously studied cetacean species and suggesting the absence of a species-specific core microbiome. Epidermal microbiome composition differed significantly between populations and across sampling years. Comparing the microbiomes of Bristol Bay individuals of known health status revealed 11 ASVs associated with potential pathogens that differed in abundance between healthy individuals and those with skin lesions or dermatitis. Molting and non-molting individuals also differed significantly in microbial diversity and the abundance of potential pathogen-associated ASVs, indicating the importance of molting in maintaining skin health.

Conclusions

We provide novel insights into the dynamics of Alaskan beluga whale epidermal microbial communities. A core epidermal microbiome was not identified across all animals. We characterize microbial dynamics related to population, sampling year and health state including level of skin molting. The results of this study provide a basis for future work to understand the role of the skin microbiome in beluga whale health and to develop health indices for management of the endangered Cook Inlet beluga whales, and cetaceans more broadly.

Introduction

Although marine mammals are in constant contact with a diverse and abundant seawater microbial community [1], they host their own unique [2], temporally stable [3], and species-specific [4, 5] epidermal microbial communities. These microbes may play an important role in or be indicators of marine mammal health and ecology [6, 7]. As an emerging field of study, marine mammal microbiomes considers these previously overlooked interactions between hosts and their microbiomes [8], and explores how spatial and temporal variability in microbial communities may help to explain the way marine mammal hosts interact with their environment [9].

For many mammals, their associated microbial communities frequently correlate with host health [10]; epidermal microbiomes are similarly thought to be correlated with individual health and disease. In some cases, the presence or absence of skin diseases may depend on the ability of the host microbial community to fight off pathogens [11]. Skin microbiome research has largely focused on humans and humanized model systems [12,13,14], and much less is known about the skin microbiomes of wildlife. Certain microorganisms may be useful indicators of health, which could help wildlife managers better understand the health of endangered or recovering individuals or populations. However, answering high-resolution questions about microbiomes and health first requires a baseline understanding of the core microbial community associated with a species, and environmentally driven spatio-temporal variability in community composition.

Beluga whales (Delphinapterus leucas) are a highly social and gregarious species, forming groups of two to 25 individuals, and, in some populations, aggregating in larger groups (hundreds to thousands) during the summer months [15, 16]. Their complex social structure includes dynamic group affiliations, both with close kin and non-kin, forming large-scale communities comprising both sexes and all age groups [17]. Populations occupying areas of high winter ice cover undergo long-range seasonal migrations to ice-free areas; populations outside the high Arctic do not necessarily undergo similar long range movements, but their movements are often closely linked to available prey, such as salmonids returning to breed in certain river systems [18]. Although fish makes up a large portion of their diet, they are also known to eat benthic invertebrates from shallow estuaries, as well as deep-sea squid and other invertebrates [19]. Longer migrations by arctic beluga whale populations may be related to their apparent need for warmer water temperatures for their annual summer molt [20, 21], similar to other polar populations of marine mammals [22,23,24], which are thought to be important to maintaining skin health. The annual molt occurs in contrast to most non-polar cetaceans, in which epidermal turnover is a continuous process, and is thought to be partly due to colder water temperatures that slow the breakdown of intercellular lipids and allow superficial cells to remain attached for longer [20]. In beluga whales, microorganisms have been visually observed to be abundant in dead or damaged tissue that remains attached to the surface before the annual molt [20]. The epidermal microbiome of a beluga whale is likely to reflect its complex life history, annual molting cycle, and highly social behavior, as well as the rapidly changing arctic and sub-arctic environments in which it lives.

Bristol Bay and Cook Inlet are both shallow estuarine bays with strong tidal influences, each comprising multiple large river systems that are spawning grounds for several large runs of salmonids, and are home to the two most southerly beluga populations in Alaska. Both populations remain resident to these tidal estuaries throughout the year and feed on flooded tidal plains in the summer, moving out of the riverine and shallow nearshore areas and into deeper water when ice begins to form [25,26,27,28,29]. In both populations, Pacific salmon (Oncorhynchus spp.) make up the main portion of the diet during May through October, supplemented in differing proportions by smelt, cod, and invertebrates (primarily shrimp) [19]. Bristol Bay is remote and has remained relatively free of anthropogenic impact over the years, aside from a flourishing salmon fishery [30], while Cook Inlet has seen an increase in human-caused disturbances due to rapid industrial growth along the coast, which is home to approximately half of Alaska’s population, an international airport, the largest vessel port in the state, and a military base, and is downstream of Alaska’s largest agricultural region [31,32,33]. These activities have led to increased ocean noise and degraded water quality within Cook Inlet, as well as potential impacts on prey availability, among other threats [27, 30, 33, 34].

While the size of the beluga population in Bristol Bay is approximately 2000 individuals and has been steady or increasing [35,36,37], the population of beluga whales in Cook Inlet, Alaska, declined severely in the 1990s, attributed to unrestricted hunting by Alaska Natives [38, 39]. It has not begun to recover despite a regional moratorium on hunting beginning in 2006 and was listed as Endangered under the U.S. Endangered Species Act in 2008. Estimated to comprise 1300 individuals in the early 1970s, it now numbers fewer than 300 individuals and was estimated to be declining at more than 2% each year for the period 2008–2018 [40]. Cook Inlet belugas face a number of possible threats, including contamination, pathogens, noise, habitat degradation, ship-strikes, and declines in available prey. Understanding this population’s barriers to recovery requires a holistic understanding of the population’s structure, ecology, and health. These aims may be aided by an understanding of their microbiota.

In this study, we characterize and compare the epidermal microbiota of two beluga whale populations: the endangered population in Cook Inlet and the nearest neighboring, purportedly healthy population, in Bristol Bay. This is the first study of the microbial communities associated with beluga whales. Our aims are 1) to describe the core epidermal microbial community associated with Alaska beluga whales, thus providing insight into the beluga-specific microbiome, 2) to understand variability in epidermal microbial communities, and factors that may affect that variability, and 3) to examine links between beluga epidermal microbial communities and beluga health that may improve our understanding of population-level health and aid the conservation of endangered populations.

Methods

Sample collection

Samples were collected between 2012 and 2018 in Cook Inlet (n = 42) and 2008 and 2016 in Bristol Bay (n = 44). Most samples from Cook Inlet were collected during biopsy surveys in 2016–2018, with a small number of additional samples from recent strandings (n = 5). Biopsy samples were collected from the upper flank near the dorsal ridge, using a remotely deployed biopsy dart to extract a section of skin and blubber tissue approximately 8 mm in diameter and up to 20 mm in length. Samples were transferred to a liquid nitrogen dry shipper immediately following collection and stored at − 80 °C prior to processing. In Bristol Bay, tissue samples were collected from whales that were captured near shore (n = 44) following the protocol outlined in Goertz et al. [41]. Individual whales longer than 250 cm, excluding mom-calf pairs, were first followed in an 18 ft. aluminum skiff until they entered shallow water, at which point the individual was corralled using a net held between the first boat and an additional boat. During capture the animal was restrained using a tail rope, hoop nets, slings, and additional staff as needed. Each individual was held for less than 2 h, during which time it underwent a visual health examination and sex determination, measurements, biological sampling, and satellite tag attachment. Skin and blubber biopsies up to 8 mm in diameter and less than 2 cm deep were taken along the flank below the dorsal ridge to closely approximate the location of remote dart biopsies. Once collected, samples were stored on dry ice in the field and later transferred to liquid nitrogen before final storage at − 80 °C prior to processing.

Microbial sequencing

DNA was extracted from healthy skin tissue (no lesions) from animals in all states of health, so that the proximity of the sample to diseased tissue could be factored out of resulting analyses. Extractions used approximately 2 mm³ of skin tissue, taken from the outer epidermal layer whenever possible, and a DNEasy Blood and Tissue Kit (Qiagen). The kit was modified by first homogenizing the tissue, and then replacing the lysis buffer with the following buffer mix: 20 mM Tris HCl pH 8.0, 2 mM EDTA pH 8.0, 1.2% v/v Triton X100 [4]. Hypervariable IV region of the small subunit ribosomal RNA (SSU rRNA) genes of bacteria and archaea were amplified in triplicate for sequencing. Each PCR contained 1 μ l of DNA, 5 μ l of GoTaq 5X Flexi buffer (Promega), 14.75 μ l of H₂O, 2.5 μ l of a 25 nM MgC l₂ solution, 0.5 μ l dinucleoside triphosphate (dNTP; Promega), and 0.25 μ l of a 5 units/ μ l GoTaq DNA polymerase solution per sample, along with 200 nM each primer (515FY and 806RB), which included unique forward and reverse barcodes for each sample [42, 43]. Amplification and barcoding occurred in the same PCR step. PCR cycle number was optimized on a per-sample basis, ranging from 35 to 38 cycles. The amplification protocol was 2 min at 95 °C, 35 to 38 cycles of 20 s at 95 °C, 15 s at 55 °C, and 5 min at 72 °C, followed by 10 min at 72 °C. The triplicate PCR products were pooled for each sample and visualized on a 2% agarose gel with 1% 10,000X SyberSafe dye, and bands were excised from the gel and purified using a MinElute gel purification kit (Qiagen). The gel-purified product was quantified using the Qubit assay (Invitrogen), and barcoded products were pooled in equimolar concentrations for 2 × 250 bp MiSeq sequencing (Illumina) at the Georgia Genomics and Bioinformatics Core, University of Georgia. A negative control was included with each extraction batch and PCR amplification run. These negative controls were carried through to gel-visualization, at which point no product was observed. To evaluate the effectiveness of the amplifications and sequencing runs, two mock communities were amplified and included with each sequencing run: 1) an equimolar community, and 2) a staggered ribosomal RNA operon (BEI Resources, NIAID, NIH, Manassas, VA, USA, as part of the Human Microbiome Project: Genomic DNA from Microbial Mock Community B, HM-782D and HM-783D, v5.2 L). A total of 19,404,299 raw reads were generated from 102 samples in two Illumina MiSeq sequencing runs, including 4 mock communities in two sequencing runs, 8 duplicates (identical samples from the same extraction that were independently amplified with unique barcodes), and 4 samples with no associated metadata, which were removed prior to data analysis. Sequences were deposited in the NCBI’s Short Read Archive under accession numbers SAMN15503410 - SAMN15503499.

To test sequencing accuracy, the program dada2 v1.12 [44] inferred 28 sequences in the mock community, 21 of which were matches to the expected reference sequences. The 7 sequences that did not match the reference sequences were compared to the NCBI database using BLAST. The most common similarities were with uncultured bacterial sequences, although one matched E. coli and another matched Rhodobacter sphaeroides. The non-reference sequences were not found in any of the samples in the dataset, so it is unlikely that these sequences were the result of contamination in the lab.

Sequence processing

The resulting sequences were processed using the dada2 pipeline, implemented in R v3.6.1 [45]. Due to the high similarity between duplicates, one duplicate from each identical pair was chosen to represent duplicated samples (Supplemental Figure S1). Sequences were first trimmed to remove low-quality regions at the tails by cutting sequences once the Phred score fell below 10. Next, sequences were removed if they had more than two expected errors per sequence, calculated based on the Phred quality score:

$$ EE=\sum \limits_{n=1}^{seqlenth}\left({10}^{\left(-Q/10\right)}\right) $$

Filtered forward and reverse sequences were then dereplicated based on the estimated error rate for each possible nucleotide transition, unique sequences were identified using pseudo-pooling, and paired reads were merged to generate full denoised sequences. Denoised sequences were clustered into amplicon sequence variants (ASVs), and chimeras were removed. ASVs were favored over more traditional OTUs for increased precision, taxonomic resolution, and reproducibility [46]. After quality filtering with dada2, 12,147,026 reads remained. At this stage, samples with fewer than 10,000 reads (n = 8) were removed from the dataset (mean reads = 122,163.2 ± 75,922.46). Additionally, the 5 samples from stranded individuals were determined to have significantly different microbial communities compared to non-stranded individuals in Cook Inlet, based on a PERMANOVA of Bray-Curtis dissimilarity metrics [47], and were removed from the dataset prior to further analysis. Taxonomy was assigned using the Silva rRNA sequence database (version 132) [48]. Non-target sequences included 16 ASVs that could not be classified by the Silva database; these were determined by an NCBI search to be cetacean, fungal, algal, or parasite and manually removed from the dataset. Removing eukaryotic sequences caused one additional sample to fall below the minimum threshold of 10,000 reads (sample DL8–8), and it was also removed from the final dataset. Quality pruning of the dataset resulted in a final dataset comprising 72 samples (33 in Cook Inlet, 39 in Bristol Bay).

Data analysis

Geographic and temporal variability in the core microbial community, as well as variability driven by sex and age, were assessed using a Principal Coordinates Analysis (PCoA) of Bray-Curtis dissimilarity, implemented in the package phyloseq [49], and Permutational Multivariate Analysis of Variance (PERMANOVA) implemented using the package vegan [50]. Observed richness, as well as Shannon and Simpson alpha diversity metrics [51, 52], were also calculated using phyloseq; similarly, beta diversity [53] within each population was estimated using pairwise Bray-Curtis dissimilarity among individuals within that population.

Correlation between epidermal microbial communities and individual health status was tested on a subset of the samples in the study, comprising individuals live-captured for sample collection in Bristol Bay. At the time of capture, qualitative health assessments were performed by a veterinarian from the Alaska SeaLife Center (CECG). Based on these health assessments, animals were designated as either healthy (n = 14) or exhibiting evidence of skin disease (n = 15). Because many marine pathogens remain undiscovered, we looked for possible pathogenic microbes by identifying ASVs in the same genera as known pathogens [54]. We tested for differences in abundance of these potential pathogens between healthy and diseased individuals using an analysis of differential representation of ASVs, implemented in DESeq2 [55].

Health assessments for Bristol Bay animals also identified individual molt status; using this information we stratified all healthy individuals into two categories based on whether they were molting (n = 4) or non-molting (n = 10). Epidermal microbial communities of molting and non-molting individuals were compared using alpha diversity indices, estimated using the package phyloseq. We further looked for specific microbes from the entire dataset, as well as specific microbes from the potential pathogen dataset, with differential abundances between the two groups using DESeq2.

Results

Sequence clustering of the 7,323,872 quality, target SSU rRNA gene sequencing reads (average per sample = 95,115, range of 16,296 - 238,307) from Bristol Bay and Cook Inlet skin microorganisms generated 13,532 ASVs among the 72 samples. The median percent per sample of non-target, eukaryotic reads removed from analysis was 0.028%; 9 individuals had > 1% non-target, eukaryotic reads. A total of 75.5% of target ASVs were classified to family, and 50.55% to genus.

Archaea were represented by 158 ASVs, primarily from the classes Nitrososphaeria and Thermoplasmata. Archaea were found in 52 samples, in both Cook Inlet and Bristol Bay, and comprised an average of 1.4% of the skin microbial communities. One outlier, sampled in Cook Inlet in 2018 (Sample DL18–6, Supplemental Table S1), had a microbial community with 42% archaeal composition. Archaea comprised < 5% of all other samples. The remaining ASVs classified as Bacteria.

Beluga core epidermal microbial community

Beluga whale microbial communities were dominated by the phyla Bacteroidetes and Proteobacteria, which made up an average of 11 and 51% of individual community composition, respectively. No individual ASVs or classified genera were present in all samples included in the study, but two families were present in all samples: Moraxellaceae (mean 18%) and Burkholderiaceae (mean 5%). Individual traits, including sex and age-class, had a significant but small effect on microbial community composition when all samples were considered (age-class PERMANOVA on Bray-Curtis dissimilarity; R² = 0.053, p-value = 0.017, sex PERMANOVA R² = 0.02, p-value = 0.035).

Geographic and temporal variability in microbial community

No individual ASVs were present in all individuals in either the Cook Inlet or Bristol Bay population. In Cook Inlet, no classified genera were present in all individuals. However, in Bristol Bay, three genera were present in all individuals: Acinetobacter, Flavobacterium, and Pseudomonas. Those genera each comprise 1309, 798, and 833 ASVs, respectively, within the Bristol Bay population.

Epidermal microbiome alpha diversity differed significantly between the two regions in three major indices: observed richness, Shannon, and Simpson indices (PERMANOVA R² = 0.28, p = 0.001, Fig. 1), with higher diversity in the Bristol Bay population across all three indices. Microbial community composition also differed significantly between regions, based on a Principal Coordinate Analysis (PCoA) of Bray-Curtis dissimilarity (Fig. 2a) and PERMANOVA testing (R² = 0.055, p-value = 0.001). Beta diversity (dispersion) was high within each region, indicating a large degree of dissimilarity in individual microbiomes within each of the two populations, and was significantly higher in Cook Inlet than in Bristol Bay (one-sided Wilcoxon p-value < 0.001), with greater variability in beta diversity in Bristol Bay than in Cook Inlet (Fig. 2b). Bacteroidetes and Firmicutes were common on beluga whales from Bristol Bay, while individual communities from Cook Inlet were dominated by Proteobacteria (Fig. 2c).

Fig. 1

Observed richness (number of ASVs), Shannon, and Simpson indices of alpha diversity among epidermal microbiomes in Bristol Bay (BB) and Cook Inlet (CI) beluga whales

Fig. 2

Comparison of epidermal microbial communities in Bristol Bay and Cook Inlet beluga whales. a PCoA based on Bray-Curtis dissimilarity comparison of the epidermal microbiota. Multivariate normal and t-distributional data ellipses are shown in dotted and solid lines, respectively. b Beta diversity within each regional population, estimated by calculating Bray-Curtis dissimilarity between individual microbial communities within populations. c Proportional abundances of top 20 ASVs, organized by phylum, found on individual whales in Bristol Bay and Cook Inlet

Inter-annual variability was apparent among years within both populations (Fig. 3), with high inter-individual variability within years. Bristol Bay skin microbial communities were dominated by Proteobacteria, including genera such as Acinetobacter in 2008 and SAR11 clade Ia in 2016, with a shift toward Bacteroidetes (Flavobacterium sp.) and Firmicutes (Staphylococcus sp. and Clostridium sp.) in the interim years. We did not have a sufficiently long time series to look for a similar pattern in the Cook Inlet population, but in that population differences were apparent between 2017 and 2018 at the sub-phylum level within the Proteobacteria. When region and year are considered in the same PERMANOVA, year has a strong effect on variability in epidermal microbial communities (R² = 0.19, p = 0.001), and the geographic effect of region is diminished, i.e. R² falls from 0.055 to 0.039 (p = 0.001). Within Cook Inlet, epidermal microbiomes differed significantly among the three locations where samples were collected (PERMANOVA R² = 0.097, p-value = 0.001), as well as among age classes (PERMANOVA R² = 0.073, p-value = 0.036), although these effects may be conflated with annual effects (Fig. 3b). Within Bristol Bay, the magnitude of effect of location was greater than in Cook Inlet among samples in this study, but it was non-significant (PERMANOVA R² = 0.15, p-value = 0.058). The effect of age class in Bristol Bay was similar to that in Cook Inlet (PERMANOVA R² = 0.069, p-value = 0.039). In this population as well, the effects of sampling location, age class, and sampling year may be conflated (Fig. 3c).

Fig. 3

Annual variability within each regional population. a Variation in proportional abundance of top 20 ASVs associated with individual whales in Bristol Bay (top row) and Cook Inlet (bottom row), organized by genus. Each bar represents an individual sample, grouped by year along the x-axis. The bottom two plots are PCoA analyses based on Bray-Curtis dissimilarity comparison of epidermal microbiota, with color by year, shape by sampling location, and label by sex, in b Bristol Bay and c Cook Inlet

Comparing microbial communities to the database of human and animal pathogens developed previously [54] at the genus level revealed 2,313 ASVs associated with pathogenic genera in both populations. The most abundant potentially pathogenic genera across all samples were Moraxella and Acinetobacter. The composition of pathogenic microbes varied with year (PERMANOVA R² = 0.29, p-value = 0.001), but not between the two regions.

Epidermal microbial variability and individual health status

Epidermal health status was assessed for 29 individuals live captured and released in Bristol Bay. Overall, healthy individuals and those with skin disease did not exhibit significant differences in the proportion of potentially pathogenic ASVs present (t-test p-value = 0.14) or composition of potential pathogens (PERMANOVA R² = 0.027, p-value = 0.922). However, testing for differential abundances of each potential pathogen between diseased and healthy individuals revealed 11 genera with significantly different abundances (using a p-value threshold of 0.05) between the two groups (Fig. 4). Of those, two ASVs were more abundant in diseased animals: Klebsiella sp. and Psychrobacter sp.

Fig. 4

Log2 fold change (i.e. log-of-the-fold-change) in abundance of potentially pathogenic genera that differed significantly between healthy and diseased individuals. The magnitude and direction of change is displayed as the change in abundance of microbes on diseased individuals when compared to healthy individuals (i.e., genera enriched in diseased animals are positive values and those enriched in healthy individuals are negative values)

Molting and non-molting animals in Bristol Bay exhibited differences in overall microbial community diversity, with a higher diversity of potentially pathogenic ASVs in individuals that were not molting (PERMANOVA R² = 0.34, p = 0.019, Fig. 5). It is important to note the small sample size of molting individuals (n= 4), which often decreases precision in distributional means. An analysis of differential abundance of potentially pathogenic microbes between molting and non-molting individuals isolated 12 genera with significantly different abundances (p < 0.05) between the two groups (Fig. 6), five of which were more abundant in non-molting individuals.

Fig. 5

Observed richness (number of ASVs), Shannon, and Simpson indices of alpha diversity among epidermal microbial communities in actively molting and non-molting animals live-captured in Bristol Bay

Fig. 6

Log-fold-change in abundance of potentially pathogenic genera that differed significantly between molting and non-molting individuals. The magnitude and direction of change is displayed as the change in abundance of microbes on non-molting individuals when compared to molting individuals

Discussion

No core epidermal microbial community in belugas

Unlike previous studies of cetacean microbiomes [2, 4], the results of this study do not present strong evidence for a species-specific core microbial community in Alaskan beluga whales. No ASVs were found on all whales included in the study. While this may be attributed partially to the long-term isolation of these two populations [56], it was also true that no ASVs were found in all whales within either of the two populations. The variability in epidermal microbiome could be affected by a number of factors (discussed further below), including social group membership, annual or seasonal variability in environmental factors or individual behavior, or preferred habitat within each region.

Individual and population health could also be important drivers of increased microbial diversity and the lack of a core microbiome in this study, as per the Anna Karenina Principle (AKP) in microbiomes, which posits that ‘all healthy microbiota are all alike; each unhealthy microbiome is dysbiotic in its own way’, i.e. that healthy individuals will have similar microbiota while diseased individuals will each deviate from the “healthy microbiome” in a unique way, such that the very presence of disease can increase microbiome beta diversity stochastically [57]. In human microbiomes this principle has been shown to hold true for approximately 50% of diseases, while the remaining diseases exhibited either anti-AKP effects (25%) or effects unrelated to AKP (25%) [58]. In this study, beta diversity was high in both populations, and higher overall in Cook Inlet than in Bristol Bay. However, at a population level it is difficult to interpret beta diversity as an indicator of disease, as individual diseases have a nearly equal likelihood of correlation with increased beta diversity vs. decreased or equal beta diversity in individual host microbiomes; thus the composition of diseases in a population may have more effect on beta diversity than the abundance of diseases in that population.

Seawater control samples were not available as part of this study; therefore, it is important to consider that some of the documented microbiome variability may be contributed by seawater microbiota. Previous studies indicate that the contribution from seawater microbes is likely limited, as marine mammals harbor microbial communities that are distinct from their surrounding environments [2, 4, 9]. However, the collection of seawater control samples as part of future studies will allow for the determination of whether and how much variability in the epidermal microbiome of belugas is affected by seawater-associated microbes.

Only 50% of ASVs were classified to genera, highlighting our paucity of knowledge when it comes to marine mammal-hosted microbial communities, and the need for continued study of the relationships between these animals and their microbiota. Proteobacteria, Bacteroides, and Firmicutes were the most common phyla found in this study; of these, the first two are common in both cetacean and seawater microbiomes, and Firmicutes has also been documented in marine mammal oral and rectal microbiomes [2]. The lowest taxonomic level found to be common across all samples was that of family: Moraxellaceae and Burkholderiaceae were present on all samples in both populations in this study. Moraxellaceae is a family of gammaproteobacteria, including known psychrotrophic species that would be well-suited to the polar habitat of beluga whales [59, 60]. The family includes both pathogenic and commensal species, but the majority of ASVs from the family Moraxellaceae that were documented in this study remain unidentified and unstudied, therefore it is unknown whether they are pathogenic. The second family occurring across all samples, Burkholderiaceae, is a broadly diverse family represented by over 50 genera in this study. This family includes multiple examples of genera that tolerate or thrive in extreme cold environments, including Janthinobacterium [61, 62] and the psychrophile Polaromonas [62, 63]. The family also includes a number of genera more commonly associated with wastewater or contaminants, such as Extensimonas [64] and Sphaerotilus [65,66,67], which are found in industrial wastewater and activated sludge, and were both documented on a small number of individuals from each population. In Bristol Bay, these animals were sampled near Dillingham; in Cook Inlet they were sampled in Eagle River in Knick Arm, and in Chickaloon Bay. While it is possible that these bacteria are part of the core epidermal microbial community, they are more likely indicative of the whales’ interactions with contaminated nearshore waters in both regions. The presence of these genera in the epidermal microbiome should not be overlooked; future studies should consider the extent and concurrence of wastewater microbiota on beluga whales and their preferred habitats, as well as the potential effects to beluga health of these microbiota or associated contaminants.

While no genera were common to all sampled individuals in Cook Inlet, three genera were found on all individuals in Bristol Bay: Acinetobacter, Flavobacterium, and Pseudomonas. Pseudomonas has previously been documented among delphinids and phocids, and may be a cosmopolitan-type associate of marine mammals [5]. The genus contains species with a wide range of metabolic diversity and niche space, and includes several known pathogens of humans, non-human animals, and plants [68]. Further research is required to determine the role of this genus in marine mammal microbiota, and whether the organisms in this genera have a similar function across all host relationships. Acinetobacter has also been documented in humpback, grey, and pilot whale skins [69, 70], is common in seawater, and is able to live in a wide range of temperatures and environmental conditions [71]. Coral-associated species of Acinetobacter have been shown to produce chemicals that inhibit biofouling [72]. However, this genus also harbors potentially pathogenic species [54]; the abundance of this genera warrants further research into to functional roles of beluga-associated species within the Acinetobacter, and whether those species are harmful or beneficial. Finally, Flavobacterium has been previously identified in the core epidermal microbiome of humpback whales [69]; the genus also includes several pathogens of freshwater and saltwater fishes [73, 74].

It is worth noting that Archaea, which have never before been documented in marine mammal microbiomes despite some use of primers that could amplify members of this domain, were present in nearly 70% of the samples in this study, although generally in low proportions. The paucity of information on Archaea makes it difficult to hypothesize the cause for the outlier with 42% archaeal composition, which was an adult male sampled in Cook Inlet in 2018; however, there is no indication that it was caused by contamination. The domain Archaea is probably most well-known for comprising extremophilic species, but have been found to be part of the microbiota of many organisms [75,76,77], and are especially numerous in marine environments [78, 79]. They are an important part of the epidermal microbiota of humans [80], in similar proportions to that found here on beluga whales [81]. In humans, they are thought to be important to nitrogen turnover and skin health [81]; they may have a similar role in the epidermal microbiota of beluga whales and other marine cetaceans. There are no known pathogenic or parasitic Archaea; they are more often mutualist or commensals.

Geographic and temporal variability in epidermal microbiome

Spatial and temporal variability are common in host-associated microbiomes [69, 82,83,84]. The majority of microbiome research is still focused on humans and humanized model systems, and that research shows geographic and temporal variability in microbiomes even among healthy adults [85]. Characterizing this type of variability in a species, which is likely at least partially driven by environmental factors, will improve our understanding of a species’ interactions with its ecosystem. Further, it allows scientists to isolate and understand microbiome variability that may be related to changes in individual or population-level health, providing key baseline data for developing monitoring tools critical for evaluating individual and population-level responses to environmental change.

Beluga whale microbiotas differed both spatially and temporally, indicating a strong environmental effect on microbial community composition that should be further investigated. However, these results should be interpreted with caution due to minimal annual overlap in sampling between the two geographic populations. Future research aiming to collect samples from both regions in the same year and season would greatly enhance our ability to detect temporal and regional effects on microbiome variability.

Alpha diversity was higher in Bristol Bay than in Cook Inlet, a phenomenon that could be driven in part by endogenous population factors, such as the greater long-term abundance and density of beluga whales in Bristol Bay. When individual beluga whales are considered as habitat space for microbial populations, the existence of less habitat space in Cook Inlet may lead to increased local extinction rates, thus decreasing microbial diversity. Similarly, the higher density of beluga whales in Bristol Bay may lead to increased opportunities for interactions among microbial communities, thus increasing microbial transfer and decreasing variability among individuals. This is reflected in the beta diversity metric, which was lower in Bristol Bay than in Cook Inlet. However, beta diversity was generally high within both populations, indicating high variability among individuals - some of this is likely accounted for by other patterns detected in this study, such as annual variability within populations, molt stage, and individual health status.

Bray-Curtis dissimilarity metrics indicate strong divergence between the two populations, which could be driven by a number of factors. Genetic evidence suggests the Cook Inlet population is demographically isolated from Bristol Bay and other populations in the Bering Sea [56]. Migration between the two populations is minimal over recent (ecological) timescales, according to a Bayesian analysis of eight nuclear loci that estimate a recent migration rate from Bristol Bay to Cook Inlet of 0.004 (95% confidence interval: 0–0.03) [56]. This value approximates one individual (95% confidence interval: 0–8 individuals) over the last two generations, or roughly 58 years, based on a population of 300 individuals and a generation time of 29 years [86], although this estimate is rough due to uncertainty in all three estimates (migration rate, generation time, and population abundance), which compound in this calculation. Similarly, it is estimated from population genetic data that in the last ~ 60 years the Bristol Bay population may have received six immigrants from Cook Inlet (95% confidence interval 0–22 individuals) based on a population of 2000 individuals [35]. With connectivity this low between the two populations, their microbial communities are almost entirely isolated, which facilitates both random and adaptive divergence between them. The differences between regions could also be explained by environmental differences between habitats, including natural differences such as oceanographic or benthic properties of the two bodies of water, as well as anthropogenic differences such as nutrient and contaminant loading.

Observed annual variability within populations could also be attributed to a number of potential drivers, including endogenous factors, annual changes in environmental conditions, or anthropogenic disturbance. In this study, we documented an effect of age on microbial community composition both between the two regions and within the Cook Inlet population. Sex was also found to have a small-but-significant effect, similar to age, both of which may reflect behavioral or endogenous differences among groups (e.g. group associations, hormone levels and/or inherent biochemical differences, diet preferences). Males and females, for example, are thought to prefer different habitats [87]. If that is the case, it is possible that sex-based habitat preferences may be confounded with that of socially driven habitat preferences. Previous studies have found that water temperature can be an important driver in epidermal microbiome composition [69], and may have an effect here as well. To determine the relative contribution of endogenous and environmental factors to microbial community variability, it will be important to conduct studies that vary one factor while sampling randomly across the others, as it is not feasible in wild populations to hold the other factors constant.

Intra-population variability may also be socially driven, either by affecting the likelihood of sampling an individual in a given year, or by affecting the individual’s preferred habitat and, by proxy, its epidermal microbiome. Within Cook Inlet we observed strong differences in epidermal microbiome associated with sampling location. Individuals migrate around the inlet throughout the year, so it is unknown whether this correlation represents rapid turnover of epidermal microbiome driven by local environmental conditions, or sub-population structure driven by differences in ecological behaviors. If the former, it would be the first documentation of rapid turnover in cetacean epidermal microbiome in response to local environmental conditions, and would suggest strong variability in water quality and habitat characteristics throughout Cook Inlet. The latter possibility is more common. Ecological behaviors such as migration routes and habitat use have been shown to vary by social group membership in this species [88, 89], which could affect both annual and spatial variability within Cook Inlet. Further, this species is thought to have a highly complex, stable social structure [17, 90, 91], which may affect individual microbial communities in a number of ways, both directly (e.g. influencing patterns of transmission among individuals) and indirectly (e.g. driving habitat use and therefore exposure to specific microbes). The large range in beta diversity metrics (0.38–1) supports the hypothesis of socially driven transmission of epidermal microbiomes. It is possible that those individuals with more similar microbiomes are closely related or associate more often. Future work exploring correlations between individual microbiota and kinship, or measures of social affiliation, may indicate an important link between social structure and microbial communities, and ultimately help us understand the pathways of disease transmission in this species.

Epidermal microbial variability and individual health status

One of the goals in understanding variability in the epidermal microbial communities of cetaceans is to develop a minimally-invasive, rapid method to determine the health status of individuals [9]. Among humans, where the field of host microbiome research has been concentrated for the last decade, emerging results indicate the potential power of this research as a health indicator in more elusive species. Microbiome indicators of health outperform genome-wide association studies by nearly 20% [92], and certain microbial taxonomic groups have been linked with an increased probability of mortality [93]. In some cases, specific taxa have been linked with common diseases, such as autism [94], allergies [95], and obesity [96, 97].

Examining the differential abundance of potentially pathogenic species revealed 11 ASVs that may be key indicators of individual health status. The two ASVs that were more abundant in diseased animals were classified as Psychrobacter and Klebsiella. The first, Psychrobacter sp., was nearly 10 times more abundant in diseased animals than in their healthy counterparts. Psychrobacter are gram-negative psychrophiles often found in cold marine environments [98,99,100], whose endotoxins have been shown to cause an inflammatory response in macrophages, and which are known to cause human infections including endocarditis and peritonitis [101]. Yet Psychrobacter are abundantly associated with the skin microbiome of healthy humpback whales [9, 69], thus indicating that much is still unknown about how these microbes interact with their hosts. The second, Klebsiella sp., was > 20 times more abundant in diseased animals. Klebsiella display a wide range of specialized adaptations that allow them to be present in nearly all habitats and niches [102]. They are opportunistic pathogens in humans, with species known to cause a variety of internal diseases and soft-tissue infections [103]. Recent research outside of humans indicates Klebsiella may comprise species that infect non-human animals, as well [104]. Identifying and studying the functional role of the particular microbes living on beluga skin could be an important step toward understanding beluga health, but further work is needed to understand their functional roles and test their reliability as indicators of disease in an individual or population.

Nine additional potentially pathogenic ASVs in seven genera were less abundant in diseased animals, which was not expected and may indicate that, despite their close relationship to known pathogens, the ASVs represented by these genera are commensal and may be important to maintaining healthy skin. Additional research should be directed toward isolating the ASVs identified here in order to study their taxonomy and function, and to determine their role in beluga health. It is important to note that the sample size used for this analysis was moderate (n= 29), and samples were collected in a number of years. As shown above, variability among years was high, therefore additional variability due to sample size may either diminish or magnify differences in abundance of pathogenic microbes. In order to control for outside effects on microbial variability when testing the role of these genera in beluga health, it will be necessary to collect a large, random sample of epidermal tissue from healthy and sick individuals from the same population and year.

Despite a relatively small sample size of molting vs. non-molting individuals we found differences in microbial community diversity (Fig. 5), as well as 12 potentially pathogenic microbes that differed in abundance between molting and non-molting individuals, five of which were more prevalent in non-molting individuals than those that were molting. This finding suggests a possible role of molting in epidermal maintenance; the general reduction in bacterial diversity likely helps to remove harmful bacteria to maintain skin health in beluga whale populations, although the same limitations of the above analysis are true here, and a directed study of molting vs. non-molting individuals would be necessary to confirm the results reported here. The annual molt of beluga whales has long been thought to be important to maintaining skin health by removing damaged or dead epidermal tissue that may be colonized by abundant microbiota [20]. As this is the first study of epidermal microbiomes in a polar cetacean species that undergoes an annual molt, it is worth noting that the lack of a core microbiome may be related to the large differences found in the microbiota of molting vs. non-molting animals, and that future studies may compare microbiomes among groups stratified by molt stage in order to better understand the effect of molting on beluga epidermal microbiota.

Conclusion

This study represents an important first step toward understanding epidermal microbial communities in beluga whales, as well as their role in beluga health. Describing natural variability within the healthy Bristol Bay population, and comparing that with the endangered Cook Inlet population, may serve as a baseline for comparisons among healthy and affected populations. Additionally, isolating healthy and diseased individuals and describing differences in their epidermal microbiomes allows us to begin to develop indices for rapid, less-invasive health assessments of these animals, using tissue biopsies or swabbing rather than live captures. We found no support for the presence of a core microbial community in beluga whales, unlike other species studied to date, and found that within each population, microbial communities shifted annually; however, with only two populations compared it is possible that there is a common biota among most populations and one or both of these are anomalous. Intra-population variability could not be fully explained by this study; it is possible that it may be further explained by individual differences in ecological behaviors or social group membership. Within Bristol Bay, we identified individual microbes that differ in abundance between healthy and diseased animals, and which may be used to develop health metrics that can be applied to beluga whales in the endangered Cook Inlet population, thus potentially providing managers some measure of insight into the question of why this population has not recovered in the nearly 15 years since the 2006 hunting moratorium.