During the austral summer 2016, four sponge species [Dendrilla antarctica (Topsent, 1905), Sphaerotylus antarcticus (Kirkpatrick, 1907), Mycale (Oxymycale) acerata (Kirkpatrick, 1907), and Hemigellius pilosus (Kirkpatrick, 1907)] were collected at different locations from two South Shetland Islands and the Antarctic Peninsula (Table 1). Replicate seawater samples (n = 3, 2 l samples) were aseptically collected adjacent to the sampled sponges from all locations. Sponges were processed after collection. A sample from each specimen was taken with a sterile scalpel and rinsed several times in 0.22 μm-filtered seawater to discard loose attached microorganisms. Seawater samples were sequentially passed through polycarbonate 5 and 0.22 μm filters (MilliporeSigma, Burlington, MA, United States) to avoid clogging the ones with the smaller pore diameter. The contents on the 0.22 μm filters were used to examine the ambient bacterioplankton communities. All samples were preserved in RNAlater until further use.

DNA extracts were submitted to Molecular Research LP¹ (Shallowater, TX, United States) for amplification, library construction and multiplexed sequencing of partial (V4) 16S rRNA gene sequences on an Illumina MiSeq platform. The HotStarTaq Plus Master Mix kit (Qiagen) was used for PCR amplifications using DNA extracts as templates with universal prokaryotic 515F (Parada et al., 2016) and 806R (Apprill et al., 2015). To barcode samples, a multiplex identifier barcode was attached to the forward primer. The thermocycler profile consisted of an initial denaturation step at 94°C for 3 min; 28 cycles of 94°C for 30 s, 53°C for 40 s, and 72°C for 1 min with a final elongation step at 72°C for 5 min. Equimolar concentrations of samples (Table 1) were pooled and purified using Agencourt AMPure XP beads (Beckman Coulter) to prepare DNA library by following Illumina TruSeq DNA library preparation protocol. Sequencing was then performed according to manufacturer’s guidelines on an Illumina MiSeq.

As described previously (Thomas et al., 2016), Illumina sequence reads were processed in mothur v1.39.5 (Schloss et al., 2009). Briefly, raw reads were demultiplexed, forward and reverse reads were then joined, and sequences <200 bp and/or with ambiguous base calls were removed. Sequences were aligned to the SILVA database (release 128, non-redundant, mothur-formatted), trimmed to the V4 region, and screened for chimeras and errors. A naïve Bayesian classifier and Greengenes taxonomy (August 2013 release, mothur-formatted) was used to aid in the removal of non-target sequences (e.g., chloroplasts, mitochondria). We used the SILVA database (release 132, non-redundant, mothur-formatted) for final taxonomic assignment. The resulting high-quality sequences were clustered into operational taxonomic units (OTUs) defined by clustering at 3% divergence and singletons were removed. We used rarefaction curves (mothur v1.39.5) to plot the OTUs observed as a function of sequencing depth. To avoid artifacts of varied sampling depth on subsequent diversity calculations, each sequence dataset was subsampled to the lowest read count (mothur v1.39.5). To place the determined OTUs into a greater context, we performed BLAST searches (NCBI-BLAST-2.7.1+) of these sequences against the database of the sponge EMP project (Moitinho-Silva et al., 2017a).

To compare prokaryotic community profiles, we constructed a Venn diagram² from presence/absence OTU data. Non-metric multidimensional scaling (nMDS) plots of Bray-Curtis similarity matrices were also constructed with mothur (v1.39.5) and R (version 3.4.3; ggplot2 package) from square-root transformed OTU relative abundance data. We also constructed bubble charts in R (version 3.4.3; ggplot2 package) from OTU relative abundances to plot community dissimilarities among locations. Significant differences among sponge species and ambient seawater were assessed using a one-way permutational multivariate analysis of variance (PERMANOVA), with the factor source (all sponge species vs. seawater). Significant differences among sponge species were assessed with a one-way PERMANOVA, with the factor source (D. antarctica, S. antarcticus, M. acerata, and H. pilosus). Differences between sponge species and locations were assessed using a two-way PERMANOVA, with the factors source (D. antarctica and S. antarcticus), location (Deception Island, Rothera, Half Moon Island, O’Higgins), and an interaction term. Pairwise comparisons were subsequently conducted for all significant PERMANOVA results. Permutational multivariate analysis of dispersion (PERMDISP) was used to detect differences in homogeneity (dispersion) among groups for all significant PERMANOVA outcomes. All multivariate statistics were performed in R (version 3.4.3; with adonis2 and betadisper functions from vegan v2.5-6 package).

