Introduced ascidians harbor highly diverse and host-specific symbiotic microbial assemblages

Many ascidian species have experienced worldwide introductions, exhibiting remarkable success in crossing geographic borders and adapting to local environmental conditions. To investigate the potential role of microbial symbionts in these introductions, we examined the microbial communities of three ascidian species common in North Carolina harbors. Replicate samples of the globally introduced species Distaplia bermudensis, Polyandrocarpa anguinea, and P. zorritensis (n = 5), and ambient seawater (n = 4), were collected in Wrightsville Beach, NC. Microbial communities were characterized by next-generation (Illumina) sequencing of partial (V4) 16S rRNA gene sequences. Ascidians hosted diverse symbiont communities, consisting of 5,696 unique microbial OTUs (at 97% sequenced identity) from 44 bacterial and three archaeal phyla. Permutational multivariate analyses of variance revealed clear differentiation of ascidian symbionts compared to seawater bacterioplankton, and distinct microbial communities inhabiting each ascidian species. 103 universal core OTUs (present in all ascidian replicates) were identified, including taxa previously described in marine invertebrate microbiomes with possible links to ammonia-oxidization, denitrification, pathogenesis, and heavy-metal processing. These results suggest ascidian microbial symbionts exhibit a high degree of host-specificity, forming intimate associations that may contribute to host adaptation to new environments via expanded tolerance thresholds and enhanced holobiont function.

