Samples were collected from six different locations around the island of Curaçao (Caribbean Sea) (Figure 1) in January 2017 with an average seawater temperature of 28–32°C. Three seagrass species were sampled: the native T. testudinum and H. wrightii and the invasive H. stipulacea. At each location, three replicate specimens of the species present in each site were collected along with environmental samples (i.e., sediment and seawater), resulting in a total of 112 samples (see the description in Figure 1). All seagrasses were collected in the shallow subtidal between depths of centimeters to a few decimeters below sea surface with a few meters distance between samples. The different species were mostly collected from monospecific meadows, though different species were always sampled close to each other with distances between meadows within 10 m. Only at the entrance of Spanish Water did the species occur in a sympatric meadow. Individual specimens were uprooted and washed free of sediment, epiphytes, and anything loosely attached using seawater from the sampling location. Excess seawater was shaken off and leaves and roots were separated and directly preserved in Xpedition^TM Lysis/Stabilization Solution (Zymo Research, California) and preserved at −20°C until further extraction. Collected sediment samples had a volume of roughly 1 ml and represented the sediment at the root depth of the respective seagrass. The portion of sediment, for each sample, was scooped using sterilized 2.5-ml tubes. Seawater bacterial communities were collected by filtering 0.5 L of seawater from the water column above the seagrass meadow, using 0.2-μm filters. The filters were preserved similarly to seagrasses and sediment samples. Environmental samples were collected as controls (1) to compare their bacterial community with those from the seagrasses and (2) to discard environmentally related operational taxonomic units (OTUs) when determining seagrass core OTUs.

Total genomic DNA was extracted from all samples using the Quick–gDNA kit (Zymo Research^TM) according to the manufacturer’s protocol for “Solid Tissue Samples” (page 4 of the manual). For all the samples, including sediment (0.25 g) and seawater filters, in the lysis step, we have used tungsten beads and an automatic homogenizer (for 2 min, at a frequency of 20/s) for a more efficient mechanic lysis. For bacterial community characterization, the nearly complete 16S rRNA genes were amplified using the universal primers 27F and 1492R with the following changes to the original protocol (Lane, 1991): initial denaturation at 95°C for 2 min, 35 cycles of denaturation at 95°C for 20 s, annealing at 55°C for 20 s, and extension at 72°C for 90 s, with a final extension at 72°C for 3 min. The 25-μl reaction mixture contained 250 μM dNTPs, 0.6 μM of each primer, 1 × 2 PCR buffer mix, 2 μl of template DNA (samples were diluted to a final concentration of about 10 ng μl^–1 per reaction), and 0.3 μl of Taq polymerase (Advantage^® R2 Clontech^TM). PCR products were cleaned using the ExoFastAP enzyme following the protocol of the manufacturer Thermo Scientific^TM. Amplified DNA was sent to Molecular Research (MR DNA, Shallowater, TX, United States) where a nested PCR was performed prior to sequencing. The modified 8-bp key-tagged primer 799F along with the reverse primer 1193r, covering the regions V5–V7 of the 16S rRNA and amplifying a fragment of ∼400 bp, were used to avoid chloroplast cross-amplification (Bodenhausen et al., 2013). PCR conditions were as follows: 95°C for 3 min, 10 cycles of 95°C for 20 s, 50°C for 30 s, 72°C for 30 s, and a final elongation of 72°C for 3 min. The samples were pooled in equal proportions based on their molecular weight (calculated based on the size of the amplicon) and DNA concentrations (using Qubit Invitrogen^®) and purified using calibrated Agencourt^® AMPure^® XP beads. DNA libraries were prepared by following Illumina TruSeq DNA library preparation protocol and paired-end (2 bp × 250 bp) sequencing performed at MR DNA¹ (Shallowater, TX, United States) on an Illumina MiSeq system following the manufacturer’s guidelines.

