The Microbiome of the Worldwide Invasive Ascidian Didemnum vexillum

All multicellular organisms, including ascidians, host diverse microbial communities that are essential for their evolution. The global invader Didemnum vexillum is a colonial species native to Japan with two main genetic clades, A (the only invasive) and B, which provides a unique opportunity to assess if the microbiome remains stable in the colonization process or shifts according to local environment. We have analyzed, using 16S amplicon sequencing, the microbiome of 65 D. vexillum colonies from 13 populations worldwide including the two clades in the native area, plus samples from a congeneric species and seawater from one of the localities. We found 3,525 zero-radius operational taxonomic units (ZOTUs) in D. vexillum, belonging to 36 bacterial and 3 archaeal phyla. The microbiome of this species had a markedly different composition from surrounding seawater and from the congeneric species. For the globally invasive clade A, we found 3,154 ZOTUs, and 8 of them were present in all colonies, constituting a core microbiome with high abundance (69.57% of the total reads) but low diversity (0.25% of the total number of ZOTUs). The variable component was quantitatively much less important but comprised a highly diverse assemblage. In a multiple regression model, global microbiome structure correlated with differences in temperature range across localities and also with geographic distances, pointing to horizontal acquisition of the symbionts. However, the ascidian may have a strong capacity to select and enrich its microbiome, as we found that the most abundant ZOTUs from tunic samples had low abundance in seawater samples from the same locality. The microbiome structure also correlated with the genetic distances between colonies obtained in a previous genome-wide analysis, suggesting some potential for vertical transmission. In geographically restricted comparisons, temperature and genetic makeup, but not geography, explained microbiome structure. The combination of a quantitatively dominant core component and a highly diverse variable fraction in the microbiome of D. vexillum can contribute to the success of this global invader in different environments.


INTRODUCTION
All multicellular organisms can be considered holobionts hosting diverse microbial communities that are essential for the evolution of the organism (Uriz et al., 2012;Zilberg-Rosenberg, 2016, 2018;Leite et al., 2018;Pollock et al., 2018;Cleary et al., 2019). The ubiquity and importance of bacteria in the living world have been recently showcased by the application of new genomic technologies (McFall-Ngai et al., 2013). Holobionts rely on their microbiome for vital processes (Mandrioli and Manicardi, 2013), which led to the concept of the hologenome (the collective genomes of the holobiont) as a unit of selection (Zilber- Rosenberg and Rosenberg, 2008;Rosenberg and Zilberg-Rosenberg, 2018). Numerous studies used different methods to determine how the microbiome provides mutualistic benefits to the host. Some of these benefits are nutrient fixation or improved metabolism (Newton et al., 2008;Barott et al., 2011;Kellogg, 2019), protection against predation or other competitors (Barott et al., 2011;Kwan et al., 2012), and adaptation to diverse environments (Rosenberg et al., 2010;Mortzfeld et al., 2016;Ziegler et al., 2017;Morrissey et al., 2019;Wilkins et al., 2019).
Ascidians are a group featuring many notorious introduced species around the world (Lambert, 2007;Shenkar and Swalla, 2011;Zhan et al., 2015). They have several characteristics of successful colonizers, such as a rapid growth rate, high fecundity, and short time to maturity (e.g., Rius et al., 2009;Pineda et al., 2013;Casso et al., 2018). They are also known to adapt to a range of temperatures, salinities, and pollution levels (Naranjo et al., 1996;Pineda et al., 2012;Nagar and Shenkar, 2016;Rocha et al., 2017). Colonial ascidians can also reproduce asexually by fragmentation and fuse to form chimeras with zooids of different genotypes, which may enhance the colony's adaptive capacity (Ben-Shlomo, 2017;Casso et al., 2019a). For introduced marine species in general, it has been pointed out that the microbial symbiont community can contribute to the survival and adaptation of the holobiont during the invasion process. Several mechanisms have been proposed, such as heavy metal detoxification (Aires et al., 2016), production of toxic chemicals that can increase competitive ability on a new habitat (Amsellem et al., 2017), or enhanced nutrient fixation (Arnaud-Haond et al., 2017). Likewise, broadly introduced species face diverse temperature regimes, and microbial symbionts have been shown to play a role in adaptation of invertebrates to temperature changes (e.g., in corals: Ziegler et al., 2017;Osman et al., 2020). In this context, the tunic ascidian microbiome has been suggested to enhance to the adaptive capacity and invasion success of this group (Evans et al., 2017;Novak et al., 2017;Dror et al., 2019). (2002) is an invasive colonial ascidian, native from the NW Pacific, with two main mitochondrial clades, named clade A and clade B, of which only clade A has spread in temperate seas worldwide (New Zealand, East and West coasts of North America, Atlantic and Mediterranean coasts of Europe, and South East of Russia: Lambert, 2009;Stefaniak et al., 2012;Ordóñez et al., 2015;Zvyagintsev et al., 2016). This clade shows a strong population structure, with three main invasive genetic groups and a high genetic differentiation between populations (Casso et al., 2019b). Environmental conditions are different between the main regions of distribution of the species and even between localities in the same region (e.g., Mediterranean and Atlantic localities in Europe). Therefore, the ascidian faced different environmental stressors during the introduction events; thus, quick responses such as those mediated by their microbiota could enhance adaptation and invasive potential. In the present study, we (i) describe and compare the microbiome of D. vexillum at different levels (species, mitochondrial clades, localities, and invasive genetic groups); (ii) identify the core microbiome of the invasive clade; and (iii) evaluate the relationship between microbiome composition and several environmental and genetic variables. We had the unique opportunity to assess individually the relationship between host genetics and microbiome as the same colonies analyzed in the present study were previously genotyped at a panel of genome-wide loci (Casso et al., 2019b).