The diversity metrics OTU richness (S), the Shannon-Weaver diversity index (H′), and the Chao1 index were all significantly different (P < 0.05) among the three ascidian species for the overall data partition (Supplementary  Table S1). Pairwise comparisons for the overall data partition indicated that P. zorritensis exhibited significantly greater observed (S) and expected (Chao1) microbial richness than P. anguinea, which exhibited a significantly greater richness than D. bermudensis (Supplementary Table S1). P. zorritensis exhibited a significantly greater Shannon diversity (H′) than D. bermudensis, but was not significantly different from P. anguinea at the overall level (Supplementary Table S1). No significant differences in the evenness or Simpson diversity indices were detected among the microbial communities of the three species at the overall level (Supplementary Table S1).
Symbiotic microbial community structure. An assessment of microbial community structure among the different sources revealed that communities clustered in response to source (Fig. 3), with significant differences between the microbial communities of D. bermudensis, P. anguinea, and P. zorritensis, as well as ambient seawater (PERMANOVA, P < 0.01), and 78% of the community structure variation explained by source. Pairwise comparisons revealed significant differences (PERMANOVA, P < 0.015) between microbial communities of each host source, and ambient seawater. Within each source, seawater exhibited less variability in microbial community structure among replicates (average similarity = 79.0%) compared to the microbiomes of D. bermudensis, P. anguinea, and P. zorritensis (average similarities = 55.9%, 64.0%, and 57.3%, respectively); however, no significant differences in dispersion were detected among sources (PERMDISP, P = 0.221).
Examining the ascidian microbial communities exclusively, community structure was again significantly different among species (PERMANOVA, P < 0.01), with host species explaining 72.6% of the variation in overall community structure. Partitioning of the ascidian microbial community data into abundant (>0.1% relative abundance) and rare (≤0.1%) components resulted in 751 abundant OTUs that accounted for 96.9% of the total sequences and 4,945 rare OTUs that accounted for 3.1% of the total sequences, post singleton removal. Both the abundant and rare OTU data partitions indicated significant differences in microbial community structure  (PERMANOVA, P < 0.01), with source explaining 75% of the variation in the abundant community structure, but only 18.6% of the variation in the rare community structure. Multiple pairwise comparisons revealed significant differences between the microbial communities of all ascidian host sources, for all three data partitions ( Table 1). The permutational multivariate analysis of dispersion across groups was not significant for either the overall or abundant OTU data partitions (PERMDISP, P = 0.86 and P = 0.818, respectively), but was significant for the rare OTU data partition (PERMDISP, P < 0.01), with all subsequent pairwise comparisons likewise proving significant (Table 1). Additionally, for the rare OTU data partition, greater dissimilarity was observed between replicates within each species (Fig. 4c) compared to the overall or abundant data partitions (Fig. 4a,b).
The congeneric species Polyandrocarpa zorritensis and P. anguinea exhibited more similar microbiomes compared to the distantly related host D. bermudensis (Fig. 4a-c), though the differences observed between these closely related ascidians were still significant ( Table 1). P. anguinea and D. bermudensis were the most dissimilar, with an average dissimilarity of 91.9, 91.9 and 92.1% for the overall, abundant OTU and rare OTU datasets, respectively. D. bermudensis and P. zorritensis had an average dissimilarity of 76.5, 75.4 and 91.8% across the three dataset partitions respectively, while P. anguinea and P. zorritensis had an average dissimilarity of 75.1, 74.4 and 82.5%. SIMPER analyses identified the top OTUs responsible for the observed dissimilarity between ascidian microbiomes ( Table 2). For the overall microbial community and the abundant OTU data partition, a small number of symbiont OTUs (four or five) contributed to the majority (>50%) of the observed dissimilarity among hosts ( Table 2). In contrast, the differences between the rare microbial communities of the three species or sources were driven by the presence or absence of a great number of microbial taxa (569 to 831 OTUs).
Core symbiont community composition and diversity. Core symbiont communities were identified for each ascidian species (i.e. OTUs shared by all five replicates) and for all ascidian hosts (i.e. OTUs shared by all five replicates from all three species). The core community in D. bermudensis consisted of 147 OTUs that comprised 93.3% of all sequences, with the top 15 most prevalent core OTUs accounting for the vast majority (84.8%) of all sequences (Supplementary Table S2 Table S2). Across all ascidian hosts, 103 universal core OTUs were identified. Although  the relative abundance of these universal core OTUs varied among hosts (Table 3), these 103 OTUs represented just 1.81% of all detected ascidian OTUs, but 79.1% of all sequences. In fact, the top 15 of these 103 OTUs represented 68.9% of all sequences and corresponded to one archaeal (Crenarchaeota) and 14 bacterial symbionts. Of the bacterial OTUs, seven were affiliated with Alphaproteobacteria, four with Gammaproteobacteria, two with   Table 2. Taxonomic classification and relative abundance of symbiont OTUs contributing to the observed dissimilarity in microbial community structure (50% cumulative, SIMPER analyses) between D. bermudensis (DB), P. anguinea (PA), and P. zorritensis (PZ). *K kingdom, P phylum, C class, O order, F family, G genus, S species.
Bacteroidetes, and one was unclassified at the phylum level (Table 3). This one unclassified OTU represented 98.3% of all unclassified bacterial sequences in D. bermudensis and 94.0% of unclassified sequences for all three species together.