The Q25 sequence data derived from the sequencing were processed using the MR DNA ribosomal and functional gene analysis pipeline (see Text Footnote 1, MR DNA, Shallowater, TX, United States). Sequences were depleted of primers; short sequences <150 bp and sequences with ambiguous base calls were removed. Sequences were quality filtered using a maximum expected error threshold of 1.0 and dereplicated. The dereplicated or unique sequences were denoised; unique sequences identified with sequencing or PCR point errors were removed, followed by chimera removal. The demultiplexing step was performed using the “split_libraries.py” Qiime command and default parameters. Subsequent data processing was done using QIIME version 1.9. (Caporaso et al., 2010) and clustered into OTUs at >97% similarity using open-reference OTU picking with the UCLUST algorithm against the SILVA reference database (version 132²). After selection of one representative sequence per OTU, the database was aligned. Taxonomy was assigned using the SILVA database (version 132). From the resulting OTU table, eukaryotic organelle sequences (i.e., chloroplasts and mitochondria) and unassigned sequences were removed, as well as rare OTUs (singletons and doubletons). The table used for downstream analysis is from now on referred to as OTU table.

One of the samples (root from H. wrightii sampled in Ascension) resulted in a very low number of sequences (7,217), compared to the remaining 111 samples (median: 65,124), likely due to DNA extraction, PCR, or sequencing problems. Hence, we decided to omit it from our final dataset. For statistical purposes, samples were normalized by rarefaction to the minimum number of sequences (29,021), per sample, and used in alpha and beta diversity analysis. Rarefaction was performed to adjust for differences in library sizes across samples to aid comparisons of alpha and beta diversity. The rarefaction method chosen (rarefaction to the minimum number of sequences) allows comparability with previous studies (using that same method). For all the other compositional analysis (differential analyses and core community), non-rarefied data were used so we would have the real proportion of the OTUs/taxa analyzed without the “randomly selected” factor of the rarefaction method.

In order to predict the functional profile of the seagrass-associated bacterial communities, a taxonomic profile of KEGG Orthology (KO) pathways was obtained from our SILVA-based OTUs (version 132) using the Tax4Fun2 library in R (Wemheuer et al., 2020). Within the Tax4Fun2 library, the functional prediction was obtained using the runRefBlast function with the database mode set to “Ref100NR” and the path_to_otus argument set to the OTU representative sequences fasta file. The function performed a next neighbor search against the 16S rRNA reference database using the BLAST algorithm. Subsequently, the makeFunctionalPrediction function was run. This function summarized the results of the next neighbor search performed during the runRefBlast in a relative abundance table of KO pathways per sample. The path_to_otu_table argument was set to the OTU table. The min_identity_to_reference argument was set to 0.95, meaning that OTUs passing a similarity threshold of 95% to their next neighbor were considered for further analysis. The specified cutoff level results in predicted function for a given OTU at an approximately generic level. This level was chosen as metabolic capabilities among prokaryotes have shown to vary also below family level (Goldford et al., 2018).

The calculated diversity indices (Observed OTUs and exponentiated Shannon diversity—the order in which the p-values are presented below) showed differences among seagrass species (F_2,73 = 6.904, p = 0.002/F_2,73 = 4.987, 0.003) and tissues (F_1,73 = 92.487, p = 0.001/F_1,73 = 49.072, 0.004). Only the Observed OTUs index showed differences among locations, which depended on the seagrass species and tissue (Species × Location F_5,73 = 2.22, p = 0.034; Tissue × Location F_5,73 = 2.985, p = 0.017).

Differences among species depend on the tissue analyzed for both indexes (Species × Tissue F_2,73 = 20.933, p = 0.001/F_2,73 = 25.318, 0.001). Looking at the observed OTUs index, the invasive H. stipulacea only shows a higher diversity when comparing their root community with that of T. testudinum (Supplementary Table 5 and Figure 2). However, the exponentiated Shannon results show that the invasive roots’ community diversity is higher than both native species (Supplementary Table 13 and Figure 2), while for the leaves, H. wrightii bacterial diversity is higher than that of H. stipulacea (Supplementary Table 13 and Figure 2). All the p-values for the α-diversity statistical analysis are in Supplementary Tables 1–13.

Multivariate analysis disclosed a clear separation of bacterial communities for the different factors (species, tissue, and location) and their interactions (Table 1 and Figure 3). Environmental bacterial community structure from sediment and seawater differed strongly from each other (p = 0.001, Figure 3) as well as from the seagrass-associated bacterial communities (p = 0.001, Figure 3). Within the seagrasses, there was a strong composition differentiation between roots and leaves (Table 1 and Figure 3), explaining about 11% of the total variation. In addition, the bacterial community structure differed among the three seagrass species (Table 1, Figure 3, and Supplementary Table 14), explaining close to 7% of the total variation. The interaction of Tissue and Species explained another 11% of the variation. So, seagrass species and tissue formed the most important structuring factors of the seagrass microbiome composition. As a random factor, location contributed little to the total variation (∼3%). Pairwise tests showed differentiation among all locations within all three seagrasses (p = [0.001–0.027]).

