Samples were collected from two adjacent Antarctic fur seal breeding colonies [Freshwater beach (FWB) and the Special study beach (SSB), Figure 1] at Bird Island, South Georgia (54°00'024.8 S, 38°03'04.1 W) during the austral summer of 2018–2019. Ten breeding females and their pups were captured from each colony. Methods used to restrain, and sample captured individuals followed protocols that have been established and refined as part of the long-term monitoring and survey program of the British Antarctic Survey (BAS). Each mother and her pup were captured concurrently around 2–3 days postpartum (in December). The mothers were captured with a noosing pole and held on a restraint board, while the pups were captured with a slip noose and restrained by hand. Biometric data (weight, length, span, and girth) were gathered from the mothers and their pups; the pups were also sexed.

Two microbial samples were collected from each individual. Skin microbial communities were sampled by rubbing a flocked swab (FLOQSwabs 552C, Copan Italia S.p.A.) across the animal's cheek, underneath the eye and behind the snout, while gut microbial communities were sampled by swabbing inside the anus. The resulting microbial samples (n = 80) were individually stored in sterile transport tubes at −20°C during the field season and subsequently at −80°C in the laboratory. Additionally, environmental controls (i.e., swabs from the soil and air, as well as from the gloves and skin of the animal handlers) were collected (n = 6) with the aim of removing potentially contaminating sequences from our samples. A full list of the samples that were collected and analyzed in this study is available in Supplementary Table S1.

Microbial DNA was extracted from all 86 swabs using a NucleoSpin DNA Forensic Isolation Kit (Macherey-Nagel, Düren, Germany). We followed the manufacturer's recommended protocol with the exception that we eluted the DNA in 45 μL of FOE buffer that was preheated to 70°C. The library preparation and 16s rRNA amplicon sequencing were conducted by Biomarker Technologies (BMK) GmbH. A ~247-bp region of the V3–V4 region of the 16S rRNA gene was amplified using the primers 338F (5′ – ACTCCTACGGGAGGCAGCA – 3′) and 806R (5′ – GGACTACHVGGGTWTCTAAT – 3′). The polymerase chain reaction (PCR) mixture contained 2.5–4 ng of DNA, 0.3 μl of each primer (diluted to 10 μM), 5 μL of KOD FX Neo Buffer (Beijing Biolink Biotechnology Co., Ltd), 2 μL dNTPs, 0.2 μL KOD FX Neo (Beijing Biolink Biotechnology Co., Ltd), and ddH₂0 up to a total volume of 10 μL. The PCR amplification program comprised an initial denaturation step of 95°C for 5 min followed by 25 cycles of denaturation at 95°C for 30 s, annealing at 50°C for 30 s, extension at 72°C for 40 s, and a final extension step at 72°C for 7 min. The PCR products were then purified using VAHTSTM DNA Clean Beads (Nanjing Vazyme Biotechnology Co., Ltd.) and amplified in a 20 μL indexing PCR consisting of 5 μL of purified PCR product, 2.5 μL of each indexing primer (Mpp1-n1, Mpp1-n2, diluted to 10 μM; Suzhou SYNBIO Technology Co., Ltd.) and 10 μL Q5 High Fidelity 2x Master Mix. The indexing PCR amplification program comprised an initial denaturation of 98°C for 30 s followed by 10 cycles of 98°C for 10 s, 65°C for 30 s, 72°C for 30 s, and a final extension step of 72°C for 5 min. PCR products were run on a 1.8% agarose gel and quantified with the ImageJ software. 50 ng of each sample was then pooled, purified with EZNA^® Cycle Pure Kit (Omega; Beijing Hongyue Innovation Technology Co., Ltd.) and size-selected by gel electrophoresis with the Monarch DNA Gel Extraction Kit (Beijing Hongyue Innovation Technology Co., Ltd.). The library was then paired-end sequenced on an Illumina Novaseq 6,000 sequencing platform.