Sampling
A total of 65 colonies of D. vexillum from 12 locations spanning the native (NW Pacific) and the known introduced distribution of the species were sampled during 2015 and 2016 in marinas and artificial substrates (Table 1 and Figure 1). Sequences of the specimens showed that five colonies were D. vexillum clade B from Sugashima (Japan), whereas the remaining 60 colonies were D. vexillum clade A from different temperate regions of the world (Casso et al., 2019b). Thus, only in one of the two populations (Sugashima and Aomori) corresponding to the native range of the species did we find the two described mitochondrial clades in sympatry. Moreover, five colonies of a different species (Didemnum sp.), likely a new species (authors' unpublished research), were sampled in a marina in Roscoff (France) and included in the present study for comparative purposes ( Table 1). All 70 colonies were collected in marinas and other artificial substrates, at least 2 m apart from each other to avoid pseudoreplication , and preserved at −20 • C in 96% ethanol. Samples were coded by geographic region as defined in previous studies Casso et al., 2019b). Clade A introduced samples were assigned to the previously defined genetic groups (Casso et al., 2019b): Europe (EUR), West North America (WNA), and New Zealand plus East North America (NZ-ENA).
Additionally, three seawater samples of 1 L were collected adjacent to the sampling site in the Ebro Delta location ( Table 1). They were separately pre-filtered using 47-mm-diameter FIGURE 1 | Map of the localities sampled in the present study. Colors follow the palette used in the nmMDS plots. The map has been drawn with the R package "rworldmap (https://cran.rproject.org/web/packages/rworldmap/index.html)." For each location, five colonies were sampled, except in Sugashima, where 10 colonies (5 of each clade) were sampled. Geographic region (Jap, Japan; NZ, New Zealand; ENA, East North America; WNA, West North America; EUR, Europe), genetic group determined in a previous study (Casso et al., 2019b), location, coordinates, and species (mitochondrial clade is indicated for D. vexillum samples) of the analyzed colonies are listed. *Locality where the three seawater samples were collected.
polycarbonate 5 µm pore-size filters to eliminate large particles and then filtered using a total of four 47 mm diameter polycarbonate 0.2 µm pore-size filters per sample (250 mL of seawater per filter). All 0.2 µm pore-size filters were preserved in 500 µL of lysis buffer (Tris 50 mM pH 8.3, EDTA 40 mM, sucrose 0.75 M) at −20 • C.

DNA Extraction and Sequencing
From each colony, a fragment of ca. 4-16 mm 2 of distal layer of the tunic (up to 0.5 mm thick) was cleaned of zooids and used for DNA extraction. We avoided using basal surfaces in contact with the substrate to prevent epibiont contamination (see section Discussion). We used the QIAamp mini kit (Qiagen, Venlo, The Netherlands) and followed the manufacturer's recommendations, with a final elution in 80 µL of AE buffer. DNA extractions of the seawater filters were performed using both the filter and the buffer. As described in Turon et al. (2018), membranes were enzymatically digested with lysozyme, proteinase K, and sodium dodecyl sulfate. DNA was extracted with phenol:chloroform-isoamyl alcohol (25:24:1, vol/vol/vol) and chloroform:isoamyl alcohol (24:1, vol/vol). Purification and concentration of the DNA was carried out with Amicon Ultra 4 Centrifugal Filter Units-100,000 NMWL (Merck Millipore, Burlington, Massachusetts, United States).
The average concentration of DNA in samples was of 13.17 ng/µL (SE ± 1.08) in a volume of 20 µL of AE buffer. A fragment of the V4 region of the bacterial 16S rRNA gene was amplified using the F515/R806 primers (Caporaso et al., 2011) and paired-end sequenced (2 × 250 bp fragments) in a MiSeq Illumina platform, at the Research Technology Support Facility of the Michigan State University (USA). All sequences obtained were deposited in the Short Read Archive of the NCBI under Bioproject PRJNA608680.

Sequence Data Processing
Raw reads were processed using USEARCH v.10.0.240 (Edgar, 2010) and VSEARCH v.1.4.1 (Rognes et al., 2016). First, paired-ends were merged (using USEARCH fastq_mergepairs command: fastq_maxdiffs = 20, setting the maximum number of mismatches in the alignment to 20 bp), and then they were truncated to 250 bp (VSEARCH fastq_filter: fastq_trunclen = 250, fastq_maxee = 0.25). Only the remaining quality-filtered sequences were used to find the set of unique sequences (USEARCH fastx_uniques). Unique sequences were then denoised using the UNOISE algorithm (Edgar, 2016), which performs error correction and chimera removal on amplicon reads (USEARCH unoise3: unoise_alpha = 2 and minsize = 8 as per default settings), to obtain a list of zero-radius operational taxonomic units (ZOTUs) (100% identity; Edgar and Flyvbjerg, 2015;Edgar, 2016). Using an identity of 100% has been shown to be optimal for the V4 16S rDNA sequences (Edgar, 2018). Finally, we constructed the final table of ZOTU abundances by mapping the reads to the list of ZOTUs using the USEARCH otutab command with the default identity threshold of 0.97. All subsequent analyses were based on this ZOTU table.
Taxonomic assignment of each ZOTU was performed using the SILVA_132_SSURef_NR99 database with SINA v.1.2.11 (Pruesse et al., 2007) (SINA: search-max-result = 7, search-minsim = 0.7). Sequences with an alignment quality score below 90% and those identified as mitochondria or chloroplasts were removed from the ZOTU table. Reads from the four filters of each seawater sample were averaged. Finally, to standardize sequencing depths among the samples, the ZOTU table was read into R (R Core Team, 2018), rarefied 100 times to the lowest read count (n = 31,206) using the R package 1 "vegan" (Oksanen et al., 2019) and averaged to get the working ZOTU table.

Microbial Community Analysis
Microbial composition was assessed by relative read abundance and by relative number of ZOTUs of the major groups of symbionts in each of the 73 samples (i.e., 70 tunic and 3 seawater samples). Core microbial communities, defined as ZOTUs present in all samples, were identified without applying a minimum abundance threshold. Unique ZOTUs of a given locality were defined as those present only in that locality, and in at least two replicates, thus excluding singletons (i.e., ZOTUs found in only one sample). Shannon's diversity and Pielou's evenness indices were calculated for each sample 1 http://www.R-project.org/ using the R package "vegan" and averaged for each locality. Diversity and evenness measures were further compared between localities with a Kruskal-Wallis (K-W) and a Dunn test using the "dunn.test" R package (Dinno, 2017). We used the Benjamini-Hochberg (White et al., 2019) false discovery rate method to correct p-values for multiple comparisons. Rarefaction curves of ZOTU richness at increasing number of reads were obtained with rarecurve function in "vegan." Using the same package, the Bray-Curtis dissimilarity matrix was calculated on rarefied read abundance values with square-root transformation (Legendre and Gallagher, 2001) and used to perform nonmetric multidimensional scaling (nmMDS) analyses using the metaMDS function. PERMANOVA (permutational analysis of variance; Anderson, 2001) tests were done with function adonis of the "vegan" package to compare different groups of samples. When this test was significant, we performed permutational pairwise tests with package "pairwiseAdonis" (Martinez-Arbizu, 2019) to assess differences between pairs of levels of the significant factor, using again the Benjamini-Hochberg method to correct p-values. For the worldwide distributed mitochondrial clade (clade A), the dissimilarity of microbiome samples was compared within localities, between localities for each genetic group (Casso et al., 2019b) and between genetic groups using density plots and analyzed with a K-W and a Dunn test with the "dunn.test" R package as described above.

Comparison With Other Variables
For the clade A, we tested the correlation between the ordination scores of the nmMDS configuration and the temperature variables in each locality (minimum, maximum, mean, and range of temperatures), obtained from World Sea Temperatures 2 , using the envfit test in "vegan." We further explored in the clade A dataset the role of the hosts' genotypes and environmental parameters (thermal and geographic distances between localities) on the microbiome composition. Prevosti genetic distances between colonies were obtained from a previous study (Casso et al., 2019b), in which exactly the same colonies were analyzed using multiple loci recovered with a genotyping-by-sequencing procedure. Thermal distances were calculated as differences in the temperature variable that best explained the ordination found in the previous analysis (envfit). Geographic distances were calculated as the shortest distance by sea (excluding routes through the Panama Canal) obtained from SEA-DISTANCES.ORG 3 . We performed a stepwise multiple regression analysis to model microbiome structure (Bray-Curtis distance) as a function of genetic, thermal, and geographic distances. The analysis was run with the "MASS" R package (Venables and Ripley, 2002) in forward mode (adding variables sequentially) using the Akaike Information Criterion (AIC) to evaluate the fitness gain as variables are added to the model. A global analysis was performed for the complete dataset and two restricted analyses were done for EUR (Europe) and WNA (West North America) as these areas were colonized by a single ascidian genetic group (Casso et al., 2019b), allowing comparisons within geographically restricted and genetically similar populations.

Microbial Community Composition and Diversity
On average, 100,326 reads per sample were retained in the ZOTU table, with a minimum of 31,206 (subsample size for rarefaction) and a maximum of 226,123 reads. A total of 4,323 ZOTUs were identified among the 73 samples, with an average of 266.4 ZOTUs per sample (SE ± 30.78). Supplementary  Table S1 lists the ZOTUs retained, with their taxonomic affiliation, abundance per sample, and sequence. Except for seawater samples, rarefaction curves reached clear asymptotes, indicating that sequencing depth was sufficient to capture microbiome diversity per sample (Supplementary Figure S1). Most of the ZOTUs (4,236) corresponded to 36 different bacteria phyla and the remaining ones (87) to 3 Archaea phyla, which represented less than 0.1% of the reads. Among bacteria, the most abundant phylum in number of reads was Proteobacteria (79.08% of the total reads of bacteria), with Alphaproteobacteria representing 41.79% of bacterial reads, and Gammaproteobacteria, 36.77% (Figure 2). Other dominant phyla were Spirochaetes (11.25% of bacterial reads), Bacteroidetes (5.24%), Epsilonbacteraeota (1.08%), and Cyanobacteria (0.89%). The remaining phyla combined accounted for ca. 2.5% of the reads (Figure 2). One sample from Sugashima clade A showed a visibly different composition (marked with an asterisk in Figure 2), which may indicate a contamination or some artifact during processing. In terms of number of ZOTUs (Supplementary Figure S2), Proteobacteria was again the most diverse (i.e., ZOTU-rich) phylum, with Gammaproteobacteria accounting for 1,183 ZOTUs (27.89%) and Alphaproteobacteria totaling 683 ZOTUs (16.12%). The second most diverse phylum was Bacteroidetes (846 ZOTUs, 19.97%). Spirochaetes, despite their abundance in some samples, had a low diversity (21 ZOTUs, 0.5%).
When analyzing the 64 samples from D. vexillum (excluding the sample from Sugashima clade A with a different composition), we found 3,525 ZOTUs (mean = 210.40; SE ± 14.91), and the taxonomic composition was similar to that of the global dataset with the 73 samples. The ZOTUs corresponded to the same 36 bacteria phyla (3,491 ZOTUs) and 3 Archaea phyla (34 ZOTUs) as in the global dataset. Archaea phyla also represented less than 0.1% of the reads, and the same most abundant bacteria phyla were identified using both datasets (Supplementary Table S1).
Seawater samples had the highest number of ZOTUs (1,738), unique ZOTUs (31.19%), and diversity values (H = 4.810; J = 0.664) ( Table 2). Among D. vexillum samples, the maximum numbers of ZOTUs and unique ZOTUs were found in Sugashima clade B samples (n = 1,381 and 1.74% unique) and Brest (n = 1,028, of which 4.47% were unique) ( Table 2). The highest values of Shannon and Pielou Indices were found in Venice (H = 2.946; J = 0.578) ( Table 2). In comparison to D. vexillum, Didemnum sp. samples from Roscoff showed a high number of ZOTUs (n = 1,025, of which 1.85% were unique) and the lowest Shannon and Pielou Indices (H = 1.884; J = 0.329). The K-W analyses of diversity (Shannon) and evenness (Pielou) showed significant differences between sample groups (K-W χ 2 , 38.942 and 44.190,  We list the total number of ZOTUs, number (and percentage) of unique ZOTUS (present only in a given sample group but in at least two samples), and average Shannon and Pielou Indices (with standard errors within parentheses). Sample groups correspond to each locality, with Ebro Delta tunic and seawater samples and clades A and B samples from Sugashima separated. Note that each sample group has five replicates except for water, with three replicates, and Sugashima clade A, with four replicates (one sample was eliminated due to its different composition, see text).

Microbial Community Structure
The global nmMDS (Figure 3) showed differences between D. vexillum, Didemnum sp., and seawater samples. The sample from Sugashima clade A with a different composition (marked with an asterisk in Figure 2) also appeared as an outlier in the ordination results (Figure 3) and was therefore removed from further analyses. Samples from clade B showed differentiation from clade A samples (albeit with some overlap), and some substructure was detected within clade A (Figure 3). Numbers are percentages of the total number of reads for each group of samples. ZOTU numbers as in Supplementary Table S1.
The PERMANOVA analyses comparing Sugashima clade B, seawater and clade A samples showed a highly significant effect (Supplementary Table S3), and all pairwise comparisons proved significant. Seawater samples (n = 3) were compared with tunic samples (n = 5) from the same locality (Ebro Delta). Among the 832 ZOTUs present in tunic samples, 526 were also found in seawater samples. The 10 most abundant ZOTUs from tunic samples represented the 90.79% of the total reads of the tunic samples and only 1.39% of the total reads of the seawater samples. Conversely, the 10 most abundant ZOTUs from seawater samples represented the 37.51% of the total reads of the seawater samples and only 1.22% of the total reads of the tunic samples. All this pointed to clearly differentiated microbiome structure between ascidians and surrounding seawater.
To further analyze the differences between D. vexillum clades and within clade A, we decided to analyze separately the samples from Japan (the area where both clades were present). Among the 14 samples from Japan, 15 ZOTUs were shared by all samples, which represented 77.3% of the reads from these samples. The same 15 ZOTUs were shared by the nine Aomori and Sugashima clade A individuals (making up for 83.57% of the reads of these individuals), 17 ZOTUs (77.02% of the reads) were shared by Aomori and Sugashima clade B, and 25 ZOTUs were shared by Sugashima clades A and B (comprising 84.70% of the reads of these samples). The Japan nmMDS (Figure 4) clearly differentiated samples from Aomori, Sugashima clade A and Sugashima clade B. The PERMANOVA analysis (Supplementary Table S3) showed overall significant differences among the three groups, and all pairwise comparisons were also significant.
Among the 59 microbiome samples belonging to clade A, we found a total of 3,154 ZOTUs. The core community was constituted by eight ZOTUs shared by all samples, which represented 69.57% of the reads, but only 0.25% of the total number of ZOTUs from these samples. The eight core ZOTUs were all bacteria: four Alphaproteobacteria, three Gammaproteobacteria, and one Spirochaetes ( Table 3). These ZOTUs were also present in all clade B and Didemnum sp. specimens, but with decreasing abundance as genetic distance increased. Thus, these ZOTUs accounted for 50.89% of the reads of clade B specimens and only 20.10% of the reads in Didemnum sp. samples (Table 3). Finally, the eight core ZOTUs were also present in the seawater samples, where they represented only 1.39% of the reads, implying they were enriched ca. 50-fold in the clade A ascidian tunics ( Table 3).
Using the clade A dataset, pairwise Bray-Curtis distances between samples from the same locality averaged 0.512 (SE ± 0.009), between samples from different localities in the same genetic group (i.e., within Japan, WNA, NZ-ENA or EUR, as defined in Casso et al. (2019b) note that these groups coincide with geographic regions except in the case of NZ-ENA), they averaged 0.650 (SE ± 0.004), and between samples from different genetic groups, the mean distance was 0.696 (SE ± 0.002) (Figure 5). These mean distances were significantly different from each other (K-W χ 2 = 291.81, df = 2, p < 0.001, followed by pairwise Dunn tests, all p < 0.001).
The clade A microbiome nmMDS (Figure 6) showed a general pattern of sample grouping by locality, although the degree of dispersion is variable. Atlantic samples (Portsmouth and Brest) grouped together; San Francisco Bay (Richmond and Sausalito) and Seattle samples also grouped together. According to the main axis, Japanese samples (Sugashima and Aomori) grouped together. Although with more dispersion, the main axis also grouped together both Mediterranean samples (Delta and Venice), Woods Hole and Alaska (Sitka). New Zealand samples (from Nelson) appear more dispersed. The PERMANOVA analysis (Supplementary Table S3) demonstrated significant differentiation, and the pairwise tests were all significant after FDR correction except for the comparison between Richmond and Sausalito.
For the restricted datasets analyzed, in the case of Europe (Supplementary Figure S3), the nmMDS showed a clear separation of the Atlantic and Mediterranean localities, and also between localities in each sea. Accordingly, the PERMANOVA results (Supplementary Table S3) showed a significant effect of locality, and all pairwise comparisons were significant. For the North West Atlantic (NWA) dataset (Supplementary Figure S4), the nmMDS configuration showed a sharp differentiation between the two San Francisco Bay localities and those further North (Seattle and Sitka). This is confirmed in the PERMANOVA analysis, which detected overall significant differentiation. The pairwise tests showed that only the comparison of Sausalito and Richmond was not significant.

Comparison With Other Variables
Using the clade A samples (Figure 6), the environmental temperature variables (minimum, maximum, mean, and range; Supplementary Table S1) were fitted to the nmMDS ordination scores. All the relationships were significant, but the better fit was with temperature range (T min r 2 = 0.372; T max r 2 = 0.490; T mean r 2 = 0.219; T range r 2 = 0.731; all p < 0.01). The corresponding vector is plotted in Figure 6. Temperature range differences were therefore chosen to calculate thermal distances used in the following analysis.
According to the forward stepwise regression procedure, temperature range was the best variable explaining the microbiome structure (as ascertained by the AIC minimization), the second variable added was the genetic distance, and in the third place, geographic distance. The three variables were retained in the final model, and the fit was highly significant (adjusted r 2 = 0.280, p < 0.001, Supplementary Table S4).
As for the restricted datasets, in the case of Europe (Supplementary Figure S3), the envfit analysis showed that the maximal temperature provided the best fit to the ordination scores (T max , r 2 = 0.934, vector depicted in Supplementary Figure S3). The multiple regression analysis showed a best fit model including first temperature (as T max ) and, next, genetic distances, but excluding geographic distances (final model, adjusted r 2 = 0.634, p < 0.001, Supplementary Table S4).
In the NWA localities (Supplementary Figure S4), the temperature variable best correlated with the ordination scores was the mean temperature (T mean , r 2 = 0.910, vector depicted in Supplementary Figure S4). The multiple regression analysis revealed again that the best final model included temperature (using T mean ), then genetic distance as the second explanatory variable, and geographic distance was excluded from the final model (adjusted r 2 = 0.679, p < 0.001, Supplementary Table S4). Thus, contrary to what we found in the global analysis, in the spatially restricted datasets, geographic distance effects did not explain the microbiome composition once the other variables have been accounted for.

DISCUSSION
In this study, we provide the first characterization of the tunic microbiome of the invasive colonial ascidian D. vexillum. We were able to analyze the microbiome of both mitochondrial clades of this species, from exactly the same samples genetically genotyped in a previous study that investigated the population genomics of the host (Casso et al., 2019b). These samples were collected from representative localities of the whole range of distribution of D. vexillum, providing an exceptional opportunity to study the differences in the microbiota composition of the different populations worldwide and their relationship with the host genetic relatedness. We also analyzed the microbiome of a population from a closely related species (Didemnum sp.) plus seawater samples from one of the localities, for comparison.
We found that the number of ZOTUs found among D. vexillum samples (3,525 ZOTUs) represents a highly diverse microbiota comparable to that found in other ascidians and marine invertebrates (Erwin et al., 2014;Cleary et al., 2019). The large difference between abundance and diversity of Archaea and bacteria (the former representing ca. 2% of ZOTUs and 0.1% of reads) in our results is likely to be confounded by primer amplification biases downplaying the representation of Archaea (Klindworth et al., 2013;Turon and Uriz, 2020).  Given the encrusting nature of the colonies and the existence of common colonial cavities in contact with the external environment in D. vexillum, as in other colonial ascidians, it is not possible to analyze only internal tunic. In our samples, we avoided basal surfaces in contact with the substrate, but had to include distal surfaces with cuticle where seawater bacteria could be found (López-Legentil et al., 2016). However, as the fixative for the sampling of the colonies was changed several times, and colonies were rinsed during manipulation, we are confident that any remaining non-internal bacteria were truly attached to the ascidian surface. Little difference was found between the microbiome of the inner tunic of a solitary ascidian and that contained in the outer tunic including the cuticle surfaces (Blasiak et al., 2014). Likewise, similar bacterial composition was detected between tunic and tunic with cuticle in several ascidian species (Cahill et al., 2016), indicating that few prokaryotic groups occur exclusively in the exterior of the tunic.
Seawater samples showed the most diverse bacterial communities, with 1,738 ZOTUs, which is 2.5 times greater than the average number of ZOTUs in the ascidian tunics (679, Table 2). This ratio may be even higher because we had only three seawater replicates, and the sequencing depth was insufficient to capture adequately the microbiome diversity in these samples, as the rarefaction curve did not reach a clear asymptote (Supplementary Figure S1). The difference in diversity we found between seawater and tunic samples is consistent with other studies comparing seawater and microbiomes from other ascidians (Erwin et al., 2014). There was also a marked difference in composition between seawater and tunic samples, as shown in the global nmMDS analysis (Figure 3). A similar degree of differentiation was also observed between D. vexillum and the congeneric Didemnum sp. samples. This indicates that the microbiome composition is clearly different between species, even if they are closely related and living in the same kind of environment, a result found in previous studies in ascidians (Erwin et al., 2014;Cahill et al., 2016) and other marine invertebrates such as sponges, chitons, corals, or sea cucumbers Cleary et al., 2019). The reduction in abundance of the clade A core microbial community when increasing the host genetic distance (Table 3) also points to the evolution of species specific microbial communities. Samples from clade B individuals showed some differences in composition with respect to clade A specimens from the same locality (Sugashima) (Figures 3, 4). These differences may be explained by the possible speciation process occurring between D. vexillum clades A and B (Casso et al., 2019b).
The core community of D. vexillum clade A samples was highly abundant (69.57% of reads) but with low diversity (eight ZOTUs). On the contrary, the variable component has low abundance and high diversity. The microbiome is related to changes in the environment (such as temperature variables and geographic distance between localities in the global analysis), which points to horizontal acquisition of the symbionts from the seawater. However, there seems to be a selective enrichment in the host of particular bacteria. For instance, the 10 most abundant ZOTUs from tunic samples in the Ebro Delta represented the 90.79% of the total reads of the tunic samples and only 1.39% of the total reads of the seawater samples, implying that there is a considerable difference in abundance in the ZOTUs present in both tunic and seawater samples. This further suggests the role of the environment on the horizontal transfer, but with an enrichment promoted by the host. A small but abundant core and a highly diverse variable microbiome component, with marked selective enrichment from the pool of water prokaryotes, has also been found in sponges  and may contribute to the host resilience to environmental changes (Leite et al., 2018;Turon et al., 2019).
We found a relationship between microbiome structure and the genetic makeup of the colonies, both in the global analysis and in the restricted analyses of particular genetic groups (Europe and North Western America, where genetic distances, and not geographic distances, explained microbiome composition). These results indicate either a vertical transmission of some of the symbionts, or the heritability of traits related to the selection of symbionts from the seawater. We cannot, at present, resolve these factors unambiguously, but vertical transmission seems the most likely one. The presence of prokaryotes in the embryos and larvae of colonial ascidians has been noted (Moss et al., 2003;Martínez-García et al., 2007;López-Legentil et al., 2011), and transmission of symbionts via larvae has been well documented in photosymbiotic relationships between ascidians and the cyanobacterial Prochloron species (Hirose, 2015;Hirose and Nozawa, 2020). Moreover, the larvae of other species of Didemnum are known to harbor bacteria (López-Legentil et al., 2015), making vertical transmission a highly plausible trait in our case. Vertical transmission is not necessarily fixed and consistent (Björk et al., 2019), but can contribute to propagate particular symbionts. We suggest that a combination of horizontal and vertical transmission shapes the microbiome of D. vexillum, as suggested for other invertebrates such as sponges (Sipkema et al., 2015).
It is noteworthy that all temperature variables and, particularly, the temperature range in the global analysis correlate with the microbiome structure of D. vexillum. In the restricted analyses of Europe and North Western America, temperature range differences were not so large, but nevertheless maximal and mean temperature, respectively, were closely correlated with the ordination of the samples, again confirming the role of temperature in shaping the microbial community. Didemnum vexillum is defined as a temperate species, however, it can live from subtropical seas such as the Mediterranean to cold regions such as Alaska. Hibernation and estivation processes may help the species to cope with unfavorable temperatures (Bullard et al., 2007;Ordóñez et al., 2015). Our results indicate that temperature adaptation seems to be a driving force in the colonization process of D. vexillum, and the microbiome can play a role in the thermal adaptation of the species, as has been suggested in fish and corals (Kokou et al., 2018;Osman et al., 2020).
In the marine realm, many instances of successful introduction of species outside their native ranges have occurred (Carlton and Geller, 1993), which implies fast adaptation to new conditions (Rius et al., 2015). It has been pointed out that the microbial symbiont community can contribute to the survival and adaptation of the holobiont during the invasion process (Aires et al., 2016;Amsellem et al., 2017;Arnaud-Haond et al., 2017). Among other mechanisms, thermal adaptation and the production of active substances may be key for the introduction success and can be mediated by the microbiome (Dou and Dong, 2019;Osman et al., 2020). Didemnum vexillum must be able to adapt to a wide range of environments in its worldwide introduction range. Chimerism (Casso et al., 2019a) and epigenetic changes (Hawes et al., 2019) have been suggested as potential mechanisms enhancing its adaptive capacity, and the microbiome is likely to play a significant role as well. We found that the microbiome composition of the introduced clade of D. vexillum is significantly influenced by temperature, genetics of the host, and geographic distance. The combination of a quantitatively dominant core microbiome with a highly diverse variable component, likely acquired by a combination of horizontal and vertical transmission, can contribute to the success of this global invader in different environments.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the Short Read Archive (SRA) of the NCBI under Bioproject PRJNA608680.

AUTHOR CONTRIBUTIONS
XT and MP designed the research. MC, MP, and XT collected samples. MC performed laboratory work and wrote the first draft of the manuscript. MC, MT, and NM ran the bioinformatics pipeline, analyzed results, and prepared tables and figures. All authors discussed and interpreted results, and revised the manuscript.