Differentiation among microbiomes of all three seagrass species was detected in leaves and roots (Supplementary Table 17) and within every sampled bay (Supplementary Table 18). All the p-values for the β-diversity statistical analysis that are not represented in Table 1 are in Supplementary Tables 14–24.

In the leaves, the invasive H. stipulacea microbiome differentiation from the two native seagrasses correlated strongly positive with the read abundance of Alteromonadales (Figure 4A, vectors 9) and specific OTUs from the Rhodobacteraceae family (Figure 4A, vectors 2). The family Rhodobacteraceae had specific members in the leaves of each of the three seagrasses, but especially in T. testudinum. T. testudinum microbiome structure also showed a distinguished correlation with members of the Rhizobiales (Figure 4A). The bacterial groups driving root differentiation among seagrass species were very different and more diverse. Members of, for example, the families Rhodobacteraceae, Piscirickettsiaceae, and Oceanospirillaceae and the orders Chromatiales and Bacteroidales positively correlated with the invasive H. stipulacea and members of the Caldithrix are more associated with native seagrasses, especially T. testudinum (Figure 4B). Some members of Desulfobacteraceae were rather specific for the invasive seagrass root while others were so for the native seagrass roots (Figure 4B). Native seagrass roots showed a higher number of discriminating taxa vectors of sulfate-reducing bacteria (e.g., Desulfosarcina, Desulfobacteraceae, etc.) (Figure 4B).

Using Tax4Fun2, a table with the predicted gene copy numbers was created displaying a total of 365 KEGG Orthologies and their associated level 1, level 2, and level 3 pathways. The total amount of OTUs used in the prediction was 17,887, whereas the proportion of sequences not used in the functional prediction varied among samples from 0.35 to 0.91 (Supplementary Tables 29–31).

The bacterial community associated with the roots of the invasive seagrass, H. stipulacea, was highly diverse (for both indexes used), and that diversity was higher when compared to those of the native species. During the process of introduction and invasion, introduced individuals usually go through a bottleneck that is thought to strongly reduce the genetic diversity of the introduced macro-organism and a stripping of associated organisms (Rocca et al., 2019). However, other studies have also reported a higher diversity of rhizosphere microbial communities of several terrestrial plants in the invasive compared to their native range (Coats and Rumpho, 2014; Dawson and Schrama, 2016; van der Putten et al., 2016; Lu-Irving et al., 2019). Given the higher diversity of H. stipulacea compared to the two native species, we could assume that the invasive species (i) may profit from that habitat shift acquiring a more diverse bacterial community that would give it an advantage for a successful establishment and proliferation (e.g., by improved nutrient acquisition) or (ii) is naturally highly diverse and there was actually a decrease from its native range that was not enough to go as low as the diversity of the native species. These assumptions could only be tested by sampling H. stipulacea in its native habitat and comparing its α-diversity to that of in its invaded habitat. Previous results on bacterial communities associated with roots and leaves of H. stipulacea in the Red Sea (Mejia et al., 2016; Rotini et al., 2017) have shown a much lower α-diversity (Shannon) for leaves (4.05 ± 0.21–5.04 ± 0.15) and roots (4.49 ± 0.02–5.30 ± 0.46) when compared to our results (8.97 ± 0.64 and 10.35 ± 1.20, leaves and roots, respectively). This gives some strength to our first assumption. However, these authors have used 454 sequencing technology with sequencing depths much lower than those accessed by Illumina sequencing, as well as different primers, both compromising the direct comparison between ours and these two studies.