Discussion
With the rapid advances in sequencing technology, our understanding of the complex interactions between microbial symbionts and their hosts continues to expand. Here, we examined the microbial communities of three widely introduced ascidian species, as well as those of the surrounding ambient seawater, as a first step in understanding the role that microbial symbionts may play in the successful establishment of these ascidian species in North Carolina. The microbial communities of D. bermudensis, P. anguinea, and P. zorritensis exhibited significant differences in both community structure and diversity. Microbial communities in the three species were also clearly different from those of the surrounding seawater, as previously reported for other ascidian species 32,34,36 .
Over half (53.2%) of all OTUs obtained herein were present only in ascidian hosts and not detected in ambient seawater. Comparisons of microbial community structure among the three ascidian species revealed a high degree of host-specificity, with significant differences in microbial community structure detected among all hosts. Interestingly, the microbiomes of the two congeneric ascidian species were more structurally similar to one another than to D. bermudensis, providing further support for the host-specificity of the ascidian-microbe interactions described herein. Within each ascidian species, broad taxonomic comparisons (phylum-level) indicated a relatively conserved symbiont composition across replicates. At the OTU-level, replicates within the same species exhibited at least 55% average similarity for both the overall and abundant OTU data partitions. Intra-specific variation was explained primarily by the presence or absence of numerous rare OTUs, for which the average similarity within hosts ranged from 13.7% (D. bermudensis) to 25.1% (P. zorritensis).
The high degree of host-specificity in ascidian microbiomes observed herein further supports previous findings indicating that each ascidian tunic possesses unique internal features or physiological characteristics that provide favorable conditions (i.e. microniches) for supporting distinct microbial communities 32,37 . Both vertical and horizontal symbiont transmission may be involved in the establishment and maintenance of these distinct microbiomes 31 . A high degree of host-specificity indicates that some degree of vertical transmission occurs in colonial ascidians, as further supported by microscopic observations of bacteria in ascidian larvae 34,38 . Ascidians also appear to selectively accumulate specific microbes from the surrounding seawater, as rare seawater microbes have been shown to be enriched in ascidian tunics 31 . More relevant to this particular study, the ascidian Polyandrocarpa zorritensis has been shown to accumulate pollution indicator microbes directly from seawater 17 , suggesting horizontal transmission may also play a role in construction of the microbiome of introduced ascidians. Indeed, in the current study, most universal ascidian symbiont OTUs (80.3%) were also detected in the seawater, with some of these OTUs matching identically (100% pairwise identity) to sequences previously reported from seawater samples or sponge tissue within the same region 39 . Given the significant structural differences between the microbial communities of the four sources, these similarities are especially striking. Taken together, our findings and those of previous studies suggest that microbial symbionts within ascidian species are sourced both vertically via transmission from parent colonies and horizontally through uptake from the local environment, supporting an environmentally "leaky" model of vertical transmission in ascidians that has previously been proposed for marine sponges and mussels 40 . Such a strategy could greatly benefit ascidians introduced into a new habitat by arriving seeded with beneficial symbionts (e.g. those contributing to host primary metabolism), yet also being capable of acquiring locally-adapted (e.g. pollution tolerant) symbionts. Thus, both the microbial   inheritance of each introduced ascidian and the microbial residents of each susceptible habitat may represent important factors in the successful establishment of introduced species. The host-specific core communities of all three ascidian species were comprised of a small number of OTUs with disproportionate abundances, accounting for the vast majority (85 to 93%) of all sequences obtained for D. bermudensis, P. anguinea, and P. zorritensis. A prior study of the microbiome in the globally introduced ascidian Styela plicata found a similarly disproportionate representation of a few core OTUs 35 . The existence of core OTUs suggests that these particular symbionts may play an important role in the survival of a particular host species, especially OTUs with high relative abundances, as they have been linked to increased host functionality in marine invertebrates 41 . However, it should be noted that the rare microbiome has been shown in some cases to dominate certain processes, such as nutrient cycling, despite low microbial abundances [42][43][44] , and thus the functional implications of rare symbionts within the ascidian microbiome should not be disregarded. Microbial symbionts may contribute to ascidian host survival in several ways, including supplementing host nutrition (e.g. Prochloron in some didemnid species 45 ), providing defense through the production of secondary metabolites 24,27,30 , or providing antifouling capabilities 20,46 . Microbial symbionts have been linked to the success of invasive species within a terrestrial framework, with the presence of specific endosymbionts being directly related to increases in invasive plant biomass 47 . Whether microbial symbionts play similar roles in facilitating successful marine invasions is still poorly understood, and the potential benefits of those symbionts to the holobiont are not well-characterized.
The identification of 103 universal core OTUs (i.e. detected in all replicates of all species) herein is especially interesting in the context of invasiveness potential. The top 15 universal core OTUs varied in relative abundance among hosts, yet accounted for a large proportion of all sequences obtained, with a potentially dominating effect on symbiotic functions within introduced ascidians. Universal core OTUs were affiliated with archaeal and bacterial taxa that exhibit functional capabilities of potential benefit to the holobiont, including nitrogen cycling, protection against UV radiation, and heavy metal processing. The most prevalent archaeal symbiont in the universal core community was assigned to the genus Nitrosopumilus, an ammonia-oxidizing archaeal lineage 48 that has been described in association with other ascidian hosts 31,32 . Ascidians expel the majority of their nitrogenous waste in the form of ammonia 49 , and high concentrations of ammonia-rich waste within the tunic may support growth of Nitrosopumilus symbionts. In return, ammonia-oxidizers may benefit ascidian hosts by processing and removing nitrogenous waste. Similarly, two lineages of proteobacterial denitrifiers (Rhodobium orientis 50 and Novispirillum sp. 51 ) were identified as top universal core OTUs within the three ascidian species. Denitrification involves the microbial reduction of nitrate and may also contribute to successful ascidian introductions within polluted harbor systems by increasing the tolerance threshold of the holobiont to higher nitrogen pollution levels. Harbor systems, in which ascidians are most commonly introduced, are notoriously polluted due to their close proximity to anthropogenic development, contamination from high boat traffic (e.g. fuel combustion and spills 52 , leaching of anti-fouling hull paints 53 ), and their relative isolation from the open ocean. As a result, pollutants become concentrated within harbor systems, including toxic heavy metals 54 , hydrocarbons 54 , and nitrogenous pollution from sewage effluence 55 . This reduced circulation and increased pollution in harbors has been linked to permanent changes in ascidian community structure 56,57 , and localized increases in introduced ascidian species have been observed in harbor systems at the apparent expense of native ascidian populations 7 . Additional studies targeting nitrifying and denitrifying symbiotic microbial guilds are needed to assess their contribution to the invasiveness potential of introduced ascidian species in polluted harbor environments, as well as the degree to which ammonia oxidization and denitrification capabilities may enhance holobiont metabolism. Notably, the simultaneous occurrence of symbiont taxa capable of aerobic and anaerobic nitrogen transformation pathways suggests that these symbionts may occur in distinct regions of the inner tunic, or are separately maintained by anoxic and normoxic conditions generated by the process of filter feeding, as has been suggested for marine sponges 58 .
In addition to putative roles in nitrogen cycling, universal core OTUs in introduced ascidians also matched to bacterial taxa capable of carotenoid production and heavy metal processing. Two universal core symbionts, R. orientis and Reichenbachiella sp., produce carotenoids 50,59 , which have been suggested to provide UV protection in other marine invertebrate hosts 60 and could be especially critical for survival in shallow harbor environments. Vibrio shilonii, another universal core symbiont, may provide heavy metal processing capabilities to ascidian hosts, as symbionts in the genus Vibrio found within marine invertebrates have shown resistance to heavy metal toxicity 61 . More specifically, a comparative genomic analysis of Vibrio spp. identified a gene coding for the copper chaperone protein CopZ within V. shilonii, suggesting this particular symbiont could metabolize copper 62 . This heavy metal processing capability could help explain how introduced ascidians continue to thrive in polluted harbor regions, even as native species disappear 7,8 . However, the extent to which these capabilities would affect the holobiont fitness and long-term colonization success is unclear. In addition, V. shilonii (also called V. mediterranei 63,64 ) is a known coral pathogen that has been linked to coral bleaching events 65 . The presence of documented invertebrate pathogens as core symbionts in introduced ascidians provides preliminary evidence for an additional negative ecological impact of species introductions to native ecosystems. Putative pathogenic microbial symbionts have been previously described within other invertebrates 66,67 , including ascidians 68 , suggesting that invertebrate microbiomes may act as reservoirs for those pathogens. Indeed, V. shilonii was able to survive lethal winter temperatures within invertebrate hosts, which subsequently served as pathogen vectors 66 . Clearly, the full ecological impacts of species introductions on ecosystem health are not well understood and merit further investigation. As global changes in temperature regimes continue to favor introduced species to the detriment of native ones 69 , including ascidians 70 , a greater understanding of microbiome contributions to invasiveness potential and ecosystem damage will help clarify pathways for invasion success and prevention.

