Host population genetics and biogeography structure the microbiome of the sponge Cliona delitrix

Abstract Sponges occur across diverse marine biomes and host internal microbial communities that can provide critical ecological functions. While strong patterns of host specificity have been observed consistently in sponge microbiomes, the precise ecological relationships between hosts and their symbiotic microbial communities remain to be fully delineated. In the current study, we investigate the relative roles of host population genetics and biogeography in structuring the microbial communities hosted by the excavating sponge Cliona delitrix. A total of 53 samples, previously used to demarcate the population genetic structure of C. delitrix, were selected from two locations in the Caribbean Sea and from eight locations across the reefs of Florida and the Bahamas. Microbial community diversity and composition were measured using Illumina‐based high‐throughput sequencing of the 16S rRNA V4 region and related to host population structure and geographic distribution. Most operational taxonomic units (OTUs) specific to Cliona delitrix microbiomes were rare, while other OTUs were shared with congeneric hosts. Across a large regional scale (>1,000 km), geographic distance was associated with considerable variability of the sponge microbiome, suggesting a distance–decay relationship, but little impact over smaller spatial scales (<300 km) was observed. Host population structure had a moderate effect on the structure of these microbial communities, regardless of geographic distance. These results support the interplay between geographic, environmental, and host factors as forces determining the community structure of microbiomes associated with C. delitrix. Moreover, these data suggest that the mechanisms of host regulation can be observed at the population genetic scale, prior to the onset of speciation.

The microbial community composition and functional ecology of sponges might be influenced by local environmental factors, especially in LMA species, since sponges are continuously exposed to a diverse and dynamic consortium of seawater microorganisms via filter feeding. Despite this exposure, most sponge species host microbial assemblages that are distinct from those found in the surrounding seawater (Taylor et al., 2013). In some cases, specific associations are maintained by vertical transfer of microorganisms from parents to eggs and larvae (Diaz, Thacker, Rützler, & Piantoni Dietrich, 2007;Olson & Gao, 2013;Pita, López-Legentil, & Erwin, 2013;Reveillaud et al., 2014;Schmitt, Weisz, Lindquist, & Hentschel, 2007;Sharp, Eam, Faulkner, & Haygood, 2007).
Across the large biogeographic range of some sponge species, one might expect a highly variable nutritional environment (with varying composition and concentrations of particulate and dissolved organic matter and inorganic nutrients), which might influence host reliance on symbiont-derived nutrition. However, data from a limited number of species suggest that this reliance is likely species-specific (Freeman & Thacker, 2011). Some recent studies have demonstrated environmental impacts on sponge microbial communities, with some variation across habitats (i.e., intertidal vs. subtidal, inshore vs. offshore reefs, open water vs. marine lakes) (Cleary et al., 2013;Luter et al., 2015;Weigel & Erwin, 2016), seasons (Hardoim et al., 2012;White et al., 2012), and latitude (Anderson, Northcote, & Page, 2010; Marino, Pawlik, López-Legentil, & Erwin, 2017; Taylor et al., 2013). Taken together, these studies suggest that sponge-associated microbial communities are influenced by both host-specific and environmental factors (Easson & Thacker, 2014;Marino et al., 2017;Taylor, Radax, Steger, & Wagner, 2007;Thomas et al., 2016). At a global scale, microbiomes of individual sponge species exhibited relatively low within-host-species variability, suggesting that sponge tissues can form a generally selective habitat at the scale of individual host species; this trend was consistent irrespective of microbial diversity or abundance (Thomas et al., 2016). However, environmental influences on microbial community structure were not explicitly tested in this large-scale study.
Sponges reproduce by brooding their embryos or through broadcast spawning, in most cases, of fertilized eggs (Maldonado & Bergquist, 2002). Their larval dispersal is limited due to a variety of factors, including planktonic larval duration (<72 hr), transport of eggs and embryo development (<2 weeks), and limited swimming capabilities that leave them at mercy of ocean currents (Chaves-Fonnegra, Feldheim, Secord, & Lopez, 2015;Maldonado & Riesgo, 2008;Maldonado & Young, 1996). In the Caribbean and western Atlantic, the population structure of sponge species tends to exhibit a high degree of isolation, while connectivity varies in relation to life-history strategies and the speed of oceanographic currents (Chaves-Fonnegra et al., 2015;Debiasse, Richards, & Shivji, 2010;Richards, Bernard, Feldheim, & Shivji, 2016). Despite these limit- DeBiasse, Richards, Shivji, & Hellberg, 2016;Swierts et al., 2017), which makes the opportunity to couple this type of data with corresponding bacterial symbiont signals even more compelling. Thus, understanding the intricacies of how these ancient animals can successfully adapt to and colonize divergent environments is still an open question and one of the greater interests under current conditions of global climate change.
As sponges have dispersed and speciated across the globe, it is clear that they have also acquired new microbial symbionts that might have helped them adapt or acclimate to novel nutritional environments (Thomas et al., 2016). While strong patterns of host specificity have been consistently observed, the relationship between host speciation and microbial community composition remains equivocal for many species. Sponge speciation is a continuous process involving slow genetic divergence over time often resulting from reproductive isolation. At intermediate steps of host speciation, potential microbial community divergence might simultaneously occur due to new selective pressures in novel environments, but to date, this question has only been investigated in a few sponge species (Griffiths et al., 2019;Swierts, Cleary, & de Voogd, 2018). Considering that sponges tend to display a high degree of population structure with low dispersal and low connectivity among populations (Chaves-Fonnegra et al., 2015), comparisons of microbial communities among genetically distinct sponge host populations allow for an evaluation of the potential effects of subspecies host genetic divergence on microbial community structure. In the current study, we coupled sponge population genetics and Illumina-based microbiome sequencing to investigate the relative influences of host population genetics and biogeography on the microbial communities of the excavating sponge, Cliona delitrix ( Figure 1). We sampled four distinct populations of C. delitrix that were collected at 10 sites across the Caribbean and western Atlantic, with two of the populations having parapatric distributions across Florida and Bahamas reefs. We investigated the microbial taxa that were conserved across this sample set, described their presence and abundance in other sponge species and environments, and tested how geography and population structure are related to microbial community diversity and composition in C. delitrix.  Table S1).