Increasing the knowledge on the different factors shaping the bacterial community of seagrasses and unveiling some key taxa associated to them would contribute to understand the influence of the microbiome on the invasion process, especially when co-analyzing invasive and native species. Bacterial community structure was seagrass species-specific for both root and leaf communities. So far, there are only a few studies that compare the bacterial communities of different seagrass species, and their results are contradictory. For example, Roth-schulze and Thomas (2016) studied the diversity of microbial communities on different types of host surfaces (including the seagrass species Posidonia australis and H. ovalis) and found that the bacterial community structure was unique to each host. Likewise, Martin et al. (2018) found that the microbial community structure of the rhizosphere of three co-occurring seagrass species, H. ovalis, Halodule uninervis, and Cymodocea serrulata, was unique for each species. In contrast, a recent study on the phyllo- and rhizosphere of the seagrasses T. testudinum and S. filiforme found that the microbial composition of roots and leaves was different without showing species specificity (Ugarelli et al., 2018). At a rhizosphere level, a lack of species specificity for bacterial communities but rather regional differences were reported for Zostera noltii and Zostera marina in Roscoff (France) and Faro (Portugal) (Cúcio et al., 2016). In our case, we detected species specificity at the island and tissue level, but not between co-occurring or neighboring seagrass species within bays. This can be due to our relatively low replication level at this scale. In contrast, seagrass-associated communities differed across all bays within each seagrass species indicating the importance of local effects. Overall, these results show that the composition of seagrass-associated bacterial communities is influenced by different factors such as host species, “microenvironment”/tissue (leaves and roots), and, ultimately, location.

As known from other studies that compared the microbiomes associated with the different plant tissues (e.g., Lebeis, 2015; Ettinger et al., 2017; Rotini et al., 2017; Ugarelli et al., 2018), the separation between leaves and roots was very clear for all the seagrass species in this study, with certain bacterial phylotypes (taxonomic groups) discriminating each of the tissues. This factor (tissue) was the one showing the highest explanation for bacterial community variation among seagrasses (11.2% for the overall differentiation). Representing the different tissues (roots and leaves) in different PCO plots shows the differentiation among the different species better (which was also statistically supported), especially at the leaf level where native seagrass species clearly segregated from the invasive H. stipulacea. For the leaves, the “invasive cluster” has 12 discriminant vectors assigned to the order Alteromonadales. This order was also found to be exclusively associated with H. stipulacea in our core community analysis (more than 50% of those OTUs were assigned to the genus Alteromonas), being present in both leaves and roots but in a higher abundance in the leaves. Alteromonadales was also found to be one of the most abundant bacterial orders associated with H. stipulacea from the northern part of the Red Sea in two independent studies (Mejia et al., 2016; Rotini et al., 2017). Also, Weidner et al. (2000) isolated a bacterium closely related to Alteromonas macleodii from H. stipulacea leaves, which is responsible for producing PGP oligosaccharides (Ferrier and Martin, 2002). These results provide support to the possible important and structural role that this phylotype may have for this seagrass species. Besides, some Alteromonas species isolated from algae have been described as having enzymatic activities relevant to the degradation of macroalgal cell walls (Goecke et al., 2010). Having those bacteria associated with their leaves could bring an advantage to this invasive species by helping it to get rid of the algal epiphytes and mitigate their possible competitive effect for light and gas exchanges that leads to reduced plant fitness (Brodersen et al., 2015). In addition, Alteromonadales have been suggested to collaborate with Thiotrichales to degrade marine dissolved organic matter in seawater (Mccarren et al., 2010). Our results on the core bacterial community analysis, specific to each species, show that Thiotrichales also show up exclusively associated with H. stipulacea, mainly and most abundantly with leaves, which could be a clue for that concerted work between these bacteria to bring some advantage to the invasive species. Experiments in axenic cultures of this seagrass species with subsequent inoculation of these two phylotypes could help clarify the correlation between these bacteria and its positive effect in H. stipulacea. However, obtaining axenic seagrasses is a challenge by itself.