Materials and Methods
Sample collection and DNA extraction. Three species were targeted in this study: Polyandrocarpa zorritensis, P. anguinea, and Distaplia bermudensis. All three species are considered introduced in North Carolina 71 . P. zorritensis was first described by Van Name in Peru 72 and exhibits a worldwide distribution, including Panama 73 , the southeastern United States 71 , Italy 74 , and Spain 75 . P. anguinea was first described by Sluiter in Mozambique 76 , and has since been described in Brazil 77 , Panama 73 and the southeastern US 71 . D. bermudensis was first described by Van Name in Bermuda 78 , and has since been observed in Puerto Rico 79 , the southeastern US 71,79 , Spain 80 , and Italy 81 .
Ascidian and ambient seawater samples were collected at two nearby sites (<0.15 km apart) in Wrightsville Beach, NC. Replicates of D. bermudensis (n = 5) were collected in June 2015 from the Bridge Tender Marina (34°13′06″N 77°48′47″W). Replicates of P. anguinea (n = 5) and ambient seawater (n = 1) were collected in September 2015 from Wrightsville Beach Marina (34°12′59.6″N 77°48′45.7″W). Replicates of P. zorritensis (n = 5) and ambient seawater (n = 3) were collected in October 2015 from the Wrightsville Beach Marina site. Ascidian identifications in the field were based on morphological characteristics, and confirmed in the lab via barcoding of a fragment of the mitochondrial gene cytochrome c oxidase subunit I (COI; see below). All samples were collected at <1 m depth and at least 5 m apart to minimize the probability of sampling clones. Samples were housed separately in Whirl-Pak bags containing ambient seawater. Once in the laboratory (less than 10 km away from the sampling sites), each ascidian sample was dissected in half. One half was immediately fixed in 10% formalin and filtered seawater for morphological analysis, and stored at room temperature. The other half of each colony was rinsed with filtered seawater and stored in absolute ethanol at −20 °C for molecular analysis. Ambient seawater samples (500-ml each) were concentrated onto 0.2-μm filters and immediately frozen at −80 °C until analysis. Prior to DNA extraction, ethanol-preserved ascidian samples were dissected under a stereomicroscope into zooid and tunic fractions for host barcoding and symbiont characterization, respectively. DNA extractions of the ascidian zooids, inner tunics (i.e., tunic tissue not in contact with either the surrounding seawater or the zooids), and seawater filters were performed with the DNeasy Blood and Tissue kit (Qiagen).

