Respiratory microbiota of humpback whales may be reduced in diversity and richness the longer they fast

Humpback whales endure several months of fasting while undertaking one of the longest annual migrations of any mammal, which depletes the whales’ energy stores and likely compromises their physiological state. Airway microbiota are linked to respiratory health in mammals. To illuminate the dynamics of airway microbiota in a physiologically challenged mammal, we investigated the bacterial communities in the blow of East Australian humpback whales at two stages of their migration: at the beginning (n = 20) and several months into their migration (n = 20), using barcoded tag sequencing of the bacterial 16S rRNA gene. We show that early in the fasting the whale blow samples had a higher diversity and richness combined with a larger number of core taxa and a different bacterial composition than later in the fasting. This study provides some evidence that the rich blow microbiota at the beginning of their fasting might reflect the whales’ uncompromised physiology and that changes in the microbiota occur during the whales’ migration.

The airways of cetaceans [1][2][3][4][5][6][7] and other mammals (humans 8,9 , horses 10, dogs 11, cats 12, mice 13 harbour a large diversity of bacteria. Studies focusing on airway microbiota of humans [14][15][16][17] and equines 10 have identified a close relationship between the composition of microbial communities and the airways' physiological state. A compromised respiratory system typically correlates with an altered microbiota 17, as it changes the rate of immigration and elimination of bacteria as well as the growth conditions within the airways. The cause that compromises the airways determines the type of change taking place in the microbial communities. In humans, certain conditions like advanced chronic obstructive pulmonary disease (COPD) 18, chronic rhinosinusitis (CRS) and pneumonia result in a decrease of community richness and in the case of pneumonia also in the dominance of few or even a single opportunistic pathogen, like Staphylococcus aureus or Pseudomonas aeruginosa 19 . In contrast, other conditions like asthma typically cause an increase in bacterial diversity 20 . Dickson et al. 17 summarised these points by adapting a famous Tolstoy quote: ' All healthy lungs are alike,every unhealthy lung is unhealthy in its own way' . Dickson et al. 21 and Dickson et al. 17 have suggested a model, characterizing the interplay between respiratory microbiota and their host not only as dynamic and continuous, but also as bidirectional. In other words, any changes taking place within the respiratory system are capable of activating the so called dysbiosis-inflammation cycle leading to shifts of both immune response and microbiota.
Several studies investigating the airway microbiota of whales 2,3,7 and dolphins 1,[4][5][6]22 analysed their exhaled breath condensate or 'blow' . Unlike the studies on humans, the cetacean-focused studies exclusively characterized the blow microbiota of the studied specimens and did not investigate a potential correlation of blow microbiota and overall health. While this would be possible for captive cetaceans, it is often difficult to gain reliable health parameters on whales and dolphins in their natural habitat 2,23 . To gain insight into the dynamics of cetacean Scientific RepoRtS | (2020) 10:12645 | https://doi.org/10.1038/s41598-020-69602-x www.nature.com/scientificreports/ airway microbiota in response to physiological challenges, we focused here on East Australian humpback whales (HW), Megaptera novaeangliae, (breeding population E1 as defined by the International Whaling Commission 24 during their annual migration. With a journey that can exceed 8000 km, HW undertake the longest seasonal migration of all mammals [25][26][27][28] . The East Australian HW leave their feeding grounds in Antarctica at the end of the austral summer to travel north along the East Australian Coast, until they reach their low latitude breeding and mating areas on the Great Barrier Reef 29 (Fig. 1). In late winter, they return south to feed in the Southern Ocean. Chittleborough et al. 29 have reported that the whales' food intake is marginal during the breeding season, thus resulting in a period of fasting of at least 4 months per year. In addition to the energy requirements for the migration, female HW also spend massive resources on gestation and lactation 30, whereas males invest large amounts of energy in competing for females 31 .
We hypothesise that an extended period of fasting compromises the whales' physiological state, which results in a change of the microbial community composition of the airways. We therefore compared the blow microbiota of HW at the beginning of their northern migration along the East Australian Coast (HumpbackNM), where they were at the beginning of their fasting, to those after approximately 3-4 months of fasting during the southern leg of their migration (HumpbackSM).

Results
Dataset overview. We collected 20 blow samples from HW, seven seawater (SeawaterSM, environmental controls) and six air (technical controls) samples in Hervey Bay, Southern Queensland, Australia, during August 2017. The HW enter Hervey Bay on the southern leg of their migration a few weeks after leaving their breeding grounds at the Great Barrier Reef and 3-4 months after leaving 33,34 their feeding grounds in Antarctica. We compared their blow microbiota based on 16S rRNA gene sequencing with the blow microbiota from 20 different individuals in the same population of HW along with six air samples collected by Pirotta et al. 3 off the coast of Sydney, New South Wales, Australia, during May and June 2017, about 1 month after the whales left their feeding grounds in Antarctica, and 26 seawater samples (SeawaterNM) obtained from the same area as Pirotta's HumpbackNM samples 35 .
We detected a total of 11,573,157 raw 16S rRNA gene sequences, which were clustered into 8838 zero-distance operational taxonomic units (zOTUs). We removed 336 zOTUs as they had an overall relative abundance of less than 0.0001% and deleted 465 zOTUs from the dataset, which were present and highly abundant in most air samples (technical controls). The resulting dataset of whale and seawater samples contained 8,037 zOTUs with a mean of 78,827 reads per sample (sd = 116,046). The rarefaction curves (Supplementary Fig S1) Table 2 and Supplementary Table S1. We also calculated and compared the same alpha diversity values with unrarefied counts (Supplementary Table S2) and received very similar results (Supplementary Table S1).
The beta diversity analysis based on the Bray-Curtis dissimilarity coefficient of unrarefied relative abundances and ordinated in a non-metric multidimensional scaling (NMDS) plot ( Fig. 2) showed group-specific and relatively tight clustering of the HumpbackSM, SeawaterSM and SeawaterNM replicates, whereas the replicates of HumpbackNM samples did not form a distinct cluster. An overlap between whale blow and seawater samples, especially in the case of HumpbackNM and SeawaterNM was observed. The PCoA plot of the weighted and unweighted generalized UNIFRAC distances (Supplementary Figures S2 and S3) displayed a tighter clustering of the HumpbackNM samples, whereas the HumpbackSM samples were more spread out. To determine if the composition of bacterial communities was significantly different between the four groups (SeawaterNM, Hump-backNM, SeawaterSM, HumpbackSM), we fitted negative binomial models to each zOTU. We detected a significant difference between the bacterial communities of HumpbackNM and HumpbackSM (sum-of-LR = 16,626, p = 0.001), of HumpbackSM and SeawaterSM (sum-of-LR = 32,792, p = 0.001), of HumpbackNM and Seawa-terNM (sum-of-LR = 86,687, p = 0.001) and of SeawaterNM and SeawaterSM (sum-of-LR = 35,528, p = 0.001). The negative binomial models also identified those zOTUs that accounted for the majority of differences in the microbiota across the samples of HumpbackSM, HumpbackNM, SeawaterSM and SeawaterNM. The heatmap in Fig. 3 shows the prevalence of the most abundant 50 of those zOTUs. When applying the negative binomial models to the zOTUs of the blow samples of HumpbackSM and HumpbackNM only, we identified 311 zOTUs as significantly different between the two whale groups (Fig. 4). These zOTUs accounted for 15% (HumpbackSM) and 33% (HumpbackNM) of total relative sequence abundance and mostly belonged to unclassified genera. The first and second most abundant genera that could be classified were Corynebacterium (HumpbackSM: 1%; HumpbackNM: 2%) and Helcococcus (HumpbackSM: 1%; HumpbackNM: 2%).
Whales at the beginning of their fasting had a diverse and abundant core. Another useful tool to determine a change in microbial community composition across hosts or time is the identification of core taxa [42][43][44][45][46] . We used an 80% threshold (taxa present in 80% of all sampled individuals) to determine the core microbiota within the blow of HumpbackNM and HumpbackSM samples at two taxonomic levels (zOTU and genus) ( Tables 3, 4). To apply a conservative estimate of core we excluded those zOTUs that were commonly found in seawater, as we were unable to tell if these were contaminants of the blow samples or true residents of the airways. The HumpbackNM samples contained 55 core zOTUs that accounted for a total relative read abundance of 14%. 13 zOTUs were disregarded as core due to their common association with seawater. In contrast, the HumpbackSM samples did not share a single core zOTU that was not also present in seawater. The large majority of core zOTUs (42 out of 55) of HumpbackNM were most similar to 16S rRNA gene sequences detected in the mouth of bottlenose dolphins (e.g. GenBank accession: KC258936.1) 1 ( Table 1). Out of these zOTUs, 27 belonged to the class Gammaproteobacteria and more precise classification was not possible for most of these sequences. Another 13 core zOTUs of HumpbackNM were similar to sequences found on skin and nares Table 1. Number of reads and alpha diversity parameters (samples size, richness, Shannon index, Chao1 and ACE) of rarefied data (rarefied number of reads = 3432) of samples of 20 humpback whale (HumpbackNM) and 26 seawater (SeawaterNM) samples at the beginning of their fasting and of 20 humpback whale (HumpbackSM) and 7 seawater (SeawaterSM) samples at a later stage. Numbers are after deleting putative technical contaminant zOTUs from the dataset.

Species
Sample size (n) www.nature.com/scientificreports/ (KF104811.1) and in the mouth of humans and other mammals (JN713454.1). The 42 dolphin-related zOTUs accounted for a total relative abundance of 12%. The HumpbackNM samples had 20 core genera with a total relative abundance of 18%, whereas HumpbackSM had only one belonging to the genus Geobacillus, accounting for 5% (Table 4). The most abundant core genera of HumpbackNM with relative read abundances of more than 1.5% each were Arcobacter, Corynebacterium, Enhydrobacter, Helcococcus and Tenacibaculum. The two groups of HW did not share any core genera. The findings that HumpbackNM had a large number of core zOTUs and genera and therefore overall more diverse blow microbiota than HumpbackSM was also observed at class level (Supplementary Fig S4).
In addition, we determined those zOTUs and genera as 'overall' core that were shared across at least 80% of all 40 whale blow samples (HumpbackNM and HumpbackSM). The only zOTU matching these criteria was seawater-associated and therefore disregarded. Hence, the 40 whale blow samples did not share any overall core zOTUs. The same applied to the genera. Whereas Pseudomonas and Corynebacterium were shared by 29 (73%) and 28 (70%) whale blow samples, respectively, none of the genera were harboured by 80%.
Potential marine mammal-specific pathogens in HW blow. In the blow of the HW at the late stage of their migration, we found eight out of a total of 146 genera that were previously identified as pathogens in marine mammals. These eight genera accounted for a total relative abundance of 7%. The HW early in their migration harboured ten potential marine mammal pathogens that included those found in their conspecifics later in their migration. These ten genera out of a total of 294 accounted for a relative abundance of 6%. The potentially pathogenic genera shared by both whale groups included Corynebacterium, Pseudomonas, Staphylococcus, Fusobacterium, Mycoplasma, Streptococcus, Acinetobacter, Stenotrophomonas [47][48][49] . Table 5 lists the medical conditions caused by these genera and the affected species.

Discussion
We showed that the microbiota of whale blow samples were significantly different to those of their surrounding seawater (Fig. 2, Table 2). While some similarities between communities exist, these findings confirm previous reports by Bik et al. 1, Apprill et al. 2 and Raverty et al. 4, who also found the bacterial community composition of cetacean blow to be different to seawater. Moreover, we conclude that the bacterial communities of blow samples of HW at a later stage of their fasting were significantly altered compared with those at the beginning of their migration (Figs. 2, 3, 4). Different techniques in sample collection, DNA extraction, PCR cycling as well as two  www.nature.com/scientificreports/ separate sequencing runs of the two groups of samples might have contributed to the differences seen in the bacterial composition of the blow. We applied appropriate statistical techniques (e.g. subsampling, PIT-trap resampling) to minimise the impact of these factors on the alpha-and beta-diversity measurement determined here. In addition, the bacterial sequences of the seawater samples taken from the two sampling sites, but processed in the different manners outlined above, were similar in their alpha diversity indicating that sample processing and sequencing runs might overall have a negligible impact for our study. We therefore provide an indication that the HW at the beginning of their fasting exhibited a significantly higher richness and diversity in their blow microbiota (Table 1) than those from a later stage of fasting. We postulate that the relatively rich and diverse whale blow microbiota at the beginning of their fasting reflected an uncompromised, and therefore stable physiological state of the airways. The whales left their feeding grounds in Antarctica only a few weeks before the sample collection of Pirotta et al. 3 and Owen et al. 51 presented evidence of feeding activity of HW off the south coast of New South Wales. Consequently, the whales had likely just entered a fasting state and hence were in their "nutritional prime". In contrast, the comparably low bacterial diversity, richness and number of core taxa of whales later in their migration is correlated to their prolonged fasting. As HW lose 25-50% of their body weight during their annual migration 34,52,53 , and are therefore susceptible to exhaustion of their energy stores before they resume feeding 54, the whales' physiological state and hence their immune system may be compromised. This state may cause a change in the composition of their airway microbiota. Studies on mice and humans showed that the airways did not need to be acutely infected to show changes in their microbial composition 13,14 . Even in healthy or subclinically infected individuals, a shift in immunological parameters resulted in a change of the diversity of microbiota.
The loss of energy stores, metabolic demands on the host and impacts on physiology and immune status due to prolonged fasting may have become increasingly relevant in recent decades. While the population of East Australian HW is growing by almost 11% per year and supposedly approaching their carrying capacity 55, the population density of the whales' main prey, Antarctic krill, are decreasing. Atkinson et al. 56 estimated a 70% decrease in krill numbers since the 1970s. Growing whale numbers together with shrinking prey density could potentially exacerbate the whales' efforts to renew their energy stores in the Antarctic feeding grounds, leaving them more vulnerable to the consequences of fasting-related exhaustion in the future.
The number and abundance of core zOTUs and genera was higher in the blow of whales at the beginning compared to the later stage of their fasting (Tables 3, 4). This observation is correlated to the reduced diversity and richness of the whales in the late fasting stage. Apprill et al. 2 did an analysis of core taxa in the blow of Northern hemisphere HW and found 25 core maximum entropy distribution (MED) sequences in 100% of samples. While they did not look at 80% core microbiota thresholds and their MED sequences are not directly comparable to the zOTUs produced in this study, Apprill et al. 2 likely sampled non-fasting individuals. They performed sample collection in the North Atlantic north of Cape Cod, Massachusetts, USA, in July, and in the North Pacific near Vancouver Island, British Columbia, Canada, in August and off the coast of Washington State, USA, in September. These regions are part of the high-latitude summer feeding areas of the HW that reach up to the Arctic 57 . In autumn, the whales start making their way to tropical waters near the Caribbean 57 . Consequently, sample collection by Apprill et al. 2 took place in the feeding grounds of the whales right before the beginning of the migration to their breeding grounds and when energy stores were expected to have filled up.  www.nature.com/scientificreports/ We found 13 zOTUs in the blow of at least 80% of whales at the beginning of their fasting that were associated with seawater. A similar finding was described in Apprill et al. 2 where three out of 25 MED sequences were related to seawater. It is possible that these zOTUs were mixed in with the blow when the whale exhaled and thus were contaminants. However, according to Apprill et al. 2 seawater regularly enters the airways of whales between breaths and thus may represent 'seawater lavages of the upper respiratory tract seeded with condensed exhalation'. Therefore, these seawater-related zOTUs may also be true inhabitants of the airways, either of transient or more permanent nature. However, as we were not able to distinguish between contamination and 'seawater lavage' at this point, we decided not to include these 13 zOTUs into our conservative estimate of the core.
In the blow of HW at the early and late stage of their migration, we detected several genera that include members previously reported as marine mammal pathogens (Table 5). Yet, the significance of these findings is unclear. Most bacterial genera commonly contain non-pathogenic 'species' alongside pathogenic ones. As 16S rRNA sequencing as used in this study rarely allows the assignment of 'species' identity, we are unable to estimate the true pathogenic capacity of the detected zOTUs. In addition, both groups of whales harboured a large percentage of zOTUs that belonged to unclassified genera (74% in the whales at the beginning of their migration and 55% in the whales later on). The pathogenic potential of these unclassified genera is completely unknown. Therefore, we cannot determine if one group of whales carried a larger ratio of potentially pathogenic bacteria than the other.
In conclusion our study provides some evidence that migration-associated fasting and physiological stress might correlate with a general shift and loss of diversity and richness in the airway microbiota of HW. Such a change in the bacterial community could be the direct result of a weakened immune system, which in turn may influence the rate of elimination of bacteria as well as the growth conditions within the airways enabling certain bacterial groups to proliferate. The analysis of the diversity and composition of HW blow might therefore represent a viable approach to assess the physiological state of cetaceans. In the future, the use of additional techniques to determine body condition including photogrammetry [58][59][60] , measuring cortisol levels 61,62 from the whales' blow and blubber and determining adipocyte size and number in external blubber 63 will contribute to even better understand the correlation between blow microbiota and the physiological state of the whale. Table 3. Taxonomy, relative abundance and details of environment of core zOTUs in the blow of 20 humpback whales at the beginning of their fasting (HW-NM) and 20 humpback whales at a later stage (HW-SM), listed in descending order of the zOTUs' relative abundance. An 80% threshold (core taxa present in 80 of individuals in each group) was used. Whereas HumpbackNM samples showed a large number of core zOTUs, HumpbackSM had none that were not associated with seawater.     www.nature.com/scientificreports/ We collected the samples from the vessel of a commercial whale watching operator at Hervey Bay (25.2882° S, 152.7677° E), Queensland, Australia, in August 2017 (n = 20). The sampled HW's were of unknown sex and age. Images were taken of each sampled whale to compare significant features between sampled individuals (e.g. wounds, scars and pigmentation) to ensure that each whale was only sampled once. The blow collection device followed the design of Hogg et al. 64 and Acevedo-Whitehouse et al. 7 . For each sampling effort, six blank sterile petri dishes (14 mm diameter) were attached to the Perspex plate and a long pole. When a whale surfaced within the range of the pole near the boat, the Perspex plate was positioned about 50 cm above the blowhole of the exhaling whale to collect the droplets in the blow. Seven seawater samples were obtained by filtering 500 ml sea water collected at the blow sampling sites through a Sterivex filter unit (0.22 µm, EMD Millipore Corporation, Billerica, USA). To obtain air samples, the sampling device was exposed to air for 30 min in the absence of whale blow. The petri dishes were subsequently swabbed and then processed in the same way as the blow samples. The air samples acted as technical controls during the DNA extraction process. After sample collection the petri dishes were removed from the Perspex plate and immediately swabbed with sterile cotton swabs (Interpath Services PTY LTD, Heidelberg, Baden-Württemberg, Germany). Each cotton swab had an absorption capacity of about 70 μl of sample material and per sample one cotton swab was used in the downstream process. Face masks were worn at all times while handling the samples to avoid contaminations. The cotton swabs were stored in sterile tubes and chilled on ice until the return of the boat to the harbour. In the field lab, the shaft of the cotton swabs was trimmed and the tip transferred into a 2 ml cryovial. About 300 µl of TE-buffer (10 mM Tris-Cl pH 7.5; 1 mM EDTA) were added to each cryovial to stabilize bacterial DNA. The samples were then stored at-8 °C during transport and then at − 20 °C until processing. All sequences (the HumpbackNM blow and air samples collected by Pirotta et al. 3, the SeawaterNM samples Pirotta et al. 3 used and the HumpbackSM blow, air and seawater samples collected in this study) were analysed together. An initial quality check was performed with FastQC 67 . Paired-end reads were processed with USEARCH (version 10.240) 68 . The reads were merged and low-quality sequences (maximum number of expected errors > 2 and more than 1 ambiguous base) and those shorter than 440 bp were removed. Primers were also removed. Processed sequences of all samples were dereplicated and unique sequences were denoised and de novo clustered into zero-radius operational taxonomic units (zOTUs) with 100% similarity using the unoise3 algorithm 69 . This is the same clustering approach utilized for amplicon sequence variants (ASV) 70 . The difference between zOTUs and ASVs does not lie in the applied similarity-threshold across sequences but in the bioinformatics program. While Callahan et al. 70 used the program DADA2 for sequence processing and 'species' assemblage, we worked with USEARCH 68 .
A de novo chimera removal was included in the unoise step. Afterwards, remaining chimeric sequences were removed using the uchime2 algorithm 71 in high confidence mode with the SILVA database (version 132) as reference dataset 72 . Subsequently, processed sequences were mapped onto zOTU sequences to calculate the presence and relative abundance of each zOTU in every sample using the otutab command with maxrejects and maxaccepts options disabled. Representative zOTU sequences were assigned a taxonomy using the SILVA rRNA sequence database, release 132 (www.arb-silva .de) 72 and the ribosomal database project (rdp), release 11 (https ://rdp.cme. msu.edu) 73 . The taxonomical information derived from the SILVA rRNA sequence database was exclusively used to identify and delete zOTU sequences derived from mitochondria and chloroplasts.
Data analysis. zOTUs that were identified as Archaea, chloroplasts or mitochondria were deleted from the dataset. We identified 457 zOTUs that were present in at least two air samples. As the air samples served as technical control, we considered those zOTUs as technical contaminants and deleted them from the dataset. We used the package phyloseq (v1.24.2) to perform a rarefaction analysis to test if a complete representation of the blow and seawater microbiota was achieved given the observed sequence sampling depths. Bacterial alpha and beta Scientific RepoRtS | (2020) 10:12645 | https://doi.org/10.1038/s41598-020-69602-x www.nature.com/scientificreports/ diversity of seawater and whale samples were assessed using the package vegan (v2.5-5) for community ecology analysis 74 . To determine alpha diversity, we calculated richness, Shannon-Wiener diversity index, Chao1 and ACE species estimator for whale and seawater samples using rarefied counts (lowest number of reads: 3432) to account for the difference in sampling depth (Table 1) and unrarefied counts (Supplementary Table S2). To visualize beta diversity, we used the Bray-Curtis dissimilarity coefficient and weighted and unweighted generalized UNIFRAC distances 75 of unrarefied relative abundances. To determine if the composition of microbial communities of the two groups of whales and their corresponding seawater samples were significantly different from each other, a typical approach is to use PERMANOVA, the Permutational Multivariate Analysis of Variance 76 . However, the implicit assumptions made by distance-based analyses such as PERMANOVA about the meanvariance relationship of the outcome are often unrealistic in community composition data. As demonstrated by Warton et al. 77, this misspecification can lead to the confounding of location and dispersion effects, and therefore incorrect conclusions. For this reason, we instead followed the approach recommended by Warton et al. 77 and used mvabund (v4.0.1) 78 to fit negative binomial regression models to each zOTU, with the log of total sequence counts per sample included as an offset. The sum of likelihood ratio statistics was used as a community-level statistic to compare models with 'migration' and 'whale blow/seawater' as an explanatory factor to an interceptonly model. Using mvabund, the sum of likelihood ratio statistics and statistical significance was evaluated with anova.manyglm using pit-trap resampling 79 . Further details are provided in the supplementary materials (S3, S4). Lokmer et al. 80 applied a similar approach for the analysis of bacterial community composition in oysters. Statistical analysis of microbial community results was performed using R statistical software (v 3.5.1) (https :// cran.r-proje ct.org/).

Data availability.
Sequence data of the blow (HumpbackSM), seawater (SeawaterSM) and the according air samples are available in the NCBI Sequence Read Archive under BioProject accession no. PRJNA521078.