We calculated three indices of alpha diversity in mothur v1.39.5 (Schloss et al., 2009) to evaluate community richness, diversity, and evenness: observed OTU richness, the inverse of Simpson index of diversity, and the Simpson index of evenness. One-way analyses of variance (ANOVA) was used to detect differences in diversity metrics among the species from Deception Island (D. antarctica, S. antarcticus, M. acerata, and H. pilosus). Two-way analysis of variance (ANOVA) was used to detect differences in diversity metrics by the factors source (D. antarctica and S. antarcticus), location (sampling sites) and an interaction term, followed by pairwise comparisons for any significant factor. All data that did not meet the statistical assumptions was transformed accordingly (log-transformation for inverse of Simpson index). The univariate statistics were performed in R (version 3.4.3; ANOVA function from car package).

We were interested in OTUs that were abundant (i.e., >0.1% relative abundance) and widespread among sponge host individuals (i.e., 90% prevalence), so we performed a core microbiome (i.e., with relative abundances >0.1% and present in all species) analysis in R (version 3.5.3, package Microbiome). We also analyzed the dataset for patterns in relative abundances of OTUs within categories (e.g., sponge vs. seawater, among locations). For this purpose, we removed from the dataset rare OTUs (<0.1% relative abundance) and OTUs with a low incidence across samples (detected in <2 samples). We used the Mann-Whitney-U test (or Wilcoxon rank sum test) with FDR p-value correction to identify significantly different patterns in OTU relative abundance among hosts and locations using QIIME 1 (Caporaso et al., 2010).

The V4 region of the 16S rRNA gene was sequenced on an Illumina MiSeq platform and a total of 4,398,237 reads were obtained after denoising and quality filtering with a library depth ranging from 39,883 to 154,480 reads. As we had four replicates per species and location in most of the cases, we discarded those samples (n = 5) with the lowest number of reads (≤51,407) to have at least three replicates per sponge and site. The only exception was S. antarcticus from O’Higgins Station, from which we had two replicates. This may represent low replication, however, those samples clustered together with most samples of the same species and we reached the same results without the O’Higgins location. Therefore, we decided to include them in further analyses. To avoid artifacts of varied sampling depth, we rarefied our libraries to the lowest read count after removing the previous samples from the dataset (n = 52,637; Supplementary Figure S1). Twenty-eight bacterial and three archaeal phyla were detected in the 11,187 OTUs recovered from seawater and sponge samples, which were predominantly affiliated to the phyla Proteobacteria and Bacteroidetes (Supplementary Figure S2). Of these, 4,619 OTUs were recovered from D. antarctica, 3,438 OTUs from S. antarcticus, 1,381 OTUs from M. acerata, and 1,490 OTUs from H. pilosus. Seawater exhibited greater richness with 6,071 OTUs, 2,511 of which were shared with either D. antarctica or S. antarcticus or with both species, and 1,240 were shared with M. acerata and/or H. pilosus (Supplementary Figure S3). Seawater from Rothera, and especially from Deception Island contained the lowest amounts of Gammaproteobacteria (Figure 1).