In the native species, the leaf microbiome also forms statistically different clusters but with some level of mixing among some T. testudinum samples with the H. wrightii cluster. In that mixed cluster, we can find several discriminating vectors corresponding to the different differential bacteria, showing a lower level of specificity than that of the invasive species. The cluster that gathers most of the Thalassia samples shows the order Rhizobiales as a differential vector of this species that are known to fix nitrogen. Bacteria from the order Rhizobiales were consistently found in high relative abundances in root and leaf communities of phylogenetically diverse plant hosts and hence were considered part of the core plant microbiota (Garrido-Oter et al., 2018). The only vector that was common to all the seagrass species belonged to the family Rhodobacteraceae. Regardless, some of the OTUs assigned to the Rhizobiales order were found exclusively associated with Thalassia samples; our core analysis also shows this order as one of the most abundant orders common to all the three species. In other studies, these bacteria were found consistently and abundantly associated with leaves of several seagrass species as H. stipulacea (Mejia et al., 2016; Rotini et al., 2017), Thalassia hemprichii (Ling et al., 2015), T. testudinum, S. filiforme (Ugarelli et al., 2018), and Z. marina (Ettinger et al., 2017). These bacteria have been described as “dominant and ubiquitous primary surface colonizers in coastal waters” (Dang et al., 2008). Some Rhodobacterales groups are known to produce antibacterial compounds that act as a deterrent for other bacteria (Dang et al., 2008; Ettinger et al., 2017), which might explain their persistence and adaptation to different seaweed species and environmental conditions (identified as a major coastal biofilm epiphyte in the Atlantic and Pacific Oceans; Dang and Lovell, 2000).

As for the roots, although statistics show significantly different bacterial communities among the different species, there are only two different groups of differentiator vectors discriminating (pointing toward) the native vs. the invasive species. Most of the discriminating vectors associated with the roots of the native seagrass species were assigned to anaerobic sulfate reducer (SR) bacteria, which are known to contribute to the detoxification of the root area (e.g., Desulfosarcina, Desulfococcus, and Desulfarculaceae), and are also capable of carbon decomposition (Kleindienst et al., 2014). The invasive Halophila root differentiation vectors are mostly organic carbon (e.g., Oceanospirillales, Myxococcales, and Bacteroidales) and hydrocarbon (e.g., Oleibacter) degraders usually found associated with seagrasses roots and rhizosphere (Crump and Koch, 2008; Ettinger et al., 2017; Gribben et al., 2017; Ugarelli et al., 2017, 2018; Cúcio et al., 2018; Martin et al., 2018) and the order Chromatiales was the only SR bacterial vector. The results suggest that, most likely, different bacterial communities (especially the core microbiome associated with native species vs. invasive species) are performing the same functions that have already been shown for the core microbiome of some coral species (Hernandez-Agreda et al., 2018). However, the order Oceanospirillales (in our specific case belonging to the Halomonadaceae family) is also a halotolerant/halophilic taxon (Franchini et al., 2014) and might play a role in H. stipulacea adaptation by conferring some salt tolerance to the host plant (Ruppel et al., 2013). That evidence is strengthened by the fact that the order Oceanospirillales, in our core analysis results, appears exclusively and in high abundances, not only in H. stipulacea roots but also in leaves.

As discussed above, we found a species-specific pattern among the bacterial communities of the three seagrass species. Additionally, we would also expect to find a more similar bacterial community structure between the native species, which evolved together in the same environment, compared to the invasive species that only recently arrived in that environment (Aires et al., 2015). Yet, our results were rather ambiguous; the core community analysis showed a higher similarity between the invasive H. stipulacea and the native H. wrightii sharing 27 OTUs (corresponding to 14 genera) against the native species which only shared seven OTUs (corresponding to 5 genera), while our differential analysis shows that native species have fewer differential genera (97) than invasive vs. native species (H. stipulacea vs. T. testudinum = 312, H. stipulacea vs. H. wrightii = 138). However, it is important to highlight that the depicted Venn diagram was drawn using the core OTU table and not the total OTU table. The core microbiome, besides its ubiquity, has been described as having essential functions, to a certain type of organism community, ranging from nutrition to protection against diseases (Shade and Handelsman, 2012). As such, the core bacteria have also been described as potential symbionts due to their tightness and persistence across different species (Hernandez-Agreda et al., 2018). That could mean that sharing a higher number of that small group of highly persistent OTUs, Halophila and Halodule, may reflect more similar needs (e.g., nutritional, adaptive, and physiological) that could be considered essential to the plants, in contrast with the two natives that might share more environmentally responsive community (Hernandez-Agreda et al., 2018). One could also argue that these two species could have exchanged more bacteria during the invasion process but, in our sampling, Halophila and Halodule were only sampled together in two sites and always together with Thalassia, which makes it less likely to have a higher exchange opportunity than Halophila and Thalassia.