Ascidian genetic barcoding.
To confirm morphological identifications of collected ascidians, a ca. 600-bp segment of the mitochondrial COI gene was amplified from zooid extractions using the universal LCO1490 and HCO2198 primers 82 or the ascidian-specific primers Tun_forward and Tun_Reverse2 83 . PCR amplifications were conducted with a total volume of 25 μl, including 5 pmol of each primer, 2X MyTaq HS Red Mix (Bioline), and 1 μl (ca. 10 ng) of template DNA. The thermocycler program consisted of an initial denaturation at 95 °C for 1 min; 35 cycles of 95 °C for 15 s, 45 °C for 15 s, and 72 °C for 10 s; and a final extension at 72 °C for 1 min. Sequence reactions were conducted with BigDye Terminator v. 3.1 (Applied Biosystems) and the same forward and reverse primers utilized in the initial amplifications. PCR products were purified via isopropanol precipitation or BigDye Xterminator (Applied Biosystems) and sequenced on an ABI 3130xl genetic analyzer (Applied Biosystems) or an AB 3500 genetic analyzer (Applied Biosystems), both available at UNCW Center for Marine Science. Raw sequence reads were processed in Geneious version 8.02 84 by aligning forward and reverse reads to create consensus sequences. These consensus sequences were compared to sequences within the GenBank database using a nucleotide-nucleotide BLAST search 85 (BLASTn) and archived in GenBank under the accession numbers MF034529 to MF034540.
Next-generation sequencing of microbial symbionts. To characterize microbial communities in ascidian tunics and seawater samples, a ca. 300 bp fragment (V4 region) of the 16S rRNA gene was amplified and sequenced from tunic and seawater extracts using the universal bacterial/archaeal forward primer 515f and reverse primer 806r 86 . Initial PCR amplifications were used to confirm the viability of DNA extractions, with a thermocycler program consisting of an initial denaturation at 95 °C for 2 min; 35 cycles of 95 °C for 15 s, 50 °C for 15 s, and 72 °C for 20 s; and a final extension at 72 °C for 2 min. DNA extracts were sent to Molecular Research LP for amplification, library construction and multiplexed sequencing of partial 16S rRNA gene sequences with the 515f and 806r primers on an Illumina MiSeq platform. Raw sequence data were deposited in the Sequence Read Archive of NCBI (accession no. SRP106072).
Next-generation sequence data processing. Raw sequences were processed in the mothur software package 87 following a modified version of the Illumina MiSeq SOP pipeline 88 as described previously 39 . Briefly, raw sequences were quality-filtered and aligned to the SILVA reference database (v119), putative chimeric sequences were removed via self-reference searching with UChime 89 , sequences were classified using a naive Bayesian classifier and bootstrap algorithm for confidence scoring 90 based on the improved Greengenes taxonomy 91 , and nontarget sequences (chloroplasts, mitochondria, and eukarya) and singletons were removed from the data set. High-quality sequences were assigned to operational taxonomic units (OTUs) in mothur based on 97% sequence identity and the average neighbor clustering algorithm, and the taxonomic classification of each OTU was determined by majority consensus 87 . To standardize sampling depths (i.e. number of sequence reads) among the different replicates, each data set was subsampled to the lowest read count (n = 24,595) from the final shared file, and all subsequent data analyses were based on the resulting subsampled data sets.
Microbial community diversity analysis. To compare microbial community diversity among the four sources (three ascidian species and ambient seawater), alpha diversity metrics for OTU richness and evenness were calculated in mothur, including observed OTU richness (S), expected OTU richness (Chao1), the Simpson evenness index (E1/D), the inverse Simpson index (D) and the Shannon-Weaver diversity index (H′). Analyses of variance (ANOVA) were used to statistically compare the diversity indices among the four sources, and Tukey's honest significance difference (HSD) tests were performed for multiple pairwise post hoc comparisons of means.
SCieNTifiC REpoRts | 7: 11033 | DOI:10.1038/s41598-017-11441-4 Microbial OTU overlap and core symbiont communities. Venn diagrams were constructed in mothur to visualize OTU overlap among sampling sources. Core OTU communities were identified for each host ("host-specific" core) and across all ascidian hosts ("universal" core). Host-specific core communities included OTUs present in all replicate samples of each host species. Universal core communities included OTUs present in all replicates of all host species. Representative sequences from the 15 most abundant core OTUs (host-specific and universal) were analyzed using BLASTn to identify related sequences in GenBank and compare corresponding sources. Sequence comparisons were based on highest percent identity matches and identified sources.
Microbial community structure analysis. To assess host-specificity and compare the structure of microbial symbiont communities among the four sources, Bray-Curtis similarity matrices were created based on OTU relative abundances in PRIMER (version 6.1.11) and visualized in multi-dimensional scaling plots. Permutational multivariate analyses of variance (PERMANOVA) were used to statistically compare the structure of microbial communities among the four sources, including multiple pairwise comparisons for significant main test PERMANOVA results. Pairwise comparisons were corrected via the Benjamini-Yekutieli (B-Y) false-discovery rate control and an experiment-wise error rate of α = 0.05 92 . Permutational multivariate analyses of dispersion (PERMDISP) were performed to verify that significant PERMANOVA results stemmed from structural differences, not unequal dispersion of variability, among groups. Sequence data were divided into abundant and rare-OTU partitions using a 0.1% relative abundance threshold 43 , resulting in a cutoff value of 24 sequences ("abundant" OTUs contained >24 sequences, "rare" OTUs contained ≤24 sequences). Bray-Curtis similarity matrices were created for each data partition and PERMANOVA were used to test for differences in the structure of abundant and rare members of the microbiome.