The taxonomic composition of prokaryotic communities recovered from ambient seawater, and from S. antarcticus, M. acerata and H. pilosus sponge hosts, were significantly different, while D. antarctica presented a microbial content quite similar to what we found in seawater (Figure 1). Two more phyla were detected in sponges compared to the surrounding seawater (Supplementary Figures S4, S5). However, if we discard phyla with low sequence abundances (i.e., <0.1% abundance), sponges and seawater harbored five bacterial and one archaeal phyla. The prokaryotic communities harbored by S. antarcticus, M. acerata, and H. pilosus were enriched for Gammaproteobacteria (>60, >80, >65% of the prokaryotic community on average, respectively) and Thaumarchaeota (>20% in S. antarcticus and H. pilosus, >8% in M. acerata) compared to seawater (<42 and <0.5%, respectively). Comparatively, D. antarctica hosted less Gammaproteobacteria (<54%) and its associated Thaumarchaeota were almost absent (<0.5%). Contrastingly, prokaryotic communities in the first three sponges were depleted in members of Alphaproteobacteria (<5% in all three species) and Bacteroidetes (<8% in all three species) compared to seawater (>24 and >30%, respectively) and D. antarctica (>26 and >17%, respectively). The inter-individual variability of the taxonomic composition depends on the host. While H. pilosus and S. antarcticus harbored a quite stable microbial signature, D. antarctica exhibited greater inter-individual variability and one of the specimens of M. acerata showed an enrichment for Thaumarchaeota (Figure 1).

Statistically significant differences in community structure (PERMANOVA) were detected among S. antarcticus, M. acerata, H. pilosus and D. antarctica and seawater microbes (F_4,39 = 10.563; P = 0.001). The four species significantly differed from each other (Pairwise comparisons; P < 0.004) except for M. acerata and H. pilosus (P = 0.100). Symbiont communities from seawater sources exhibited no overlap with sponge species in the multidimensional space, and all sponge species occupied distinct regions of the nMDS plot (Figure 2). In addition, a significant interaction between host species (S. antarcticus and D. antarctica) and location occurred (PERMANOVA, F_3,18 = 2.422; P = 0.010), though we could not detect significant pairwise differences in community structure after p-value correction. Dispersion analysis revealed equal variability within S. antarcticus and D. antarctica prokaryotic communities regardless of location (P > 0.05 in all comparisons), but microbiomes of the latter species from O’Higgins were more variable (P = 0.046; Figure 2).

Larger mean values of richness, diversity (i.e., inverse Simpson diversity index), and evenness in symbiont communities from seawater compared to host species were observed (P < 0.001 in all pairwise comparisons; Table 2). When we analyzed the sponges from Deception Island, the microbiome of all the species seemed to be equally richer and diverse with similar evenness except for the comparison between H. pilosus and M. acerata, where the former species presented more OTU richness (P = 0.048). Comparing D. antarctica and S. antarcticus from the four locations studied, a two-way ANOVA detected significant differences between hosts and locations for species richness, diversity and evenness (P < 0.03 in all cases). D. antarctica harbored a richer microbiome but less diverse than S. antarcticus. Additionally, the sponge prokaryotic communities from Half-moon Island had greater species richness and diversity compared to those from Rothera (pairwise comparisons P < 0.05). Although there was an effect of location on the microbiome evenness, differences among pairs were not detected.

The abundance of shared OTUs (n = 2,893) between sponge-associated and seawater prokaryotic communities was calculated and just 2.2% presented relative abundances over 0.1%. Those few OTUs (n = 65) accounted for 93 and 89% of the total relative abundance of sponge-associated and seawater prokaryotic communities, respectively, which meant that sponge-specific OTUs (64%) fell within the “rare biosphere” (<0.1% relative abundance).

