On the way to specificity - Microbiome reflects sponge genetic cluster primarily in highly structured populations

4412 | Molecular Ecology. 2020;29:4412–4427.wileyonlinelibrary.com/journal/mec 1 | INTRODUCTION Sponges (phylum Porifera) are early-divergent metazoans that represent a major part of the marine benthic fauna across the world's oceans, providing fundamental services to the ecosystems where they live (Bell, 2008; Taylor et al., 2007). Despite their relatively simple body plan, sponges are known for hosting complex, dense, diverse, and highly specific microbial communities (Thomas et al., 2016). In some cases, these microbial associates comprise as much as 90% of the sponge volume and can contribute significantly to host metabolism and biochemical repertoire (Taylor et al., 2007; Webster & Thomas, 2016). Each sponge species harbours a specific symbiotic community, resulting from the combination of two proposed acquisition mechanisms: environmental (horizontal) Received: 3 March 2020 | Revised: 21 August 2020 | Accepted: 28 August 2020 DOI: 10.1111/mec.15635 ORIGINAL ARTICLE On the way to specificity - Microbiome reflects sponge genetic cluster primarily in highly structured populations Cristina Díez-Vives1 | Sergi Taboada2,3 | Carlos Leiva1,4 | Kathrin Busch5 | Ute Hentschel5 | Ana Riesgo1,6 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Molecular Ecology published by John Wiley & Sons Ltd 1Department of Life Sciences, The Natural History Museum, London, UK 2Departamento de Ciencias de la Vida, EUUS Marine Biodiversity Group, Universidad de Alcalá, Alcalá de Henares, Spain 3Departamento de Biología (Zoología), Universidad Autónoma de Madrid, Facultad de Ciencias, Madrid, Spain 4Department of Genetics, Microbiology and Statistics, Faculty of Biology, University of Barcelona, Barcelona, Spain 5GEOMAR Helmholtz Centre for Ocean Research Kiel, Research Unit Marine Symbioses, Kiel, Germany 6Department of Biodiversity and Evolutionary Biology, Museo Nacional de Ciencias Naturales de Madrid (CSIC), Madrid, Spain Correspondence Cristina Díez-Vives, Department of Life Sciences, The Natural History Museum of London, London, UK. Email: cristinadiez[email protected] Funding information European Union’s Horizon 2020 Research and Innovation Actions, Grant/Award Number: 679849; Marie SkłodowskaCurie Individual Fellowships, Grant/Award Number: 796011 Abstract Most animals, including sponges (Porifera), have species-specific microbiomes. Which genetic or environmental factors play major roles structuring the microbial community at the intraspecific level in sponges is, however, largely unknown. In this study, we tested whether geographic location or genetic structure of conspecific sponges influences their microbial assembly. For that, we used three sponge species with different rates of gene flow, and collected samples along their entire distribution range (two from the Mediterranean and one from the Southern Ocean) yielding a total of 393 samples. These three sponge species have been previously analysed by microsatellites or single nucleotide polymorphisms, and here we investigate their microbiomes by amplicon sequencing of the microbial 16S rRNA gene. The sponge Petrosia ficiformis, with highly isolated populations (low gene flow), showed a stronger influence of the host genetic distance on the microbial composition than the spatial distance. Host-specificity was therefore detected at the genotypic level, with individuals belonging to the same host genetic cluster harbouring more similar microbiomes than distant ones. On the contrary, the microbiome of Ircinia fasciculata and Dendrilla antarctica - both with weak population structure (high gene flow) - seemed influenced by location rather than by host genetic distance. Our results suggest that in sponge species with high population structure, the host genetic cluster influence the microbial community more than the geographic location. KEYWORDS host-specificity, microbiome, sponges | 4413 DÍEZ-VIVES Et al. acquisition of microbes, and parental (vertical) transmission (Björk et al., 2019; Sipkema et al., 2015; Taylor et al., 2007). This host species-specific nature of sponge-associated microbial communities is now well-established, with most studies reporting the “host identity” (i.e., host species) to be the single strongest influence on the composition of the associated microbial community (O’Brien et al., 2019; Reveillaud et al., 2014; Thomas et al., 2016). The sponge phylogeny (“host relatedness”), however, is not unequivocally linked with the microbial similarity. Sponge relatedness has been associated with microbial diversity and community composition in some studies (Easson & Thacker, 2014; Schöttner et al., 2013; Souza et al., 2017; Thomas et al., 2016). Contrastingly, other studies reported closely related sponges harbouring different microbial communities (Easson & Thacker, 2014; Schmitt et al., 2012). Host-specificity patterns below species level are even less understood, because in almost all cases, research has focused on describing the intraspecific variability of sponges with respect to environmental drivers of variation, and among these studies, results have been inconsistent. In some studies, intraspecific microbial communities differed over spatial and temporal scales, location, nutrient concentration or habitat (Anderson et al., 2010; Burgsdorf et al., 2014; Fiore et al., 2013; Luter et al., 2015; Turque et al., 2010; Weigel & Erwin, 2016). In contrast, other studies reported stable microbial communities over spatial scales, different seasons and depths (Björk et al., 2013; Erwin et al., 2012, 2015; Hentschel et al., 2002; Pita et al., 2013; Reveillaud et al., 2014; Simister et al., 2013; Taylor et al., 2005). Studies including more than one sponge species in different locations also found contrasting findings for the different species included, proposing different strengths in host-symbiont interactions for different sponge species (Cleary et al., 2013; Lee et al., 2009; Thomas et al., 2016), or an effect of the genetic variability of the sponge host (Taylor et al., 2005). Recently, few studies have revealed an apparent influence of intraspecific host genetics in structuring the microbial communities in single sponge species from the Caribbean (Easson et al., 2020; Griffiths et al., 2019; Marino et al., 2017). However, in the Indo-Pacific, microbial variation was predominantly related to geography as opposed to host genetic groups (Swierts et al., 2018). The role of subspecies host genetic divergence in determining the intraspecific variability of the microbial community, therefore, has yet to be largely determined for marine sponges. In other organisms, including animals and plants, the impact of host genetic variation on the microbiome composition has already been reported, including the human gut microbiome (Davenport, 2016; Kolde et al., 2018; Spor et al., 2011), mouse lines (Benson et al., 2010), or plant genotypes (Bouffaud et al., 2014; Qian et al., 2018). In sponges, host genotype variability is only known for some species with the appropriate markers for fine resolution (see review Pérez-Portela & Riesgo, 2018). In fact, classic analyses with ribosomal or mitochondrial DNA markers are not sensitive enough to detect such the levels of genetic variability. In turn, microsatellites and single nucleotide polymorphisms (SNPs) are powerful genetic markers that provide resolution of population structure and therefore genetic variability at local and global scales for sponges (Leiva et al., 2019; Pérez-Portela & Riesgo, 2018). Using these techniques to investigate population connectivity, moderate to high gene flow has been detected for sponges (Chaves-Fonnegra et al., 2015; Giles et al., 2015; Leiva et al., 2019; Riesgo et al., 2016; Taboada et al., 2018), and very rarely, true genetic isolation (low gene flow) of different populations has been reported (Riesgo et al., 2019). Whether this variation in gene flow and genetic divergence has any impact on the composition of the microbiome has never been tested appropriately. The goal of the present study was to determine the effect of both intraspecific genetic variation and geographic location (over 3,000 km in the Mediterranean Sea, and 740 km in the Southern Ocean) on the microbial community structure and composition of three marine sponge species. We characterized the microbiome compositions of two Mediterranean demosponges (Ircinia fasciculata and Petrosia ficiformis) and one Antarctic demosponge (Dendrilla antarctica), which present different levels of gene flow and a priori microbiome acquisition strategies. The connectivity of the sponge hosts had previously been assessed with microsatellite markers for the Mediterranean sponges (Riesgo et al., 2016, 2019) and with SNPs for the Antarctic sponge species (Leiva et al., 2019). We tested whether a possible intraspecific host-specificity signal is present both across geographic space and host genetic clusters. Our study contributes to refine our understanding of the relationship between host speciation and microbial community composition in sponges, one of the oldest animal phyla on the planet. 2 | MATERIALS AND METHODS 2.1 | Sponge sampling Three different sponge species were used in this study to examine host sponge genetic distance and associated microbial community dissimilarities. Sponge species were collected at different times and locations, and therefore, were not directly compared with each other (Table S1). A total of 168 individuals of Petrosia ficiformis, a high microbial abundance (HMA) sponge that displays exclusively horizontal transmission of symbionts (Lepore et al., 1995; Maldonado & Riesgo, 2009), were collected in shallow waters of the Mediterranean Sea along 17 locations during three sampling campaigns in July–August of different years (see details of collection in Riesgo et al., 2019). For the second sponge, Ircinia fasciculata, also a HMA sponge that in turn displays vertical transmission (Björk et al., 2019), 166 individuals were collected also in shallow waters of the Mediterranean Sea, along 11 locations in July–August of four different years (see Riesgo et al., 2016 for details). Finally, 62 individuals of Dendrilla antarctica, low microbial abundance (LMA) sponge with horizontal transmission (Koutsouveli et al., 2018), were collected from Antarctic shallow waters in seven locations along the Antarctic Peninsula in one single campaign during the 2015–2016 Austral summer (see Leiva et al., 2019). All sponge species were preserved in absolute ethanol that was replaced with fresh ethanol at least three times within 48 hr and stored at −20°C until further processed. 4414 | DÍEZ-VIVES Et al. DNA was extracted with the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) following the manufacturer's instructions with a minor modification concerning overall cell lysis time (that is, incubation was conducted overnight) and the final DNA elution step (performed twice using 50 μl of buffer EB each time). 2.2 | Genetic distances of the host Genetic clusters were assigned to individuals based on microsatellite data sets for P. ficiformis and I. faciculata (Riesgo et al., 2016, 2019), and from single-nucleotide polymorphisms (SNPs) for D. antarctica (Leiva et al., 2019) using a Bayesian clustering approach in STRUCTURE 2.3.4 (Pritchard et al., 2000), that calculates population allele frequencies and then assigns individuals to populations probabilistically (see references Leiva et al., 2019; Riesgo et al., 2016, 2019 for details of the analyses). Then, Euclidean genetic distances among individual sponge samples were calculated with GENODIVE version 2.0b23 (Meirmans & Van Tienderen, 2004) for the microsatellite data sets, and using the dist function in R v.2.14 (R Core Team, 2019) for the SNP data set. Finally, population differentiation between pairwise sampling sites and genetic clusters was also estimated with GENODIVE using the FST statistic and an infinite allele model (IAM). Significance of FST values was analysed with 20,000 permutations. 2.3 | 16S rRNA amplicon sequencing For P. ficiformis and I. fasciculata, we targeted the V3–V4 hypervariable regions of the 16S rRNA gene, while the V4 hypervariable region was used for D. antarctica. The V3V4 region was amplified using a one-step PCR with the following conditions: 98°C for 30 s, followed by 30 cycles of 98°C for 9 s, 55°C for 1 min, 72°C for 1.5 min, and a final elongation at 72°C for 10 min. We used the primer pair 341F (Muyzer et al., 1993) and 806R (Caporaso et al., 2011) in a dual-barcoding approach (Kozich et al., 2013). Verification of PCR-products was accomplished by electrophoresis on an agarose gel. Normalisation and cleaning was done with the SequalPrep Normalization Plate Kit (Invitrogen). Afterwards products were pooled equimolarly and sequenced on a MiSeq platform using v3 chemistry (2 × 300 bp) at the University Kiel, Germany (https:// www.ikmb.uni-kiel.de/). The V4 region of the 16S rRNA gene was amplified using general bacterial primers 515F-Y (Parada et al., 2016) and 806R (Apprill et al., 2015), with the Illumina adapter overhang sequences in both primers. These primers contain degenerated bases to remove the previous bias against Crenarchaeota/Thaumarchaeota, and the Alphaproteobacterial clade SAR11. We used the PCRBIO HiFi Polymerase (PCR Biosystems Ltd) under the following conditions: 95°C for 3 min, followed by 25 cycles of 95°C for 20 s, 60°C for 20 s and 72°C for 30 s, after which a final elongation step at 72°C for 5 min was performed. DNA amplification was done in duplicates, and PCR products were checked in 1% agarose gel to determine the success of amplification and the relative intensity of bands. PCR products were purified with AgencourtAMPure XP Beads (Beckman Coulter Inc.), and libraries prepared with the Nextera XT DNA Library Preparation Kit (Illumina Inc.). An equimolar pool of DNA was generated by normalizing all samples at 4 nM for sequencing. Next generation, paired-end sequencing was performed at the Natural History Museum of London (https://www.nhm.ac.uk/) on an Illumina MiSeq device using v3 chemistry (2 × 300 bp). 2.4 | Read processing, taxonomic assignment and core ASVs Raw paired reads were imported into Mothur (v.1.41.3), and an adaptation of the MiSeq SOP protocol was followed (Kozich et al., 2013). Briefly, primer sequences were removed and sequence contigs built from overlapping paired reads. The merged amplicon sequence lengths were ca. 458 bp and 298 bp for the V3–V4 and V4 regions, respectively. Sequences with >0 N bases or with >15 homopolymers were discarded. Unique sequences were aligned against the Silva reference data set (release 132), and poorly aligned sequences removed. Unoise3 (Callahan et al., 2016), which is implemented within mothur, was used for denoising (i.e., error correction) of unique aligned sequences, to infer amplicon sequence variants (ASVs), allowing one mismatch per 100 bp (Oksanen et al., 2018). Any singletons remaining at this stage were removed. Reference based chimera checking was conducted using UCHIME with the Silva reference data set and parameter minh = 0.3. ASVs were classified using the Silva database v.132, with a cutoff value of 80. ASVs classified as eukaryotic-choroplast-mitochondria or unknown were discarded, this represented less than 0.002% of sequences for the Mediterranean samples and 0.4% in the Antarctica data set. Sequences that remained unclassified further than to Kingdom bacteria or archaea, accounted for 3.1% of Mediterranean samples, and 13% Antarctic samples. Any sample with less than 1,000 sequences was discarded. Community sampling efficiency was examined using rarefaction curves. Description of the microbial community was done using the total number of ASVs transformed to relative abundances within each individual. Furthermore, for alpha diversity analyses, samples were rarefied to 5,000 sequences for P. ficiformis and I. fasciculata (discarding 19 and 29 extra samples that did not reach the minimum, respectively), and a minimum of 11,000 sequences for D. antarctica (two samples were excluded). The core microbiome was determined on the rarefied data sets using two definitions: ASVs that were present in 100% of the samples at any abundance, and in 80% of samples at any abundance. 2.5 | Statistical design and analysis To analyse the influence of genetic and geographic effect on the microbial community, both continuous and discrete variables were | 4415 DÍEZ-VIVES Et al. included for statistical analyses. Continuous genetic distances among hosts were calculated as Euclidean genetic distances. Samples with 0 Euclidean distance (i.e., clones) were discarded from further analyses. Continuous geographic distances were calculated as kilometres between sampling sites based on the GPS coordinates of each site. Discrete genetic groups were based on the host genetic clusters defined previously (see Leiva et al., 2019; Riesgo et al., 2016, 2019), and discrete geographic groups were designated at sampling collection level (locations). Correlation between community composition (Bray–Curtis dissimilarity) and both continuous distances (i.e. host Euclidean distance and geographic distances) were tested using Mantel and partial Mantel tests implemented in the R package vegan. The Mantel approach computes Pearson correlation between continuous distances, with significance based on 999 permutations of the distance matrix. Correlations were analysed globally (including all samples), using partial Mantel test controlling for the effect of a third variable, and also on each location separately and each genetic cluster. Measures of ASV richness, Shannon index, and inverse Simpson's index were calculated using the rarefied samples in R v.3.6.1. These metrics were compared among genetic clusters using analyses of variance (ANOVA). Pairwise comparisons were conducted using TukeyHSD. Beta diversity was calculated using the Bray-Curtis dissimilarity coefficient. ASVs were filtered by a relative abundance >0.01% in at least 5% of samples, leaving 1,557, 2,032 and 1,059 ASVs for P. ficiformis, I. fasciculata and D. antarctica respectively. The relative abundances were then log2 transformed prior to calculation of Bray-Curtis dissimilarities. These dissimilarity matrices were visualised using Principal coordinates analysis (PCoA) using “cmdscale” in vegan v. 2.5–6 (De Cáceres & Legendre, 2009). We compared distances among genetic clusters and locations (discrete factors) by permutational multivariate analyses of variance (PERMANOVA) using “adonis” in vegan and a type II of sums of squares for partitioning terms in unbalanced designs. When samples were collected in different years, this factor was added into the analysis. For P. ficiformis, beta diversity analyses were re-run after transforming the data to presence absence (and using Jaccard distances), and at genus level (ASVs abundances belonging to the same genera were aggregated), to see the effect on alternative data sets. Beta diversity analyses were initially performed on the entire data sets, however, since genetic, geographic and year effect may be confounding, we repeated analyses for each independent location, genetic cluster and year of sampling. In this case, locations with insufficient host genetic variation were excluded, meaning that only locations with more than one identified genetic cluster and with at least three replicates (for P. ficiformis and I. fasciculata) were considered. In D. antarctica we reduced the number to two replicates due to the lower number of samples in the data set. This resulted in the informative locations BLA, LIG and NAP for P. ficiformis; CALA, CRO, NAP, TOSS, CAB and ESC for I. fasciculata; and CIE, KG, ADE and HM for D. antarctica. For the analysis of each host genetic cluster separately, we included only genetic clusters present in more than one location with 2 or 3 replicates, similarly as before. These were Pf3, Pf4, Pf5, Pf6, Pf8 for P. ficiformis; If1, If2, If4, If5 for I. fasciculata; and Da1, Da2, Da3, Da4 for D. antarctica. Moreover, to disentangle genetic and geographic effects, these selected locations and genetic clusters were also analysed together. 2.6 | Indicator species in Petrosia ficiformis Indicator species analysis (identification of species associated with or indicative of groups of samples) was conducted using the R package indicspecies (Dufrene & Legendre, 1997) to identify microbial taxa characteristic of each genetic cluster in P. ficiformis. This analysis assesses the strength of the relationship between ASVs abundance and different host genetic clusters by comparing ASVs abundance in microbiotas of one genetic cluster to their abundance in the others. Enrichment values were calculated for each indicator as a log-transformed ratio of each two groups. The Indicator Value index (Dufrene & Legendre, 1997) is the product of two components: (a) the “the specificity or positive predictive value” as the probability that the surveyed ASVs only belongs to the target genetic cluster; and (b) the probability of finding the ASVs in all samples belonging to the genetic cluster, called “the fidelity or sensitivity” component. The statistical significance of this relationship is tested using a permutation test. Multiple pairwise comparisons were corrected based on the Benjamini–Yekutieli false discovery rate control using “p.adjust” function of the stats library package in R. 3 | RESULTS 3.1 | Genetic distances and genetic clusters of sponge hosts Petrosia ficiformis and Ircinia fasciculata were collected across the Mediterranean Sea, including 168 and 163 samples, respectively. Dendrilla antarctica included 62 samples from the Southern Ocean (Figure 1, Table S1). The data set of P. ficiformis was grouped into seven genetic clusters (Pf2–Pf6 and Pf8 as reported in Riesgo et al., 2019), five genetic clusters (If1–If5) for I. fasciculata (Riesgo et al., 2016), and five genetic clusters (Da1–Da4 plus an unclustered group, i.e., individuals with multiple genetic clusters assigned in which none was dominant over the others) for D. antarctica (Leiva et al., 2019). Euclidean genetic distances between individuals for each sponge species can be found in Tables S2–S4. We identified one clone in P. ficiformis and three clones in I. fasciculata which were discarded from the following analyses, as well as the unclassified group of D. antarctica. Fixation indices (FST) between genetic clusters were the largest in P. ficiformis ranging from 0.048 to 0.276, followed by D. antarctica (from 0 to 0.156) and I. fasciculata (from 0 to 0.078; Table S5). These FST values can be considered as low, moderate and high gene flow, respectively (Figure 2), when compared to other sponges (Pérez-Portela & Riesgo, 2018). 4416 | DÍEZ-VIVES Et al. In P. ficiformis, the PCA on the Euclidean distances showed that the first and second principle coordinates explained 24% of the total variation among the samples (Figure 3a). This ordination revealed a clearer clustering of sponge individuals by genetic group rather than by location (i.e., less overlapping of groups). In addition, and in concordance with having the lowest FST values for P. ficiformis (i.e., 0.048), genetic groups Pf4 and Pf8 were found to be more closely related than with the rest, presenting overlapping clustering. The PCA in D. antarctica explained 26% of the total variation (Figure S1a), showing separation by genetic cluster but not by location, which suggests a stronger host genetic effect over spatial effects. In the case of I. fasciculata, the PCA explained 11.8% of the total variation among the samples (Figure S1b). This ordination did not show clear clustering neither by genetic cluster nor by location. 3.2 | Microbial composition A total of 107,913, 66,524, and 18,165 unique ASVs were found among all samples of P. ficiformis, I. fasciculata, and D. antarctica, respectively. Abundance values ranged from 1,057 to 115,327 sequences per sample (Table S6). Rarefaction curves showed good representation of amplicon sequences present in D. antarctica, but a number of samples did not approach asymptotes for the two other species, suggesting more species would be observed with greater sequencing effort (Figure S2a). In P. ficiformis, the observed ASVs were assigned to 35 phyla and 354 genera. The most abundant phyla were the Chloroflexi with 28.1% mra (mean relative abundance), Proteobacteria (25.4% mra) and Acidobacteria (10.6% mra; Figure S2b). Dominant classes FIGURE 1 Distribution map of sampling sites over the Mediterranean Sea for P. ficiformis and I. fasciculata and along the Antarctic Peninsula for D. antarctica. The colours of the genetic clusters follow those seen in the original papers (see main text). Location names correspond to: CAR, Carboneras; CART, Cartagena; BLA, Blanes; FEL, Sant Feliu; ULL, Ullastres; MRS, Marseille; LIG, Liguria; NAP, Naples; SLO, Slovenia; SCRO, South Croatia; JECRO, Jelsa Croatia; CRE, Creta; ISR, Israel; TAR, Tarifa; ALI, Alicante; CALA, Calafat; CAB, Cabrera; ESC, Escala; TOSS, Tossa; CAIA, Caials; COR, Corsica; and CRO, Croatia in the Mediterranean Sea; and ADE, Adelaide Island; PAR, Paradise Bay; CIE, Cierva Cove; DEC, Deception Island; HM, Half Moon Island; KG, King George Island; and OH, O'Higgins Bay in the Antarctic Peninsula [Colour figure can be viewed at wileyonlinelibrary.com] West Mediterranean East Mediterranean 70ºW 65ºW 60ºW 55ºW 65ºS Larsen Ice Shelf Weddell Se a ADE CIE DEC HM KG OH PAR CALA ESC ALI BLA TOSS CAIA COR CAB NAP TAR CRO If1 If2 If3 If4 If5 Genetic Cluster Da1 Da2 Da3 Da4 Unc Genetic Cluster D. antarctica I. fasciculata ISR CRE JECRO SCRO SLO NAP LIG MRS FEL ULL BLA CART CAR Pf2 Pf3 Pf4 Pf5 Pf6 Genetic Cluster Pf8 P. ficiformis | 4417 DÍEZ-VIVES Et al. and orders can be found in Table S7. Interestingly, the two most abundant individual ASVs belonged to the phyla Nitrospira and Dadabacteria (2.9% and 2.5% mra, respectively), which were not among the dominant phyla. Looking at the core microbiota, only two ASVs were shared across all samples, but a total of 55 ASVs were shared in 80% of all individuals, and these represented from FIGURE 2 Fixation index values (FST) for 17 sponge species taken from the literature. The three sponge species studied here are highlighted in bold. The legend shows a colour coding indicating the geographical span of the sampling used for each species [Colour figure can be viewed at wileyonlinelibrary.com] 0.0 0.1 0.2 0.3 0.4 Sponge Species FST Scopalina lo phyropoda Paraleucilla magna Plenaster craigi Stylissa carteri Ircinia fasciculata Xestospongia sp. Xestospo ngia testudinaria Paraleucilla magna Dendrilla antarctica Spongia officinalis Clathrina aurea Scopalina lophyropoda Aphrocallistes vastus Cliona delitrix Spongia lamella Crambe crambe Petrosia ficiformis sampled over: 1000 km 2000 km 100-500 km 10-100 km FIGURE 3 Ordination plots for Petrosia ficiformis showing (a) clustering of host genetic distances of all samples coloured by the corresponding location (left side) and by the assigned genetic cluster (right side), and (b) clustering of microbiome dissimilarities of all samples, coloured as in (a). Centroids are marked with their respective factor label and groups are circled with 0.7 data coverage for data ellipses [Colour figure can be viewed at wileyonlinelibrary.com] BLA CAR CART CRE FEL ISR JECRO LIG MRS NAP SCRO SLO ULL PCoA 1 [14.5%] PCoA 2 [9.6%] −1.0 −0.5 0.0 0.51.0 1.5 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 Pf2 Pf3 Pf4 Pf5 Pf6 Pf8 PCoA 1 [14.5%] −1.0−0.50.0 0.51.0 1.5 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 BLA CAR CART CRE FEL ISR JECRO LIG MRS NAP SCRO SLO ULL PCoA 1 [17.5%] PCoA 2 [9.6%] −0.3−0.2 −0.1 0.0 0.1 0.2 0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 Pf2 Pf3 Pf4 Pf5 Pf6 Pf8 PCoA 1 [17.5%] −0.3 −0.2−0.10.0 0.1 0. 20 .3 −0.2 −0.1 0.00.1 0.20.3 (b) Microbiome composition ordination (a) Host genetic distance ordination By location By genetic cluster By location By genetic cluster 4418 | DÍEZ-VIVES Et al. 18.9% to 47% of the total abundance of the microbiome (75% of samples had at least 30.5% mra). The core included nine of the most abundant phyla (Table S7), with Acidobacteria, Chloroflexi and Proteobacteria including 67%, 34% and 28% of its original abundance among the core ASVs respectively. We also checked the relative abundances of ASVs that were unique to specific collection sites (site-specific ASVs), and these represented relatively low percentages, i.e., 3.04 ± 1.16% mra, with maximum values found in ISR (Table S8). For I. fasciculata, the observed ASVs were assigned to 34 phyla, and 613 genera. The most abundant phyla were the Proteobacteria (36.3% mra), Cyanobacteria (19.4% mra) and Chloroflexi (18.7% mra; Figure S2b, Table S9). The two most abundant ASVs belonged to Cyanobacteria and Dadabacteria (16.1 and 4.4% mra, respectively). Regarding the core microbiota, one ASV was present in all samples, and 37 ASVs were present in 80% of all samples, and these represented from 7.8% to 71.8% of the total abundance of the microbiome (75% of samples had at least 41.8% mra). The core included 11 of the most abundant phyla (Table S9). The phylum Cyanobacteria represented 19.4% mra, and its most abundant ASV (17.4% mra) was among the core bacteria. The phylum Chloroflexi, however, with a similar mean relative abundance (18.7%) was not so faithfully shared among all individuals, with a core of 4 ASVs including only 3.2% mra. Proteobacterial ASVs were conserved around a 44.5% of their original abundances in the core. Site-specific ASVs represented 4 ± 3.78% mra (Table S8), which included higher values and larger variability than P. ficiformis samples. ASVs in D. antarctica were assigned to 50 phyla, and 913 genera. The most abundant phyla were the Proteobacteria (64.5% mra), Bacteroidetes (16.0% mra) and Verrucomicrobia (3.2% mra; Figure S2b). The remaining phyla had less than 1.5% mra (Table S10). The core microbiota was constituted by four ASVs present in all 62 samples, and 41 ASVs in more than 80% of the samples. The percentage these represented varied from 19.2% to 74.7% of the microbial community abundance (75% of samples had at least 50.7%). The most abundant core ASVs belonged to Proteobacteria representing 35.5% mra, which was 55.2% of the total Proteobacteria abundance in the sponge. Dendrilla antarctica site-specific ASVs represented the lowest values of the three species with 0.79 ± 0.78% mra (Table S8). Because the sponges P. ficiformis and I. fasciculata were amplified with different primers than D. antarctica, a direct comparison between the ASVs data sets was not possible. Therefore, we used the taxonomic annotation (at the genus level) to look at general differences between the species knowing the limitations that this approach represented. P. ficiformis and I. fasciculata shared most of the genera, while D. antarctica hosted different ones (Figure S2c), which was also reflected in the distance among the samples in a PCoA plot (Figure S2d). However, the specific microbial ASVs within these genera were different among P. ficiformis and I. fasciculata (Figure S2e), which is in agreement with their microbiomes being species-specific. 3.3 | Correlation of host genetic distance and spatial distance with microbial dissimilarity When comparing full data sets of host genetic distances and microbial dissimilarity, Pearson's r values from the partial Mantel tests ranged from 0.123 to 0.252 (Table S11), indicating moderate correlation signals, and large variability of the samples (Figure 4). Previous ecological simulations indicated that with moderate-strong signals, Pearson's r were ca. 0.2 on average (Mazel et al., 2018). In turn, Pearson's r values for the comparison between geographical distances and microbial Bray-Curtis dissimilarity values ranged from 0.155 to 0.407 (Figure 4, Table S11). Comparing the effects of both host genetic distance and spatial distance on the microbial community, P. ficiformis showed a stronger correlation with the host genetic distance, while in D. antarctica and I. fasciculata the correlation was stronger with the geographical distances between locations (Figure 4, Table S11). To remove the bias of the geographical location and year, we analysed the effect of phylosymbiosis in each site independently and each sampling year. In P. ficiformis, Pearson's r values ranged from −0.33 to 0.64. Four out of the 14 locations for P. ficiformis showed strong correlations (rm(MH) > .30, with p < .01) between the microbial community Bray-Curtis dissimilarity and host genetic distance, these were BLA_2013, LIG_2008 and NAP_2013, and marginally BLA_2006 (p = .03; Figure S3a and Table S11). The remaining 10 locations, however, included a dominant single genetic cluster (meaning that if more than one genetic clusters were present, the other ones had less than three replicates), and, in those cases, the genetic variation was probably insufficient to detect any correlation. In D. antarctica Pearson's r values ranged from −0.38 to 0.99. One out of the seven locations for D. antarctica (ADE) had high significant correlation (rm(MH) > .40, p < .01) between the microbial dissimilarity and the host genetic distance. Three other locations, CIE, KG and HM, presented multiple genetic clusters but did not show significant positive correlation (Figure S3a and Table S11). For I. fasciculata, Pearson's r values ranged from −0.46 to 0.20, and two out of fifteen locations (i.e. CALA and NAP) had high correlations (rm(MH) > .50, p = .001; Figure S3a and Table S11). Similarly as before, another two locations, CRO_2013 and TOSS_2010, presented multiple genetic clusters but no correlation was detected. The rest of locations included a dominant single genetic cluster, and therefore low genetic variance. Similar results were observed when considering the host genetic clusters independently and testing the effect of geographical distance over microbiome dissimilarity (Figure S3b and Table S11). All the genetic clusters including multiple locations had medium to strong correlation for P. ficiformis, while for D. antarctica only Da1 showed correlation (rm(MD) = .30, p = .002), and for I. fasciculata two of the three genetic clusters including multiple locations showed high correlation (r > .40, p < .05; Figure S3b and Table S11). | 4419 DÍEZ-VIVES Et al. 3.4 | Analysis of variance of microbial communities by genetic cluster and locations 3.4.1 | Alpha diversity The Shannon index was used to unravel whether the genetic cluster or the locations exhibited different diversity patterns. The alpha diversity ranged between 4.42–6.22, 2.14–5.66, and 2.16–4.98 for P. ficiformis, I. fasciculata and D. antarctica, respectively (Table S6 and Figure S4). Diversity was significantly different among genetic clusters of P. ficiformis (ANOVA, F12,136 = 9.73, p < .001); and I. fasciculata (ANOVA, F10,123 = 3.28, p < .001; Table S12), due to pairwise differences between Pf4 and Pf5 in the former, and between If4–If3 and If4–If5 in the later. Furthermore, the alpha diversity across locations was also significantly different for these two species (ANOVA, F6,141 = 4.03, and F4,129 = 8.88, p < .001, Table S12), which was associated with JECRO having lower diversity values compared to all other locations in P. ficiformis, and between CALA and NAP in I. fasciculata. Dendrilla antarctica showed no significant differences associated with either genetic cluster or location (ANOVA, p > .01). 3.4.2 | Beta diversity An ordination of the microbial community composition by location and genetic cluster is shown in Figure 3b and Figure S1c–d. For the Mediterranean sponges, in general, samples from the Eastern Mediterranean waters harboured different communities than the Western areas. Examples are ISR, CRE and JECRO for P. ficiformis that appeared separated from the rest of the samples (Figure 3b), and CRO and NAP for I. fasciculata that clustered away from the rest of samples (Figure S1d). In D. antarctica, samples from KG, ADE and DEC harboured more different communities than the other four locations (Figure S1c). The microbiomes of the P. ficiformis sponges showed separation of samples by genetic cluster (Figure 3b), with a noteworthy overlap between Pf4 and Pf8 samples. This pattern was also observed in the ordination plot of the host Euclidean distances FIGURE 4 Correlation plots of microbiome dissimilarity (Bray Curtis) versus host genetic distance (Euclidean distance) corrected by the spatial distance on the top row (MH|D); and microbiome dissimilarity (Bray Curtis) versus spatial distance (kilometres between locations) corrected by host genetic distances on the bottom row (MD|H), for the three sponge species. Mantel test statistic (Pearson's r) and p-values (p) are shown inside each plot [Colour figure can be viewed at wileyonlinelibrary.com] 1.01.5 2.02.5 3.03.5 4.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Petrosia ficiformis Euclidean distance (host) Bray Curtis dissimilarity (microbiome) 1.52.0 2.53.0 3.54.0 0.4 0.6 0.8 1.0 Ircinia fasciculata 10 15 20 25 30 0.2 0.4 0.6 0.8 Dendrilla antarctica Euclidean distance (host) Euclidean distance (host) 0 500 1000 1500 2000 2500 3000 0.20.3 0.40.5 0.6 0.70.8 0.9 Spatial distance (location) Bray Curtis dissimilarity (microbiome) 0 500 1000 1500 2000 2500 3000 0.40.6 0.81.0 Spatial distance (location) 0 200 400 600 0.20.4 0.6 0.8 Spatial distance (location) rm(MH|D) = 0.252, p = 0.001 rm(MH|D) = 0.123, p = 0.018 rm(MH|D) = 0.173, p = 0.003 rm(MD|H) = 0.155, p = 0.002 rm(MD|H) = 0.407, p = 0.001 rm(MD|H) = 0.235, p = 0.001 4420 | DÍEZ-VIVES Et al. by genetic cluster (Figure 3a), which is here mirrored by their microbial community. Ordination coloured by location showed a large overlap of samples. For D. antarctica and I. fasciculata in general groups were not clearly separated neither by location nor genetic cluster (Figure S1c,d). Differences in the number of observed genetic clusters in each location, and possible confounding factors occurring among the variables, prevented the statistical analyses using all samples. Therefore, in order to perform a factorial analysis of the effect of the host genetic cluster and the spatial distance on the microbial composition, we subset the data set to contain only samples that allowed meaningful comparisons (see methods). This means that for P. ficiformis we used only independent locations containing more than one genetic cluster with at least three replicates (namely BLA, LIG, and NAP), and there was a clear grouping of samples associated with their host genetic cluster (PERMANOVA test confirmed significant clustering; Figure S5a and Table S13a). Moreover, within individual genetic clusters, samples diverged according to their sampling location (Figure S5b and Table S13a). To unravel which factor contributed more to the variance, we constructed a data set including BLA, LIG and NAP, and the genetic clusters Pf3, Pf4, Pf5, and Pf8. Generally, there was a clear clustering of genetic groups of P. ficiformis regardless of location, with the expected overlap of genetic clusters Pf4 and Pf8 (Figure 5). Both, location and genetic cluster, influenced the microbial communities (Table S13b), but the host genetic cluster had a stronger effect, explaining 19.9% of the variability on the microbiome composition, while location explained 3.4% and year 3.9%, and the interaction of these factors 4.8%. In pairwise comparisons, all genetic clusters were different to each other (p = .001) except for the pair Pf4–Pf8 (p = .011; Table S13). To further test whether differences in the microbial community were driven by specific assemblages of ASVs or by differences in relative abundances of shared ASVs, we repeated the analysis using Jaccard distances on the presence/absence of ASVs (Figure S6a). We also tested the differences at genus level (by aggregating the abundance of ASVs belonging to same genus (Figure S6b). These two analyses showed that genetic cluster was still the dominant factor compared to location or year, but the effect was less strong than when accounting for relative abundances at ASVs level (i.e., 11.9% and 11.8% for the presence/absence and genus level analyses, respectively; Table S13b). Moreover, in pairwise comparisons using genus level, not only the pair Pf4–Pf8 was not significantly different, but also the genetic cluster Pf3 was not different to either Pf4 or Pf8 (Table S13c). Four locations could be tested independently for D. antarctica (CIE, KG, HM and ADE, Figure S5a). Microbial communities were not different between the different genetic clusters for any location (p > .01, Table S13a). Within individual genetic clusters, however, microbial communities from both Da1 and Da4 were significantly different among locations (Figure S5b, Table S13a). Combining CIE, KG, HM and ADE, ANOVA results showed that location had a stronger effect explaining 25.4% of variability while genetic cluster explained 13.2%, and the interaction of these factors 15.3% (Figure S1, Table S12b), but pairwise comparisons showed that only Da1 versus Da3 and Da4, and Da2 versus Da4 presented different compositions (Table S13c). In I. fasciculata, the selected data set included six locations (CALA, CRO, NAP, TOSS, CAB and ESC), and genetic clusters If1 to If4 (Figure S5a). Samples were not different by genetic clusters (p > .01, Table S13a), except for CAB and ESC (p = .001). However, in these last two locations, samples from either genetic cluster were collected in different years (Table S1) preventing the disclosure of the main factor. Among individual genetic clusters, If4 and If5 showed grouping of microbial communities by location (p = .001, Table S13a). Combining the locations (CALA, CRO, NAP, and TOSS but excluding CAB and ESC; Figure 5), PERMANOVA indicated that microbial composition was not significantly different by genetic cluster (p > .01), but it was by location (p = .001), which explained 10.7% of variability. 3.5 | Specificity of the microbiome within genetic clusters in Petrosia ficiformis Since P. ficiformis was the sponge species with the strongest specificity between microbiome and host genetic cluster, our goal was FIGURE 5 Ordination plots for selected (i.e. informative) locations and genetic clusters in the three sponge species. Samples are given symbols by the corresponding location, colours display the assigned genetic cluster, and dotted lines circle samples from different years [Colour figure can be viewed at wileyonlinelibrary.com] −0.2 −0.1 0.0 0.1 0.2 0.3 PCoA 1 [27%] PCoA 2 [9.6%] P. ficiformis Genetic cluster Pf3 Pf4 Pf5 Pf8 Location BLA LIG NAP −0.2 −0.1 0.0 0.1 0.2 0.3 PCoA 1 [29.5%] PCoA 2 [15.3%] I. fasciculata −0.2 0.0 0.2 −0.20.0 0.2 0.4 −0.4 −0.2 0.00.2 −0.3 −0.2 −0.10.0 0.10.2 0.3 PCoA 1 [22.5%] PCoA 2 [15.2%] Genetic cluster Da1 Da2 Da3 Da4 Location CIE KG HM ADE D. antarctica Location CALA CRO NAP TOSS Genetic cluster If1 If2 If3 If4 | 4427 DÍEZ-VIVES Et al. Schmitt, S., Tsai, P., Bell, J., Fromont, J., Ilan, M., Lindquist, N., Perez, T., Rodrigo, A., Schupp, P. J., Vacelet, J., Webster, N., Hentschel, U., & Taylor, M. W. (2012). Assessing the complex sponge microbiota: Core, variable and species-specific bacterial communities in marine sponges. The ISME Journal, 6(3), 565–576. https://doi.org/10.1038/ ismej.2011.116 Schöttner, S., Hoffmann, F., Cárdenas, P., Rapp, H. T., Boetius, A., & Ramette, A. (2013). Relationships between host phylogeny, host type and bacterial community diversity in cold-water coral reef sponges. PLoS One, 8(2), e55505–e55511. https://doi.org/10.1371/journ al. pone.0055505 Shaw, K. L., & Mullen, S. P. (2014). Speciation continuum. Journal of Heredity, 105(S1), 741–742. https://doi.org/10.1093/jhere d/esu060 Simister, R., Taylor, M. W., Rogers, K. M., Schupp, P. J., & Deines, P. (2013). Temporal molecular and isotopic analysis of active bacterial communities in two New Zealand sponges. FEMS Microbiology Ecology, 85(1), 195–205. https://doi.org/10.1111/1574-6941.12109 Sipkema, D., de Caralt, S., Morillo, J. A., Al-Soud, W. A., Sørensen, S. J., Smidt, H., & Uriz, M. J. (2015). Similar sponge-associated bacteria can be acquired via both vertical and horizontal transmission. Environmental Microbiology, 17(10), 3807–3821. https://doi. org/10.1111/1462-2920.12827 Souza, D. T., Genuário, D. B., Silva, F. S. P., Pansa, C. C., Kavamura, V. N., Moraes, F. C., Taketani, R. G., & Melo, I. S. (2017). Analysis of bacterial composition in marine sponges reveals the influence of host phylogeny and environment. FEMS Microbiology Ecology, 93(1), fiw204. https://doi.org/10.1093/femse c/fiw204 Spor, A., Koren, O., & Ley, R. (2011). Unravelling the effects of the environment and host genotype on the gut microbiome. Nature Reviews Microbiology, 9(4), 279–290. https://doi.org/10.1038/nrmic ro2540 Swierts, T., Cleary, D. F. R., & de Voogd, N. J. (2018). Prokaryotic communities of Indo-Pacific giant barrel sponges are more strongly influenced by geography than host phylogeny. FEMS Microbiology Ecology, 94(12), 716–812. https://doi.org/10.1093/femse c/fiy194 Taboada, S., Riesgo, A., Wiklund, H., Paterson, G. L. J., Koutsouveli, V., Santodomingo, N., Dale, A. C., Smith, C. R., Jones, D. O. B., Dahlgren, T. G., & Glover, A. G. (2018). Implications of population connectivity studies for the design of marine protected areas in the deep sea: An example of a demosponge from the Clarion-Clipperton Zone. Molecular Ecology, 27(23), 4657–4679. https://doi.org/10.1111/mec.14888 Taylor, M., Radax, R., Steger, D., & Wagner, M. (2007). Sponge-associated microorganisms: Evolution, ecology, and biotechnological potential. Microbiology and Molecular Biology Reviews, 71(2), 295–347. https:// doi.org/10.1128/MMBR.00040-06 Taylor, M. W., Schupp, P. J., de Nys, R., Kjelleberg, S., & Steinberg, P. D. (2005). Biogeography of bacteria associated with the marine sponge Cymbastela concentrica. Environmental Microbiology, 7(3), 419–433. https://doi.org/10.1111/j.1462-2920.2004.00711.x Thomas, T., Moitinho-Silva, L., Lurgi, M., Björk, J. R., Easson, C., Astudillo-García, C., Olson, J. B., Erwin, P. M., López-Legentil, S., Luter, H., Chaves-Fonnegra, A., Costa, R., Schupp, P. J., Steindler, L., Erpenbeck, D., Gilbert, J., Knight, R., Ackermann, G., Victor Lopez, J., … Webster, N. S. (2016). Diversity, structure and convergent evolution of the global sponge microbiome. Nature Communications, 7, 11870. https://doi.org/10.1038/ncomm s11870 Thompson, L. R., Sanders, J. G., McDonald, D., Amir, A., Ladau, J., Locey, K. J., Prill, R. J., Tripathi, A., Gibbons, S. M., Ackermann, G., NavasMolina, J. A., Janssen, S., Kopylova, E., Vázquez-Baeza, Y., González, A., Morton, J. T., Mirarab, S., Zech Xu, Z., Jiang, L., … Knight, R. (2017). A communal catalogue reveals Earth’s multiscale microbial diversity. Nature, 104, 11436–11524. https://doi.org/10.1038/natur e24621 Tikhonov, M., Leach, R. W., & Wingreen, N. S. (2015). Interpreting 16S metagenomic data without clustering to achieve sub-OTU resolution. The ISME Journal, 9, 68–80. https://doi.org/10.1038/ismej.2014.117 Turon, M., Cáliz, J., Garate, L., Casamayor, E. O., & Uriz, M. J. (2018). Showcasing the role of seawater in bacteria recruitment and microbiome stability in sponges. Scientific Reports, 8(1), 15201. https://doi. org/10.1038/s41598-018-33545-1 Turque, A. S., Batista, D., Silveira, C. B., Cardoso, A. M., Vieira, R. P., Moraes, F. C., Clementino, M. M., Albano, R. M., Paranhos, R., Martins, O. B., & Muricy, G. (2010). Environmental shaping of sponge associated archaeal communities. PLoS One, 5(12), e15774. https:// doi.org/10.1371/journ al.pone.0015774 Wagner, M. R., Lundberg, D. S., del Rio, T. G., Tringe, S. G., Dangl, J. L., & Mitchell-Olds, T. (2016). Host genotype and age shape the leaf and root microbiomes of a wild perennial plant. Nature Communications, 7(1), 12151. https://doi.org/10.1038/ncomm s12151 Webster, N. S., & Thomas, T. (2016). The Sponge Hologenome. MBio, 7(2), e00135–e216. https://doi.org/10.1128/mBio.00135-16 Weigel, B. L., & Erwin, P. M. (2016). Intraspecific variation in microbial symbiont communities of the sun sponge, Hymeniacidon heliophila, from intertidal and subtidal habitats. Applied and Environmental Microbiology, 82(2), 650–658. https://doi.org/10.1128/AEM.02980-15 SUPPORTING INFORMATION Additional supporting information may be found online in the Supporting Information section. How to cite this article: Díez-Vives C, Taboada S, Leiva C, Busch K, Hentschel U, Riesgo A. On the way to specificity - Microbiome reflects sponge genetic cluster primarily in highly structured populations. Mol Ecol. 2020;29:4412–4427. https:// doi.org/10.1111/mec.15635