The native and invasive species separated by a large number of differential genera (mostly enriched in the native species, relative to the native) could be explained by the fact that they did not evolve in the same environment, not sharing the same original “bacterial source” which contributed for a higher differentiation. Native species differentiation showed a more conspicuous separation by tissue than by species, which was already discussed above.

Several functional hypotheses of microbial roles can be suggested for the euryhaline capacity of H. stipulacea, a species able to invade different environments mainly due to its tolerance to a wide range of salinities (den Hartog, 1970; Oscar et al., 2018; Por, 1971; Winters et al., 2020). Unlike most of the seagrasses to which the salinity limits range from 20 to 40 PSUs (Touchette, 2007), H. stipulacea thrives in a wide range of salinities that can go up to 60 PSUs (Oscar et al., 2018). The capacity to survive and grow in those high salinities has been associated with its invasive success and advantage over the seagrasses it has been displacing [see review by Winters et al. (2020)]. The interaction between plants and their associated microbiomes play an important role in plant adaptation to new environments and to tolerate stress conditions (Hardoim et al., 2008; Rosenberg and Zilber-Rosenberg, 2011), raising the hypothesis that part of H. stipulacea’s tolerance to salinity could be attributed to its associated microbes (Oscar et al., 2018). Based on the review of Etesami and Beattie (2018), where several halophytic plants were mined for halotolerant PGP bacterial genera, we have selected genera to look for in the three studied seagrass species. The invasive H. stipulacea has shown a much higher relative abundance of most of the considered halotolerant bacteria, with PGP abilities, which may be a clue of the importance of these bacterial phylotypes for this species (Ruppel et al., 2013; Etesami and Beattie, 2018). Several studies have investigated halophyte species with the perspective of seeking for potential sources of halotolerant bacteria with PGP potential (Ramadoss et al., 2013; Sharma et al., 2016). These would be further used as an alternative strategy to, for example, enhance crop salt tolerance while boosting crop growth (Dodd and Pérez-Alfocea, 2012). The genera found to be more abundant in the invasive species have growth-promoting properties that range from N2 fixation, indole-3- acetic acid (IAA) production, siderophore production, phosphate solubilization, and 1-Aminocyclopropane-1-Carboxylate (ACC) deaminase activity. Only some will be discussed next. The genus Halomonas has been widely recognized for its important growth-promoting traits such as (i) the ability to provide phosphorus to the plant under P-limiting conditions; (ii) the production of IAA-like phytohormones, providing Fe to plants through chelation and uptake (siderophores); and (iii) ACC deaminase production, which lowers ethylene (plant stress hormone) levels and promotes plant growth and development under adverse environmental conditions (Etesami and Beattie, 2018; Gupta and Pandey, 2019). Our analysis shows a much higher relative abundance of Halomonas in the invasive species compared to the native ones. This may be an indication that these bacteria could be helping to mitigate the osmotic stress in extreme (or discrepant) salinities and, therefore, contributing to H. stipulacea invasiveness. Furthermore, other highly abundant genera associated with H. stipulacea, Salinimonas, and Salinicola, have been recognized as producing quorum sensing (QS) inhibitors (Romero et al., 2012; Tang et al., 2013). QS is one of the pathogenic bacteria processes that allow to reduce host immune response and the establishment of infection (Reading and Sperandio, 2006). Therefore, inhibition of QS could be increased in Halophila and so alleviating the disease pressure, which will be an advantage compared to the native species and concomitantly promote its successful spreading. As for the halotolerant PGP genus that occurred exclusively in H. stipulacea, Lysinibacillus, besides some of the properties mentioned above, was also found to be cadmium and lead tolerant, removing them from the soil and mitigating the heavy metal stress (Pal and Sengupta, 2019). Rhodococcus and Azospirillum were present in an extremely high proportion in the invasive species compared to the native ones (where their proportion was close to 0); these species were found to greatly increase germination of plants in high salinity conditions (Rueda-Puente, 2017). Despite the strong indications that H. stipulacea could be, indeed, taking advantage from the high abundance of salinity-tolerant PGP, only the isolation of these specific phylotypes with posterior inoculation experiments would give us a solid answer about the potential of these bacteria. We stress that our list of halotolerant bacteria only considered those reported in Etesami and Beattie reviews (Table 1). As such, we did not account for other halotolerant bacteria in our dataset, or other bacteria that could have the same or similar functions (salt tolerant with PGP abilities) and were not documented yet.