In addition to community-level metrics of diversity and structure, we performed a core microbiome analysis to investigate patterns in abundant and prevalent individual OTUs among sponge hosts. The strict core microbiome of the sponge hosts was formed by seven OTUs accounting for 50% of total relative abundance on average (0.1% of total OTUs present in sponge hosts, 4.2% of OTUs with relative abundance >0.1% in at least one sample; Supplementary Table S1). Significant sponge enrichments in four core OTUs, affiliated to Gammaproteobacteria, were detected with a mean fold-change in abundance of 44.9 ± 8.0 (±SE) with respect to seawater (mean relative abundance 0.28%). Three additional sponge core OTUs, which were affiliated to the groups Alphaproteobacteria, Bacteroidetes, and Gammaproteobacteria, were more abundant in seawater communities (fold-change 14.3 ± 4.5; Supplementary Table S2). We also determined the variable community (i.e., with relative abundances >0.1% and present in at least two species) formed by 56 OTUs that represented on average 39% of the total abundance (0.7% of total OTUs present in sponges, 33.3% of OTUs with a relative abundance >0.1% in at least one sample; Supplementary Table S1). Forty-three sponge variable OTUs that had a mean fold-change in abundance of 49.5 ± 8.2 were extremely rare in seawater (mean relative abundance <0.03%). Thirteen additional sponge variable OTUs were enriched in seawater with mean relative abundance >2.7% (fold-change 13.7 ± 2.6). We comparatively determined the core microbiome of D. antarctica and S. antarcticus. Both species harbored ten major OTUs (six out of seven sponge core OTUs and four variable OTUs) representing over 70% of the microbiome relative abundance (Supplementary Table S1). Seawater presented instead a core microbiome of 18 OTUs (including six sponge core OTUs) accounting for over 60% of the microbiome in relative abundance (Supplementary Table S1). Comparing the sponge species analyzed, four core and 17 variable OTUs seemed to be sponge-enriched for either D. antarctica (OTUs 1, 4, 9, 11, 13, 27, 33, and 45; cumulative 84% relative abundance; Figure 3), S. antarcticus (OTUs 2, 7, 23, 31, 43, 51, and 56; 54% relative abundance; Figure 3), M. acerata (OTU 6; 65% relative abundance) or H. pilosus (OTUs 3, 14, 19, 34, and 53; 71% relative abundance), although they were also detected in the other hosts and in seawater in lower frequencies (Supplementary Table S2). If we compared locations, neither D. antarctica nor S. antarcticus presented differences in their microbiome abundances among sampling sites (Figure 3 and Supplementary Table S2). However, OTU 2 decreased its relative abundance to 0.3% in S. antarcticus at Half Moon Island, whereas the proportion was over 40% on average for the rest of locations (Figure 3).

From the 8,009 OTUs present in the sponges analyzed, 56.1% could be described as host-specific (i.e., present only in one sponge species), ranging from 254 OTUs in M. acerata to 2,230 OTUs in D. antarctica. However, all of these OTUs belonged to the rare biosphere associated to the sponge host with mean relative abundance below 0.1%. Although OTUs from core and variable prokaryotic communities were detected in all the sponge species analyzed, and even in ambient seawater, the specificity for one of the hosts of some variable OTUs (38%) could be considered due to the low number of reads recovered from the other habitats (<0.02% relative abundance).

Local BLAST searches against the sponge EMP database showed that 87% of the OTUs (n = 9,782) were found among the sponge microbiome collection at high sequence similarities (Supplementary File S1). Both core and variable communities associated to Antarctic sponges had a closest relative in the sponge microbiome database with sequence identities over 97%. Only three variable OTUs had similarities between 93 and 96%. The core microbiome associated to the Antarctic sponges is also associated to other sponge hosts and habitats (Supplementary File S1).

Although there is a lot of diversity to be uncovered, we have captured all abundant microbes present in Antarctic sponges and seawater (Supplementary Figure S1). The sponge associated prokaryotic communities were less diverse than surrounding seawater as previously reported for other Antarctic sponge hosts, reinforcing the view that these sponges were composed of low microbial abundance (LMA) microbiomes (Moitinho-Silva et al., 2017b). Both biotopes presented Thaumarchaeota as the predominant archaeal phylum, which accumulated 10.8% of the reads in sponges while just 0.1% in seawater. With regard to bacteria, differences lay in the fact that sponges hosted bacteria of Nitrospirae and Nitrospinae phyla, while seawater had instead representatives of Cyanobacteria and Verrucomicrobia. However, all these phyla accumulated microbial abundances ranging from 0.1 to 1.4%.