Paired-end sequences (247 bp) were analyzed on the Galaxy web platform (http://www.galaxy.org; Afgan et al., 2016). First, adapters were trimmed using Cutadapt (Martin, 2011) and only reads with a minimum read length of 50 were retained. Quality control was performed using FastQC (Andrews, 2010). All of the samples contained adequate numbers (>10,000) of high-quality reads, and we therefore retained them all for subsequent analyses. DADA2 (Galaxy version 1.28; Callahan et al., 2016) was used for further filtering and processing the raw sequences into amplicon sequence variants (ASVs) after trimming 20 bp from the 5′ end of each read and setting the expected error to three for the forward reads and four for the reverse reads. Taxonomy was assigned to the ASVs using the SILVA database version 138.1 (Quast et al., 2013) and the phyla names were then corrected according to Oren and Garrity (2021). Using the decontam package version 1.18.0 (Davis et al., 2018) in R version 4.3.1 (R Core Team, 2022) with the “prevalence” method setting the probability threshold to 0.1, a total of 421 potential contaminant ASVs were removed for the skin and 699 for the gut. In addition, ASVs classified as mitochondrial or chloroplastic, which could not be identified at the phylum level, or which were present in only one sample (singletons) were discarded. Microbial community composition and the core microbiota were characterized for both the skin and gut samples, with the latter defined as including only ASVs present in at least 90% of the samples. This analysis was performed at the phylum level based on non-normalized ASVs. The resulting ASV tables were subsequently analyzed in R version 4.3.1 (R Core Team, 2022).

We calculated microbial alpha diversity using the Shannon diversity index, which considers both the richness and evenness of the taxa present (Allaby, 2014). Shannon diversity was computed using the phyloseq R package version 1.42.0. To control for variation in sequencing depth among the samples, the reads were rarefied to the lowest read depth, following observations of rarefaction curves (Supplementary Figure S1). However, the resulting alpha diversity metrics were strongly positively correlated with those based on the non-rarefied dataset (Spearman's rank correlation: r = 0.99, p < 0.001). Consequently, in order to maximize the available data analyzed and to avoid excluding rare ASVs, we used the non-rarefied dataset for subsequent analysis.

We aimed to investigate the effects of body site (skin vs. gut), colony (FWB vs. SSB), and age class (mother vs. pup) on alpha diversity. First, we compared the Shannon diversity index between skin and gut samples, finding them to be significantly different (Kruska-Wallis: χ(1)2 = 47.4, p < 0.001, Supplementary Figure S2). We then fitted linear mixed effect models (LMMs) with the lmerTest R package version 3.1–3 (Kuznetsova et al., 2017); separate models for the skin and gut microbial communities were run. Each model included the Shannon index as the response variable and colony and age class as fixed effects. Mother-pup pair ID was included as a random effect to control for the spatial proximity of mothers and their offspring. To explore the effect of sex on alpha diversity, we also ran linear models (LM) only for pups, which included both female and male observations. For all models, the significance of the fixed effects was assessed at a = 0.05 by running Type II and Type I analysis of variance (for the LMs and LMMs, respectively), as well as by evaluating whether the 95% confidence intervals (95% CIs) of the effect size estimates overlapped zero. The significance of the random effects was assessed using the ranova function of the lmerTest R package. The variance explained by the fixed effects (Rmarginal2) and by both the fixed and random effects (Rconditional2) was calculated using the r.squaredGLMM function of the MuMIn package version 1.47.5 (Barton, 2009). The model residuals were inspected for normality and homogeneity of variance using the performance package version 0.12.0 (Lüdecke et al., 2021).

Before testing for differences in microbial composition (beta diversity) among the samples, we applied cumulative sum scaling (CSS) normalization to our filtered dataset using the metagenomeseq package version 1.30.0 (Paulson et al., 2013). We then calculated pairwise microbial dissimilarity matrices based on Bray-Curtis dissimilarities and visualized differences among the groups using non-metric multidimensional scaling (nMDS) plots. To statistically assess whether any differences in microbial composition could be related to colony, age class or mother-pup pair, we implemented permutational multivariate analyses of variance (PERMANOVA; Anderson, 2017) with 1,000 permutations using the adonis function of the vegan R package version 2.6–4 (Oksanen et al., 2017). To evaluate whether significant results from the PERMANOVA arose due to variation in dispersion within groups instead of from differences in the mean, we performed homogeneity of group dispersion tests using betadisper in the vegan R package (Oksanen et al., 2017).

Finally, we identified differentially abundant ASVs between colonies and age classes by analyzing the filtered and non-normalized datasets using the corncob package version 0.4.1 (Martin et al., 2020). This software uses beta-binomial regressions to model the relative abundances of all the taxa while taking into account differential variability. Significance testing was performed using likelihood ratio tests and the p-values were obtained after controlling for the Type I error using Bonferroni correction (p < 0.05).

We identified 30 distinct microbial phyla for the skin and 25 for the gut (Figure 2; Supplementary Table S5). Although both body sites shared the same five most abundant phyla, their relative abundances differed (Pseudomonadota: skin = 30.1%, gut = 18.3%; Bacillota: skin = 26%, gut = 36.2%; Actinomycetota: skin = 17.3%, gut = 8.2%, Bacteroidota: skin = 14.6%, gut = 16%; Fusobacteriota: skin = 8%, gut = 15.6%). Focusing on those taxa present in at least 90% of the samples from each body site, we identified 136 core ASVs in the skin and 64 in the gut (Supplementary Tables S6, S7). Four phyla (Bacillota, Pseudomonadota, Fusobacteriota, and Actinomycetota) dominated the core microbiota of both the skin and the gut, with Pseudomonadota being the most abundant phylum in the skin core microbiota (33%) and Bacillota being the most abundant phylum in the gut core microbiota (42%).

Skin microbial alpha diversity, quantified using the Shannon index, was significantly higher than gut alpha diversity (Kruska-Wallis: χ(1)2 = 47.4, p < 0.001) and there was little overlap in the distributions of Shannon index values between the two body sites (Supplementary Figure S1). For the skin, microbial alpha diversity was significantly lower in the high-density colony (LMM_{Colony = SSB}: β = –0.74, χ² = 17.51, p < 0.001), but did not differ significantly between mothers and pups (LMM_{Ageclass = Pup}: β = −0.24, χ² = 3.08, p = 0.08; Figure 3A; Supplementary Table S8), with the fixed effects explaining almost 40% of the variation (Rmarginal2 = 0.38, Rconditional2 = 0.52). By contrast, gut alpha diversity did not differ significantly between the colonies (LMM_{Colony = SSB}: β = −0.07, χ² = 0.09, p = 0.76) or between mothers and pups (LMM_{Ageclass = Pup}: β = 0.06, χ² = 0.07, p = 0.79; Figure 3C; Supplementary Table S8), with the fixed effects explaining only a small proportion of the total variation (Rmarginal2 = 0.004, Rconditional2 = 0.01). Lastly, no differences in alpha diversity were found between male and female pups for either the skin (LM_{Sex = Male}: β = −0.51, F = 3.36, p = 0.08) or gut (LM_{Sex = Male}: β = −0.02, F = 1e-03, p = 0.98; Supplementary Table S9).

Visualizing microbial composition with nMDS plots based on Bray-Curtis dissimilarities, we observed two distinct clusters corresponding to FWB and SSB for the skin, while clustering was less defined for the gut (Figures 3B, D). This was reflected in the results of PERMANOVAs, which revealed an almost three times higher proportion of the variance explained by colony for the skin (PERMANOVA_Colony: R² = 0.23, p = 0.001) relative to the gut (PERMANOVA_Colony: R² = 0.08, p = 0.001; Supplementary Table S10). Significant differences in microbial composition between the age classes were found for both body sites (PERMANOVA_{Skin − Age class}: R² = 0.04, p = 0.007; PERMANOVA_{Gut − Ageclass}: R² = 0.09, p = 0.001) although this pattern was weaker for the skin (Supplementary Table S10). Additionally, microbial communities showed significantly greater variability among mothers for the skin (betadisper_{Skin − Age class}: F = 11.13, p = 0.002) and among pups for the gut (betadisper_{Gut − Age class}: F = 18.6, p < 0.001; Figures 3B, D).

As previously found by Grosser et al. (2019), mothers and their pups also carried significantly more similar skin microbial communities than expected by chance (PERMANOVA_{Skin − PairID}: R² = 0.45, p = 0.001). However, this pattern was not present in the gut microbial communities (PERMANOVA_{Gut − PairID}: R² = 0.42, p = 0.1). Additionally, no significant differences in microbial composition were found between male and female pups either for the skin or the gut (PERMANOVA_{Skin − Sex}: R² = 0.05, p = 0.23; PERMANOVA_{Gut − Sex}: R² = 0.06, p = 0.28; Supplementary Table S10).

Using beta binomial regression models, we identified differentially abundant ASVs between the two colonies and age classes. A greater number of differentially abundant ASVs was detected for the skin compared to the gut in relation to both colony and age class. For colony comparisons, most of the differentially abundant ASVs were more abundant in the low-density colony (FWB) and corresponded to the phyla Patescibacteria, Gemmatimonadota, Bacteroidota, and Actinomycetota. By contrast, only two ASVs were more abundant in the high-density colony (SSB), belonging to the Fusobacteriota and Bacillota (Figure 4A). Regarding age class, no differentially abundant ASVs were identified for pups, while those taxa that were more abundant in mothers included phyla typically present in diverse terrestrial and marine habitats, as well as the Campylobacterota (Figure 4B).

For the gut, only two ASVs were differentially abundant between the two colonies: an ASV that was more abundant in the high-density colony (SSB) belonged to the Camplylobacterota, a phylum associated with bite wounds in pinnipeds (Foster et al., 2020), while an ASV that was more abundant in the low-density colony (FWB) belonged to the Patescibacteria (Figure 4C). Most of the differentially abundant ASVs were detected in relation to age class for the gut. Only one ASV belonging to the Bacteroidota was significantly more abundant in mothers, with the rest of the pup-associated ASVs corresponding to phyla commonly found in the soil, freshwater and marine environments, as well as the Pseudomonadota, a taxon that is detected in human infants (Yang et al., 2024, Figure 4D).

We found that the skin and gut microbial communities of Antarctic fur seals were distinct, despite sharing five highly abundant phyla (Bacillota, Pseudomonadota, Bacteroidota, Fusobacteriota, and Actinomycetota). More specifically, the phylum Pseudomonadota dominated the skin microbial communities. This was expected as this taxon is generally detected at high abundance on the skin of cetaceans and fur seals, as well as in seawater (Grosser et al., 2019; Van Cise et al., 2020). By contrast, and in line with other studies of Antarctic pinnipeds, including other species of fur seals (Nelson et al., 2013; Medeiros et al., 2016), the gut microbiota showed a high abundance of Bacillota. This taxon, a well-known gut colonizer, is thought to aid the host in efficiently extracting energy from food, which promotes fat deposition, an essential process for pinnipeds that depend heavily on fat storage for energy intake and thermoregulation (Pacheco-Sandoval et al., 2019).

The environment in which a host lives, including both the social and physical environment, can have a major impact on host microbial communities (Tasnim et al., 2017; Callewaert et al., 2020). To investigate the effects of host and environmental factors on the Antarctic fur seal microbiota, we used a unique natural experiment comprising two adjacent breeding colonies that differ in social density but which experience ostensibly the same environmental conditions due to their close physical proximity. Focusing initially on skin microbial communities, we found marked differences between the two colonies, with alpha diversity being significantly lower in the high-density colony. This corroborates the results of a previous study of mother-pup pairs from the same two colonies that used similar approaches but only investigated the skin microbiota (Grosser et al., 2019). As these two studies were conducted 4 years apart and there was no overlap between the individuals sampled at the two timepoints, we conclude that these patterns are reproducible and hence show temporal stability. This finding is important in its own right as replication studies in the fields of ecology and evolution are uncommon (Kelly, 2006; Nakagawa and Parker, 2015; although see Tebbe et al., 2020), reflecting the general perception that research in these fields can be difficult to replicate, partly due to the challenges of controlling for confounding factors such as environmental variation in natural settings.

Clear and repeatable differences in skin microbial alpha diversity between the high- and low-density colonies are consistent with a suppressive effect of social density on microbial diversity (Grosser et al., 2019). In line with this, cortisol levels have been shown to be elevated in mothers and pups from SSB (Meise et al., 2016; Nagel et al., 2022) reflecting a potentially more stressful and competitive environment. A negative association between elevated social stress and alpha diversity has been also reported in captive hamsters (Partrick et al., 2018), while other studies have shown that high stress levels, indicated by increased glucocorticoid concentrations, have a suppressive effect on alpha diversity in humans and red squirrels (Stothart et al., 2016; Keskitalo et al., 2021). Therefore, our results lend further support to the argument that overcrowding increases social stress and reduces skin microbial alpha diversity.

By contrast, gut alpha diversity did not differ significantly between FWB and SSB. Even though studies of the impacts of social density and stress on mammalian gut microbiomes are scarce, a similar pattern was observed in a recent study of laboratory mice in which social overcrowding elevated corticosterone levels without affecting gut alpha diversity (Delaroque et al., 2021). Similarly, studies of great apes (Vlčková et al., 2018; Hickmott et al., 2022) did not find any significant effects of stress on gut microbial diversity, suggesting that these species may buffer stress through social interactions and behavioral flexibility (Hickmott et al., 2022). Overall, the lack of an obvious relationship between proxy measures of social stress and gut alpha diversity in Antarctic fur seals and other species may be a reflection of the host's ability to regulate its microbiota, maintaining homeostasis and protecting against harmful colonizers (Foster et al., 2017). In particular, gut epithelial cells play a central role by secreting antimicrobial peptides and limiting oxygen availability, thereby fostering a microbiota that supports host health and nutrition (Dabrowska and Witkiewicz, 2016).

To further investigate community-level differences, we analyzed the beta diversity of the two breeding colonies across both body sites. Consistent with our results for alpha diversity, we found pronounced differences in microbiota composition between the colonies for the skin, while these differences were considerably weaker for the gut. This alignment between alpha and beta diversity is not surprising as the composition and diversity of microbial communities are often closely linked (Grosser et al., 2019). This linkage suggests that similar regulating factors may influence both aspects of microbial community structure.

Our dataset also enabled us to test for differences in microbial communities between mothers and pups. We found no significant differences in alpha diversity between the two age classes for either the skin or the gut. Similarly, while mothers and pups exhibited significant differences in microbial community composition, these patterns were relatively weak, particularly for the skin. These findings are consistent with a previous study of Antarctic fur seals that also observed minimal differences in skin microbial diversity and composition between mothers and pups from the same colonies (Grosser et al., 2019). However, they contrast with findings from studies of other pinnipeds, such as harbor seals and southern elephant seals, where pronounced differences in the gut microbiota of adult females and pups were observed (Nelson et al., 2013; Pacheco-Sandoval et al., 2022). The similarities between Antarctic fur seal mothers and newborn pups may result from the vertical transmission of microbes at birth, during nursing and via close physical proximity, which also exposes them to similar environmental microbes (Coates et al., 2019). By contrast, Pacheco-Sandoval et al. (2022) observed pronounced differences in the gut microbiota of harbor seal mothers and pups, likely due to diet-driven changes as the pups were transitioning from milk to solid foods. To better understand the colonization and development of microbial communities during early life-stages, longitudinal studies of Antarctic fur seals pups from birth to weaning will be needed.

Additionally, we found that mother-pup pairs shared more similar microbial communities than expected by chance for the skin (as previously reported by Grosser et al., 2019) but not for the gut. As argued above, these similarities in skin microbial communities are to be expected given the close proximity of mothers and their newborn pups during activities such as suckling and resting, which can facilitate the transmission of skin microbes (Coates et al., 2019). By contrast, the development of the gut microbiota may be more strongly influenced by factors such as dietary shifts, including the transition from milk to solid food, as well as by the maturation of the immune system (Coates et al., 2019; Reese et al., 2021).

To gain deeper insights into the patterns described above, we conducted differential abundance analyses. Focusing on the effects of social density, we found that several taxa typically associated with pathogenesis were enriched in both the skin and gut microbiota of individuals from the high-density colony (SSB). Notably, Bacillota were more abundant in the skin microbiota of animals from SSB, a taxon previously linked to skin wounds and disease progression in stranded porpoises (Li et al., 2022). Overcrowding and social stress are likely contributing factors, as they are known to reduce microbial alpha diversity and induce shifts in microbial abundance that can lead to dysbiosis and the colonization of pathogenic taxa (Alverdy and Luo, 2017; Delaroque et al., 2021).

Interestingly, despite the absence of significant differences in gut alpha diversity between the colonies and only a weak influence of breeding colony on gut beta diversity, we found a significantly higher abundance of Campylobacterota in the gut microbiota of individuals from the high-density colony (SSB). This phylum has been linked to gastrointestinal diseases in humans and to pathogenesis associated with bite wounds in pinnipeds (Foster et al., 2020; Watkins et al., 2022). This observation lends further support to the argument that, while animals from both colonies maintain diverse microbiota, individuals from the high-density colony appear to experience an enrichment of pathogenic taxa. To better understand the effects of social stress on the gut microbiome, future studies should incorporate measurements of stress hormones such as cortisol to explore the relationship between stress and the gut microbiota at the individual level. Additionally, metagenomic approaches could be used to investigate the functional characteristics of differentially abundant microbes and their roles in host health and disease.

Above and beyond these patterns, we also identified several differentially abundant taxa at the low-density colony (FWB) belonging to phyla commonly found in soil, freshwater and marine habitats (i.e., Patescibacteria, Actinomycetota, Gemmatimonadota, Anandan et al., 2016; Mujakić et al., 2022; Haro-Moreno et al., 2023). The presence of a small stream at FWB likely facilitates the acquisition of these environmental microbes, which appear to be an important source of variability in microbial diversity and composition at the low-density colony. This finding aligns with a previous study of harbor seals, where populations near freshwater sources were found to exhibit higher alpha diversity (Pacheco-Sandoval et al., 2019). It emphasizes the importance of environmental sources of variation, particularly in the absence of external stressors such as overcrowding, and further supports the idea that host habitat plays a key role in shaping microbial community composition in marine mammals, alongside diet, and phylogeny (Bik et al., 2016).

Finally, regarding age class, we found that the gut microbiota of pups was dominated by the phyla Pseudomonadota, while the gut microbiota of mothers was dominated by Bacteroidota. This pattern has been previously described in human mothers and their infants (Suárez-Martínez et al., 2023; Yang et al., 2024). It reflects the natural progression of gut colonization during development, with facultative anaerobes like Pseudomonadota thriving in the oxygen-rich environment of the newborn gut. As individuals transition from milk to solid foods during weaning, the gut microbiota becomes dominated by phyla like Bacteroidota, which are able to break down complex polysaccharides (Suárez-Martínez et al., 2023). We also found that the gut microbiota of pups was enriched for environmental phyla such as Patescibacteria, Deinococcota, and Acidobacteriota, suggesting that many early colonizers originate from the environment where the individuals are born (Schmidt et al., 2019; Suárez-Martínez et al., 2023; Yang et al., 2024).