| Statistical design
To analyze and relate population genetic structure and geographic distance among hosts, both discrete and continuous population and geographic variables were included for statistical analyses (Table 1).
Discrete population groups were based on population genetic clusters of the host from Chaves-Fonnegra et al. (2015) and are treated as genetic populations in the present study. Discrete geographic groups were designated at the reef level (site of collection). Continuous genetic distances among populations were calculated as the Bruvo distance, a common method for estimating intraspecific genetic distances among samples based on population genetic data (Bruvo, Michiels, D'souza, & Schulenburg, 2004, Figure 3). Continuous geographic distances were calculated as the Euclidean distance among sampling sites based on the GPS coordinates of each site (Chaves-Fonnegra et al., 2015).

| DNA extraction and sequencing
Total genomic DNA was extracted from each sponge sample using the PowerSoil DNA Extraction Kit (MoBio) following the standard protocols of the Earth Microbiome Project (EMP; www.earth micro biome.org; Thompson et al, 2017). Extracted DNA was shipped to the University of Colorado, Bolder, CO, USA, where the V4 region of the 16S rRNA gene was amplified using the primer set 515F-806rB and then sequenced on an Illumina HiSeq 2,500 platform (Illumina) following the EMP standard protocols. Sequence processing was performed following the bioinformatics methods outlined in Thomas et al. (2016) and using the program Mothur (Schloss et al., 2009). Sequences were clustered into operational taxonomic units (OTUs), and the identity of each OTU was determined using the SILVA, GreenGenes, and RDP databases (Cole et al., 2013;DeSantis et al., 2006;Quast et al., 2012).

| Alpha diversity-all host populations
After converting raw OTU read counts to relative abundance, OTU richness (S) and inverse Simpson's diversity (D) were compared among population groups and collection sites using an analysis of variance (ANOVA; Oksanen et al., 2016).

| Beta diversity-all host populations
Beta diversity (compositional dissimilarity) was calculated among samples using the Bray-Curtis dissimilarity (BCD) calculation. BCD was calculated on two transformations of the OTU data table (relative abundance (RA-BCD) and presence-absence (PA-BCD)) to determine the relative dissimilarity associated with changes in microbial taxa abundance and taxa presence. Microbial beta diversity was compared among discrete genetic populations and geographic groups using the permuted multivariate ANOVA (PERMANOVA) function "adonis" in the R package vegan (Oksanen et al., 2016).
Pairwise differences among groups were assessed using the pairwise F I G U R E 1 Cliona delitrix is an excavating sponge in the family Clionaidae. This species has a wide geographic range that extends from Florida, USA, through the Caribbean, and along the Atlantic coast of South America (World Porifera Database) PERMANOVA function in the package RVAidememoire (Hervé & Hervé, 2014), which uses the same PERMANOVA calculations as "adonis" with a multiple comparison correction based on Benjamini and Hochberg (1995). Distance-decay relationships between community composition and the continuous genetic and geographic distances were tested using the Mantel tests and partial Mantel tests (Legendre & Legendre, 2012)

| Beta diversity-parapatric host populations
Two populations in the current study exhibit parapatry (Chaves-Fonnegra et al., 2015), while the other two populations were sampled from disparate geographic locations (Belize and Panama).
Thus, for these two geographically disparate populations, geographic and genetic effects may be confounding. To better tease

| Presence and abundance of taxa across sites and individuals
To assess the conservation of microbial associations across populations and geographic distances, we investigated the presence and abundance of taxa across collection sites and individuals. The occurrence frequency of microbial taxa was assessed across the 10 collection sites to investigate the relationship between occurrence frequency and relative abundance of taxa. We also designated a subset of taxa as "core" taxa. We defined "core" taxa as those OTUs that occurred in at least 42 of the 48 samples (~88%). This cutoff was chosen to allow possible absence of core taxa from geographically disparate samples (i.e., Panama, n = 6 or Belize, n = 4). After determining the core taxa (OTUs), we queried the presence and abundance of these core taxa in a larger microbiome dataset of sponge and seawater samples (EMP study 1740). This larger dataset was composed of samples that were concurrently sequenced by the EMP and subjected to the same postsequencing processing. This dataset of 1,227 sponge and environmental samples, which includes the samples in the current study, is described in Thomas et al. (2016).
OTUs were partitioned into four categories based on their occurrence and abundance in C. delitrix and in the larger sponge and envi- delitrix-specific (unique to C. delitrix).

| Host microbiome-final sequences and OTUs
After quality filtering, 22,348 total V4 sequences were obtained from 48 C. delitrix samples from 10 distinct sites (Table S1). The raw data for this study are available at the EMP data portal (https ://qiita. ucsd.edu/emp/study/ descr iptio n/1740; Study ID: 1740; Thomas, 2014). Before analysis, singletons and OTUs that occurred in fewer than two samples (4%) were removed, which reduced the total number of OTUs to 10,151.

| Beta diversity-discrete-geographic and genetic groups
We observed significant differences in RA-BCD associated with all discrete factors when each was considered separately (PERMANOVA; genetic population: p = .001; reef: p = .001; interaction: p = .05). However, significant overlap among independent variables was observed when all factors were considered in the same PERMANOVA (Table 2). We did observe a significant interaction between reef and genetic population (PERMANOVA; p = .05), and together, these two factors accounted for nearly half of the variance among samples (Table 2). These results may suggest overlap in the explanatory power of these factors at this spatial scale.
Pairwise PERMANOVA revealed significant differences among all genetic populations, but differences among reefs were not always significant (Table S2). For PA-BCD, we observed that a significant portion of the variance associated with each factor was associated with the presence of unique taxa (RA-R 2 /PA-R 2 ; genetic population:

| Beta diversity-continuous-geographic and genetic gradients
We observed significant correlation between RA-BCD and both geo- These partial Mantel tests revealed that geographic distance was more strongly correlated with RA-BCD than genetic distance; however, the effect of genetic distance on beta diversity remained significant in the partial Mantel test.

| Beta diversity-parapatric populations Florida Reef Tract and Bahamas
We observed significant differences in RA-BCD among genetic population groups and reefs, even when only considering parapatric populations. These two factors exhibited a significant overlap in variance (~8%) and a significant interaction (PERMANOVA; Genetic

| Presence and abundance of taxa across sites and individuals
Initially, most taxa found in the 48 C. delitrix samples occurred in a single individual ( Figure 5), which is observed by comparing OTU presence and frequency in a single reef before ( Figure 5a) and after (Figure 5c) a data cleaning step that specifically removed and surprisingly, these unique OTUs were all low abundance taxa (mean ± SD: 0.03 ± 0.01% of total community).

| Sponge-enriched OTUs
The ten sponge-enriched OTUs, all of which were Proteobacteria, comprised 3%-70% (mean ± SD: 43 ± 18%) of the total microbial community in C. delitrix. Samples collected from Panama (population D) had the lowest mean abundance for the enriched taxa (mean ± SD: 19 ± 16%). However, the Belize samples (population C; mean ± SD: 43 ± 10%) and the two populations in Florida (population A: mean ± SD: 45 ± 15%; population B: mean ± SD: 49 ± 18%) both showed higher mean relative abundances for these ten enriched OTUs (Figure 6). When we compared enriched OTU relative abundance to genetic and geographic distance variables, we found a significant independent (Mantel test; genetic dist.:  Table 4). Specifically, seven OTUs showed significant differences, one with significant variation among reef, five with significant differences among genetic populations, and one with significant differences related to both reef and genetic population (Table 4).

| Genus-and species-specific OTUs
The majority of the core OTUs were specific to the genus Cliona or the species C. delitrix when compared to a concurrently sequenced dataset of 1,227 sponge and environmental microbiomes. However, all of these OTUs were low abundance taxa and the maximum abundance of these taxa in any individual was 0.7% (Otu001869). Eleven of these low abundance OTUs appear to be conserved across some members of the genus Cliona, specifically Cliona celata. C. delitrix and C. celata were collected from opposite sides of the Atlantic and processed at separate institutions before being sent to EMP (Thomas et al., 2016). Additionally, other sponge species in this broader dataset (study 1740) were collected at the same sites as these two species,

| D ISCUSS I ON
The present study focuses on the symbiotic microbiomes of C. delitrix and represents one of the first integrative studies to link host and microbial genetics on a broad scale. Variation in the symbiotic microbial community associated with C. delitrix was related to both host population genetics (genetic group; Chaves-Fonnegra et al., 2015) and biogeography. Our analysis revealed a significant relationship between geographic distance (i.e., a distance-decay relationship) and change in microbial community composition, but only across large geographic distances (~15° latitude). Significant microbiome variation was not observed over smaller geographic distances (e.g., within the Florida Reef Tract), which is likely due to the high observed variability at the level of individual collection reef. Significant, but moderate, correlations between microbiome variation and host population genetics were observed at all geographic scales. This finding suggests that host specificity in sponge microbiomes may be more complex than previously documented and perhaps sensitive to even small, population-level genetic variation. Although considerable variation was observed among samples, several taxa were conserved across geographic and genetic groups. A total of 63 core OTUs were designated as core taxa. Ten of these core taxa were "sponge-enriched," and these taxa were consistently the most dominant members of C. delitrix microbiomes. Forty-nine core taxa were specific to the genus Cliona or to C. delitrix samples when compared to a concurrently sequenced dataset of over 1,200 sponge and environmental microbiomes. Surprisingly, these genus and speciesspecific taxa occurred at low abundance in all samples despite their persistent presence.
The interplay between host and environmental forces is a common theme in microbiome research. Some research in sponge microbiology suggested that high environmental influence and low host specificity would occur in sponges that host a low abundance of microbes (i.e., LMA species; e.g., Webster et al., 2013).
However, recent studies have observed a strong and persistent effect of host species irrespective of geography, microbial abundance, or microbial community complexity. Several studies that have observed a relationship between microbial community structure and environmental variability have also noted strong host effects (e.g., Cleary et al., 2013;Griffiths et al., 2019;Hardoim et al., 2012;Marino et al., 2017;Reveillaud et al., 2014), with most of the environmental variability being observed within species. A recent study on geographic variation in the microbiome of the sponge Ircinia campana observed significant variation in microbial communities only over large geographic gradients, but also noted the presence of location-specific microbial taxa (Marino et al., 2017).

F I G U R E 6
Bubble plot of 63 core taxa in C. delitrix samples showing variability in the mean relative abundance in core OTUs among populations of C. delitrix. Each plot represents the same 63 core taxa, and the eight most abundant taxa (in all populations) are labeled with their respective OTU identities in each panel. Each bubble represents one core OTU, and bubble size indicates the relative abundance of an OTU. Core taxa categories are represented as different colors: Blue-Environmental; orange-Sponge-enriched; green-Cliona-specific; yellow-Cliona delitrix-specific Similarly, C. delitrix microbial communities in the current study exhibited significant intraspecific variation in microbial community composition at large latitudinal gradients (~15°), while small-scale geographic effects were largely due to site-specific microbial taxa.
Large-scale geographic gradients in the current study were mostly driven by differences among samples from Belize and Panama, when compared to the eight other collection reefs that contained the two remaining genetic population groups.  observation was particularly striking when we narrowed the geographic extent of sampling to only include the two parapatric populations along the Florida Reef Tract and Bahamas. In this instance, spatial effects (distance-decay) were absent, but a moderate correlation with genetic dissimilarity remained. Previous studies have documented a connection between host genetics and associated microbial communities and shown that community structure (i.e., alpha diversity) exhibits a phylogenetic signal whereby more closely related species have more similar community diversity (Easson & Thacker, 2014;Schöttner et al., 2013;Thomas et al., 2016). However, these studies have also suggested that selection for divergent microbial community composition remains strong even among closely related species; thus, selection for specific membership in these microbiomes appears to be regulated mostly at the species level (Easson & Thacker, 2014). Our findings in this study support these previous results through both a lack of within-species differences in alpha diversity and a correlation with host genetic variation and microbial community composition. These results might indicate that as these hosts begin to genetically diverge from one another, possibly due to reproductive isolation, the composition of their associated microbial communities will begin to reflect it, even when the divergence is minimal.

| Core communities
The concept of a "core" symbiont community (taxa that are shared among individual hosts) was initially used to highlight taxa that were shared across many host species collected from different locations. For sponges, this hypothesis was proposed to support the occurrence of a uniform microbial community among diverse sponge species (Hentschel et al., 2002). As research on sponge microbiomes has progressed, the term "core" has varied greatly in its definition and use (Astudillo-García et al., 2017;Schmitt et al., 2012). In the current study, we used it to describe symbiont taxa that occurred in at least 88% of collected samples. With the increasing sequencing depth afforded by next-generation sequencing techniques, we now understand that many "core" symbiont taxa are detectable in the environment and thus likely enriched from it (i.e., horizontally acquired; Taylor et al., 2013).
Ten OTUs in the current study were both core and sporadically detected in other sponge hosts or environmental samples in the sponge microbiome project (Easson & Thacker, 2014;Thomas et al., 2016). Eight of these OTUs composed the dominant members transmitted. If these taxa are horizontally transmitted, C. delitrix juveniles must have a mechanism for selecting these specific taxa out of the vast number that exist in the environment similar to those mechanisms in other symbiotic systems (Baker et al., 2019;Mcfall-Ngai, 2014;Nyholm & McFall-Ngai, 2004).
The majority of the core OTUs were specific to either the genus Cliona or C. delitrix. These OTUs, however, all occurred at values of low relative abundance, with a maximum relative abundance of 0.7% in a single sample. The observed low relative abundance is most likely below the detection limit of many earlier methods for assessing symbiont communities (e.g., clone libraries). Low abundance of these taxa could be interpreted as indicating that these taxa play relatively minor ecological roles in the holobiont, but the persistence of these rare taxa in hosts separated by wide geographic expanses and evolutionary time suggests otherwise. Some recent research has focused on microbial taxa found at low abundance, terming them the "rare biosphere" and explored the forces driving the dynamics of the rare biosphere. In some systems, rare taxa are disproportionately active and can exhibit widely varying temporal profiles in abundance and activity (Sogin et al., 2006). Much of the research into the rare biosphere has focused on free-living systems (e.g., bacterioplankton), but in a symbiotic system, constraints on microbial abundance could be quite different due to host-symbiont interactions. The abundance of rare taxa could also be driven by top-down forces such as phage predation or host immune responses, which require taxa to remain rare to avoid detection (reviewed in Lynch & Neufeld, 2015). If this scenario explains the rare taxon abundance in the current study, it might signal that at least some of these taxa are parasitic cheaters that have found a way to persist in association with the host sponge.
Rare taxa can also represent a genetic "seed bank" of ecological potential by containing a diverse cache of metabolic machinery.
In the context of the current study, this scenario would suggest that these rare taxa are important in specialized situations, such as other life stages (e.g., larval stages) or at time points crucial to the species (e.g., during spawning). An environmental or simply stochastic explanation for the presence of these rare taxa seems implausible given their persistence across the Caribbean in C. delitrix and in some cases across the Atlantic Ocean in C. celata, as well as their conspicuous absence in other host taxa (even sympatric species), but additional research is needed to adequately address this hypothesis. Sequencing the microbiomes of C. delitrix larva and adults during spawning would likely provide additional evidence as to the role of these rare but persistent taxa.
Our study explored variability in sponge microbial communities associated with variation in sponge microbiomes and, similar to previous studies, found an interplay between host and geographic forces in structuring the microbiome of C. delitrix. Variability among C. delitrix microbiomes was apparent at large geographic scales, while a significant but moderate host population genetic correlation was observed independent of geography. Although variation among samples was observed, several bacterial taxa were consistently found in C. delitrix samples. The most dominant core taxa were also observed at low abundance in the environment, while the Cliona-specific taxa all occurred at low abundance within the host sponge. To date, microbiome divergence has been observed at the level of host species with limited evidence of divergence among populations (Griffiths et al., 2019;Swierts et al., 2018). The results of the current study add to this body of evidence and indicate that microbiome divergence can be observed at the population genetic scale, prior to the onset of speciation.

ACK N OWLED G EM ENTS
We thank all the scientists that participated in the sampling and submission of sponge specimens to the Earth Microbiome Project (EMP), and in particular, we thank Torsten Thomas and Nicole Webster for coordinating this effort. We also thank our colleagues

AUTH O R CO NTR I B UTI O N S
A. Chaves-Fonnegra designed the study and prepared the samples for sequencing. C.G. Easson and R.W. Thacker processed the sequence data, designed and implemented the data analysis, and wrote the first draft of the manuscript. C.G. Easson, A. Chaves-Fonnegra, J.V. Lopez, and R.W. Thacker contributed to the writing and editing of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequence data used in this study are available through the Qiita server at (https ://qiita.ucsd.edu/emp/ , Study ID -1740 (https ://qiita. ucsd.edu/emp/study/ descr iptio n/1740; Thomas, 2014).