Classes Gamma- and Alpha-proteobacteria dominate bacterial assemblages in association with marine sponges from different biotopes (Thomas et al., 2016; Moitinho-Silva et al., 2017a; Pita et al., 2018; Cleary et al., 2019). These classes also dominate Antarctic sponges as previously reported on Antarctic marine shallow (Webster et al., 2004; Rodríguez-Marconi et al., 2015; Cárdenas et al., 2018, 2019) and deep waters (Steinert et al., 2019) with Bacteroidetes also as important members in terms of number of OTUs recovered (n = 1,174) and in relative abundance (10%). The presence of the phylum Nitrospirae (>3% in relative abundance) in one of the sponges analyzed (S. antarcticus) is particularly noticeable. Members of this phylum are potentially involved in nitrification processes, specifically in nitrite oxidation (Radax et al., 2012). A previous study provided a preliminary description of the bacterial communities present in Deception Island waters, which were very different than the usual Antarctic microbial diversity (Angulo-Preckler et al., 2015). We found that Bacteroidetes dominated the seawater prokaryotic community, probably influenced by the special conditions of this island, an active volcano with high concentrations of geothermal elements. This suggested a correlation between the environmental microbial diversity and the geochemical composition of the island. However, the sponges studied from the same area contained a different prokaryotic signature than the environment.

The phylum Chloroflexi, and other bacterial phyla such as Actinobacteria or Acidobacteria, have been frequently found in high microbial abundance (HMA) sponges (Schmitt et al., 2011; Moitinho-Silva et al., 2017b), occasionally reaching high percentages in relative abundances (5–10%). Contrastingly, these phyla represent a much lower percentage in Antarctic sponges (0.0005–0.4%). Our results are in agreement with previous studies carried out in Antarctica (Rodríguez-Marconi et al., 2015; Cárdenas et al., 2018; Steinert et al., 2019; Moreno-Pino et al., 2020). Other widely described but less abundant phyla associated with sponges from other geographical areas are Poribacteria and PAUC34f (Moitinho-Silva et al., 2017b). These bacterial groups were not detected neither in the Antarctic sponges from the present study nor in other species previously studied in Antarctica (Rodríguez-Marconi et al., 2015; Cárdenas et al., 2018, 2019; Steinert et al., 2019; Moreno-Pino et al., 2020). These phyla have been found overrepresented in HMA over LMA sponges and can be used as “HMA/LMA indicators” (Moitinho-Silva et al., 2017b; Glasl et al., 2018). This suggests that the microbial community structure of the sponges from the present study resemble that from LMA sponges.

Besides Bacteria, the presence of Archaea associated to marine sponges from tropical to cold waters has been previously documented (Radax et al., 2012; Cardoso et al., 2013; Jackson et al., 2013; Kennedy et al., 2014; Polónia et al., 2014; Turon and Uriz, 2020). These archaeal communities were dominated by the phylum Thaumarchaeota, which can reach high relative abundances in the sponge prokaryotic community. Three of the species analyzed in this study (S. antarcticus, M. acerata, H. pilosus) harbored abundant populations of an OTU (∼10–25%) affiliated to Candidatus Nitrosopumilus. In Antarctica, the presence of Thaumarchaeota-related sequences has been reported in multiple species (Steinert et al., 2019; Moreno-Pino et al., 2020), including M. acerata (Webster et al., 2004; Rodríguez-Marconi et al., 2015). These taxa act as ammonia oxidizers not only in tropical reefs but also in cold environments (Radax et al., 2012; Cardoso et al., 2013; Polónia et al., 2015). Indeed, the genomic potential for ammonia oxidation was recently evidenced for Leucetta antarctica microbiome from Antarctica (Moreno-Pino et al., 2020). Together with the potential nitrite oxidizers detected (members of Nitrospirae), they may be considered as relevant factors in nitrogen cycling in Antarctica. However, further studies are needed to understand the functional roles these microorganisms play in Antarctic symbiosis and their contribution to the ecosystem.