From the predicted pathways found, most were not relevant for conferring a possible advantage to the invasive species (considered mostly “housekeeping functions”), and some were plant-specific and not likely found in bacteria. For example, Halophila roots have shown a higher abundance of biosynthesis of isoflavonoids; these are plant secondary metabolites known to play an important role in adaptation to their environment, both as defensive and as chemical signals in symbiotic nitrogen fixation with rhizobia (García-Calderón et al., 2020) and are not produced by bacteria. That means that there is a low representation of those pathways that could possibly bring advantage to the seagrasses, enhancing its importance as possible assets for the invasion facilitation. The putative functional traits that could possibly bring some advantage to the invasive species were differentially found in H. stipulacea roots and were lipoic acid metabolism; glycerine, serine, and threonine metabolism; and hypotaurine and taurine metabolism. Lipoic acid metabolism has been associated with protecting plants against oxidative stress and helping on detoxification (Navari-Izzo et al., 2002). Glycerine, serine, and threonine metabolism pathways have been described as important processes in plants by linking nitrogen and carbon metabolisms and helping to maintain the redox balance and energy levels in stress conditions (Igamberdiev and Kleczkowski, 2018). High levels of hypotaurine have been correlated with high sulfur oxidation and its primary function was assigned to sulfide detoxification in marine invertebrates (Rosenberg et al., 2006); this organic osmolyte may have the same role in plants. Also, it has been shown that taurine promotes the growth of wheat seedling and increased root length (Hao et al., 2004). The leaves of H. stipulacea were also found to have some differentially associated putative advantageous functions such as the biosynthesis of staurosporine that have been implied in the protection against some plant fungal and other pathogens (Park et al., 2006).

All these ideas, as well as all our discussion around the putative functions of some particular bacterial phylotypes, are merely hypotheses, based on literature research and a predictive algorithm (Tax4Fun). These hypotheses are raised as suggestions for future testing by meta-omic analyses and manipulation experiments with bacterial isolates, especially to verify if the associated bacterial community actually confers some adaptive advantage to the invasive H. stipulacea.

This study revealed that bacterial communities are different between the three seagrass species in Curaçao, where the invasive species holds a more diverse bacterial community compared to the native ones, with three times more species-specific OTUs. Bacterial community structure follows a species-specific pattern with some influence of the environment (sampling location), and that pattern holds for both leaves and roots. Despite showing significant differences among seagrasses, the root microbiome is less structured, especially for both native species, which shared a higher number of phylotypes, showing a higher influence of the surrounding sediment, which could facilitate microbiome sharing. For those communities that are tightly connected to the seagrasses, likely performing essential functions (“core community”), there is a higher similarity between the invasive H. stipulacea and the native H. wrightii, which could mean that these two species have a more similar organizational and physiological structuring (similar essential functions) relying on similar core-associated bacteria. The bacterial taxa that are, differentially, more closely associated with the invasive species, along with the highly abundant halotolerant/halophilic bacteria, have been described in the literature to possibly perform beneficial functions, such as growth promotion, epiphyte degradation, salinity tolerance, detoxification, nutrient provision, and QS inhibition, which could help to prevent diseases. Besides, predicted functional profiles have also revealed that the invasive species has a higher number of differentiating putative KEGG pathways, especially at the root level, as compared to the native species, from which only a few could be considered advantageous to the invasive species by bringing possible protection against pathogens, helping in detoxification and increasing stress tolerance. All these are predictions and should not be taken as certain but instead used to formulate a new hypothesis and open paths to a more detailed experimental investigation about seagrass microbiome function and their role in the invasion process. However, it is important to keep in mind that these differences may also be related to physiological and morphological differences among seagrass species. In order to assess this hypothesis, further analysis using meta-omics analyses should be done, linked to experimental manipulations that could provide evidence on the role of these functional differences. In the case of H. stipulacea, we suggest a future comparison of the bacterial community from its native range with that from both the Mediterranean and Caribbean habitats where this species has a contrasting invasion success. Nevertheless, this investigation provides novel insights, information, and clues on native vs. invasive seagrass microbiomes that fill in this knowledge gap and guide more in-depth investigations on how associated microbiomes might influence the invasion process.