A minimal core prokaryotic community (0.1%) was found in Antarctic sponges consisting of very few OTUs (n = 7) present in all four species but also in ambient seawater. The core prokaryotic community of sponges is rather small, pattern that has been documented from different biogeographic regions, including Antarctica (Schmitt et al., 2012; Easson and Thacker, 2014; Rodríguez-Marconi et al., 2015). We also defined a small variable community (0.7%) that consisted of 56 OTUs hosted in at least two of the sponge species analyzed with relative abundances above 0.1%. Nevertheless, 43% of these variable taxa had such low abundances (<0.001% on average) in three of the species that they could be considered as species-specific OTUs. The concept of “sponge or host specificity” needs to be revised as molecular techniques now allow deep sequencing of associated prokaryotic assemblages. A taxon we thought unique of a particular biotope (e.g., sponges, a particular host), we are able now to detect it in other habitats but at much lower numbers. Thus, the term “sponge- or host-enriched” has been introduced (Moitinho-Silva et al., 2014). In this regard, D. antarctica and M. acerata presented an enrichment of one of the core OTUs and S. antarcticus was enriched with other two core taxa representing a 45-fold change compared to seawater. The additional three sponge core OTUs were enriched in planktonic communities (Supplementary Table S2). The same enrichment pattern was found in 30% of the variable OTUs (Supplementary Table S2), which were categorized as rare (0.06% on average) in the surrounding seawater representing a 50-fold increase.

The fact that predominant OTUs were detected in all sponges from different geographic locations and also in the surrounding seawater suggests a global distribution of these microbes through the Antarctic environment. These OTUs might be adapted to their particular niches representing different ecotypes of the same microorganisms, so it is possible that they are horizontally transmitted through strong selective mechanisms (Schmitt et al., 2012; Turon et al., 2018). The majority of the sponge prokaryotic community (99%) belonged to the “rare biosphere” with mean relative abundances below 0.1%, and over 50% could be described as host-specific taxa. The specificity of host-microbe associations seems to extend beyond some dominant taxa from the variable community into the rare biosphere (Reveillaud et al., 2014).

Beyond the existence of certain plasticity of microbiomes within a particular host, the dominance of a gammaproteobacterial OTU (Betaproteobacteriales) in M. acerata was also recently reported from another location of the Antarctic Peninsula (Cárdenas et al., 2019). This spatial stability was also documented in the present study for dominant OTUs in D. antarctica and S. antarcticus across a geographic scale <700 km. These results are in agreement with a recent study in the Caribbean that found little variations in the sponge microbiome of Cliona delitrix over small geographic scales (<300 km), while a considerable geographic distance impact over a large regional scale (>1,000 km) was reported (Easson et al., 2020). To date, published data support the combination of host identity, geography, and environment as the main forces determining the structure of sponge microbiomes (Webster et al., 2010; Schmitt et al., 2012; Easson et al., 2020; Sacristán-Soriano et al., 2020b). In Antarctica in particular, the prokaryotic signature of sponges might be also related to a biogeographic partitioning of Southern Ocean microorganisms caused by the Polar Front, such as the deficit of Cyanobacteria in Antarctica (Wilkins et al., 2013), as evidenced in the present and previous studies (Rodríguez-Marconi et al., 2015; Cárdenas et al., 2018; Moreno-Pino et al., 2020). Since most collection efforts of sponges have so far explored tropical and temperate environments, this study contributes to expand our knowledge on sponge microbiome structure in polar waters.

The fact that nearly 90% of the Antarctic sponge microbiomes had a blast hit with a sequence similarity over 97% against the sponge EMP collection would represent the ubiquity of the sponge-associated microbes through species and habitats.