Metabarcoding of the bacterial assemblages associated with Toxopneustes roseus in the Mexican Central Pacific

1. Introduction

Marine microorganisms are central drivers of ocean biogeochemical cycles; they control the emission of radioactive gases and are essential for community functioning because they provide ecosystem services [1,2]. The diversity of microbial metabolisms has allowed bacteria to succeed as free-living organisms and as symbionts [3]. These are relevant in multicellular organisms considered “holobionts” because they are composed of the macroscopic host and its synergistic interdependence with bacteria, archaea, fungi, and other microorganisms [4,5].

The associations between marine invertebrates and microbial assemblages have been studied in recent years in terms of diversity, composition, and stability [6,7]. Specifically, the bacterial assemblage of marine invertebrates such as sponges and corals has been extensively characterized due to its ecological importance [8,9]. Sea urchins are considered a key group that maintains the function and stability of marine ecosystems. Because they control the abundance and distribution of other benthic species [10], sea urchins provide different microhabitats for a variety of bacteria that play a fundamental role in shaping and maintaining host homeostasis [11], and include the establishment of native microbiota, digestion, nutrient cycling, health, and immunity [12].

Over the past 20 years, next-generation sequencing has improved our understanding of the diversity of symbiont bacteria and the benefits that bacteria offer to associated organisms, especially by providing nutrients and chemical defenses to their hosts [13]. These studies characterized the bacteria and compared different compartments of the digestive system of Lytechinus variegatus [11] and Srongylocentrotus purpuratus [14]. Changes in bacterial assemblage composition support the essential role of specific bacterial taxa on sea urchin health and digestion. Bacterial specificity depends on the type of feeding of sea urchins, such as the digestion of wood fragments and macroalgae [15,16]. Similarly, other studies have focused on the diversity and activity of sulfur bacteria, as well as sulfur metabolism in the digestive tract of the sea urchin Echinocardium cordatum [17].

The Eastern Tropical Pacific includes the Central Mexican Pacific (MCP) region, which harbors important coral ecosystems that provide a favorable habitat for echinoderms [18]. In this region, the pink sea urchin Toxopneustes roseus (family Toxopneustidae) is one of the most notable and abundant echinoderms of shallow areas [19,20]. T. roseus has extensive ecological contribution; it participates in energy transfer, conducts bioerosion processes, and helps maintain the integrity and persistence of algae in the ecosystem [21,22]. Like many sea urchins of the family Toxopneustidae (i.e., T. pileolus and Tripneustes gratilla), T. roseus contains a source of valuable metabolites such as lectins and bioactive substances with mitogenic and cytotoxic capacity [23,24,25]. Therefore, the sea urchin T. roseus species has biotechnological potential for drug development.

The ecological and biotechnological importance of T. roseus, its distribution along the MCP, and its proximity to human settlements make this species a study model to evaluate the response of the bacterial assemblage to diverse environmental conditions. Therefore, this study aimed to evaluate the spatiotemporal variation of the bacterial assemblage associated with T. roseus in the MCP and its relationship with several environmental variables. The results will contribute to identifying the variations in the structure of bacterial assemblages associated with T. roseus. This knowledge improves understanding of their richness, diversity, and community similarity changes in a specific time and space. Likewise, this knowledge will help establish some bacterial taxa as indicators of the general condition of the community to which they are associated.

2. Materials and Methods

2.1. Study Area

The MCP is located in an oceanographic area limited between two biogeographic provinces, the Warm Temperate Northeast Pacific, and the Tropical Eastern Pacific [26]. This region has mixed conditions resulting from the confluence of currents of different origins [27,28,29]. Due to the geomorphological conditions of the region, the MCP is favored by the discontinuous presence of reef ecosystems with different degrees of coral development, protection status, and human disturbance. The main human disturbances include urban coastal development, fisheries, and tourism [30,31].

The study area covered both insular and coastal sites within the MCP: Nayarit, Jalisco, and Colima (Figure 1), with an approximate distance of 350 km. The sampling sites were: A) Isla Isabel National Park (21° 50′39.62″ N, 105° 52′46.96″ W), the northernmost site; B) Islas Marietas National Park (20° 41′54.1″ N, 105° 34′58.1″ W), C) Islas e islotes de Bahía Chamela Sanctuary (19°33′19.7″N 105°06′28.8″W), D) Bahía Cuastecomates-Punta Melaque (19°13′51.7″N 104°44′01.9″W), E) Carrizales (19°05′44.9″N 104°26′08.7″W) in Bahía Ceníceros, and F) Punto B (19° 5′55.21″ N, 104° 23′24.47″ W) in Bahía Santiago, both on the coast of southern boundary of the MCP.

2.2. Sampling Collection

Samples were obtained in the warm-wet (summer) and cold-dry (winter) seasons of two consecutive years: summer 2017 (August-September), winter 2018 (January-February), summer 2018 (August-September), and winter 2019 (January-February). Three individuals with visual evidence of being healthy of the sea urchin T. roseus were collected for each sampling site and season, for a total of 72 samples. The organisms were transported to land, washed with sterile seawater in an aseptic space, and preserved in cryovials with liquid nitrogen. Each cryovial included the pedicellaria, spines, digestive tract, and gonads of its respective replicate. Once in the laboratory, the cryovials were stored in ultracold conditions (-80°C) until processing.

Variables representative of habitat structure, as well as seawater quality and conditions, were determined for each sampling site and season (Table S1). Sea surface temperature (SST), dissolved oxygen (DO), and salinity (PPM) were recorded in situ in triplicate with a YSI-556 multiparameter probe. Concentrations of chlorophyll α, phosphate (PO₄), silicate (SIO₂), nitrite+nitrate (NO₃+NO₂), ammonium (NH4), total coliforms, and fecal coliforms were measured from seawater samples. Chlorophyll α concentration was estimated following Strickland & Parsons [32], while nutrients were determined following the Skalar San Plus II nutrient analyzer protocol. We determined seawater transparency with a Secchi disk to estimate the light extinction coefficient (turbidity), following the method proposed by Poole & Atkins [33].

The benthic habitat was recorded with video transects (25 m long) filmed 40 cm from the bottom. Each video transect was divided into 40 frames and 50 random systematic points (2,000 points per video) to estimate the composition and coverage of benthic organisms and substrate types. Videos were then analyzed with high-resolution monitors. Structural habitat elements were classified into hermatypic corals, hydrocorals, hydrozoans, octocorals, sponges, fleshy macroalgae, filamentous algae (turf), crustose coralline algae, articulated calcareous algae, and other sessile and vagile organisms. Coverages of other benthic components were designated as sandy substrate, rocky substrate, rubble, recently dead coral, and other benthic elements [34]. Depth (m) was recorded with dive computers, while topographic complexity was determined with the chain method [35].

2.3. Sequence Extraction and Analysis

Bacterial DNA extraction was performed with the modified TRIzol Reagent protocol (Thermo Fisher Scientific). The modifications included an additional Proteinase K treatment (20μl) and a water bath incubation (65°C) for each replicate. The extracted DNA (120 μl) was purified and individually concentrated with the GE Healthcare purification kit (Life Sciences™ illustra™ GFX™). Then, the DNA was quantified by fluorometry in Qubit^® 3.0 (Thermo Fisher Scientific).

Amplification of each replicate aggregate sample was performed by PCR using the modified 16S Metagenomics Kit (Thermo Fisher Scientific) and DreamTaq Green PCR Master Mix (2X) (Thermo Fisher Scientific). The two sets of primers included in the 16S Metagenomics Kit were used to amplify the following hypervariable regions of the bacterial 16S region: V2-4-8, and V3-6, 7-9. The final volume per reaction of the 16S Metagenomics Kit was 15 µL: 7.5 µL of Environmental Master Mix (2X), 1.5 µL of 16S Primer Set (10X), 3 µL of sample, and 3 µL of H₂O. The final volume per reaction of the DreamTaq Green PCR Master Mix (2X) was 12.5 µl: 6.25 µl of DreamTaq Green PCR Master Mix (2X), 1.5 µl of 16S Primer Set (10X), 2 µl of sample, and 2.75 µl of H₂O. Thermal cycler conditions were as follows: 10 min at 95°C, 30 cycles of 30 s at 95°C, 30 s at 58°C, 20 s at 72°C, with a final step of 7 min at 72°C. The presence and quality of amplicons was verified by electrophoresis on 1% agarose gels.

The libraries were constructed using the Ion Plus Fragment Library Kit protocol. To identify the libraries, the ends were phosphorylated, adapters were linked, and barcodes were added to each sample (50-100 ng of DNA). Libraries were quantified with real-time PCR (qPCR) to determine their concentration and calculate the equimolar dilution factor for mixing. The template was prepared with an emulsion PCR on the Ion One Touch 2 System (Life Technologies). The concentration of the template was determined by fluorometry on the Qubit^® 3.0 kit (Thermo Fisher Scientific). The enriched template was loaded onto a PGM 318TM chip, following the instructions of the Ion PGM Hi Q View Sequencing Kit protocol for 400 bp fragments. Sequencing was performed on the Ion Torrent™ Personal Genome Machine^® (PGM) of the National Laboratory of Plant Identification and Characterization (University of Guadalajara).

The Ion Reporter™ Software (ThermoFisher) platform was used to detect and remove primers and short sequences (<150 base pairs) using the 16S Metagenomics Plugin (ThermoFisher). Subsequently, the screened reads were clustered into operational taxonomic units (OTU) with 3% divergence, followed by removal of OTUs with less than ten copies. The final OTUs were taxonomically classified using the curated GreenGenes v13.5 [36] and MicroSEQ(R) 16S Reference Library v2013.1. Data analyses were performed at the family level because the finest and most complete identification of the libraries (99.3% of the sequences) was obtained under this classification. All sequences were uploaded to the National Center for Biotechnology Information repository under project number PRJNA742839.

2.4. Statistical Analysis

To evaluate the change in the structure and composition of the bacterial assemblage between years, seasons, and sites, we used a three-way experimental design with crossed factors (type I fixed-effect model) and no replication as a preliminary analysis. The season factor was removed from the design because it did not contribute to the explained variation of the model (Table S2); in this manner, we used levels to increase replication in a simpler design. Therefore, the spatiotemporal variation of the bacterial assemblage associated with T. roseus was evaluated with a two-way experimental design with crossed factors and fixed effects (type I model), expressed as follows:

$$Y=\mu +{Year}_i+{Site}_j+\left({Year}_i*{Site}_j\right)+{\epsilon }_{ijl}$$

(1)

where Y was the multidimensional response variable or matrix; The Year_i factor presented two levels (periods 2017-2018, and 2018-2019); The Site_j factor corresponded to the six levels or sampling sites (Isla Isabel National Park, Islas Marietas National Park, Islas e islotes de Bahía Chamela Sanctuary, Bahía Cuastecomates-Punta Melaque, Carrizales, and Punto B); The term Year_i * Site_j expressed the interaction between the aforementioned factors, representing the spatiotemporal variation of the model; and ε_ijl was the cumulative error of the model.

The sampling effort was evaluated with sample-based rarefactions to compare the expected bacterial family richness against the nonparametric estimators ICE, Chao 2, Jackknife 1, and Jackknife 2. Rarefaction curves were constructed as a global model for the entire sampling effort by considering replicates from all years, seasons, and sites. The curves were obtained with 10,000 random combinations without replacement with EstimateS V9.1 [37]. Then, they were plotted in SigmaPlot 12.0 (Systat Software Inc., San Jose, CA, USA).

The alpha diversity of the bacterial assemblage associated with T. roseus was analyzed with the estimation of total abundance (number of sequences) per family (N), and Hill numbers (^qD): a) Family richness (⁰D); b) Shannon’s exponential function of diversity (¹D); c) Inverse of Simpson’s dominance (²D); and d) Hill’s equity (²¹D) estimated as ²D/¹D [38]. Indices were estimated per aggregate sample of replicates, and these were compared with a multidimensional, permutational analysis of variance (PERMANOVA) based on the experimental design, using data standardized to Z-values and a Euclidean distance matrix [39].

Beta diversity (composition and abundance of bacterial families) was evaluated with another PERMANOVA, using the same experimental design and statistical evaluation method but based on a matrix in Hellinger distances [40]. The contribution of bacterial taxa to the average dissimilarity between years and sites was estimated with a similarity percentage analysis (SIMPER) with Hellinger distances at a cut-off of 65% cumulative contribution to the average dissimilarity.

The variation of environmental variables was analyzed with a third PERMANOVA, using the experimental design described before but without replication. To do this, we used the data at the site level per year standardized to Z values, with which we constructed a matrix in Euclidean distances [39]. Then, a SIMPER analysis was performed with Euclidean distances and a cut-off at 65% of the cumulative contribution of the environmental variables, using a two-way design crossed with the factor Year and factor Site.

The results of the PERMANOVA for alpha and beta diversity were visualized in ordinates constructed with principal coordinate analysis (PCO). However, the third PERMANOVA (environmental variables) results were represented in a non-metric multidimensional scaling (nMDS). All ordinations were based on the pretreatments and similarity coefficients of their respective PERMANOVA. The alpha (Hill abundance and numbers) and beta (important bacterial families selected based on SIMPER results) diversity variables were represented as vectors based on multiple correlations. Spearman correlations were used for the environmental variables (relevant variables according to the SIMPER analysis). The statistical significance of the overall and after the PERMANOVA tests was calculated with 10,000 permutations of residuals under a reduced model and type III sum of squares. When less than 100 permutations were obtained, Monte Carlo (MC) tests were used. Rarefaction curves, Hill number estimates, and ordinate analyses (PCO and nMDS) were developed in PRIMER v6.1.11 + PERMANOVA v1.01 [41].

The relationship between bacterial assemblages and environmental variables was evaluated at the sampling site level with a canonical redundancy analysis (RDA) [42]. The composition and abundance of the bacterial families previously treated with a Hellinger transformation were used as the biological matrix (Y), while the matrix of predictor variables (X) corresponded to the environmental variables previously standardized to Z-values (Table S1). The multicollinearity of the predictor variables was reduced by considering those variables that presented a Pearson correlation (r) ≤ 0.75, and a variance inflation factor (VIF) < 10. The adjustment of the RDA model was calculated as the adjusted coefficient of determination (R²_adj), based on the Trace statistic. The correlation among sites, bacterial families, and predictor variables was plotted in an ordination (triplot). In addition, a second RDA model was performed at the Site per Year level with the previous method (Table S1). This analysis was complemented by the results of the two-way SIMPER with crossed factors. The statistical significance of both RDA models was determined with 10,000 Monte-Carlo permutations under a reduced model in CANOCO v4.5 software [42].

3. Results

3.1. Characterization of the Bacterial Assembly of T. roseus

The sequencing generated a total of 2,023,663 reads, of which 99.3 % (2,009,347 reads) were identified at the family level. Sampling effort ranged between 70.5 % - 79.9 %, representing 76.4% of the expected total richness of bacterial families according to nonparametric estimators (Table S3, Figure S1). The 1,855 OTUs obtained were classified into 190 families, 80 orders, 39 classes, and 19 phyla.

At the temporal level, the phylum Proteobacteria contributed the highest to relative abundance (67.0 %-71.1 %), followed by Bacteroidetes, and Firmicutes (Figure 2A). In total, 109 families were shared between years, with 113 families recorded in the first year and 189 families in the second year of sampling. The dominant families were Helicobacteraceae (42.9 %- 46.9 %), and Vibrionaceae (6.4 %-7.1 %) (Figure 2B). Other predominant families were Burkholderiaceae (4.1 %-5.3 %), Campylobacteraceae (2.2 %-6.2 %), Flavobacteriaceae (2.6 %-5.8 %), and Acholeplasmataceae (2.1 %-3.9 %).

From a spatial perspective, the six sites were characterized by the dominance of the phylum Proteobacteria, registering relative abundances between 38.6 %-91.1 % (Figure 2C). Bacteroidetes (10.8 %-33.7 %) was the second abundant phylum particularly at Punto B (PB), Islas Marietas National Park (PNIM), and Isla Isabel National Park (PNII). Firmicutes (7.0 %-10.3 %) predominated in Bahía Cuastecomates-Punta Melaque (CUM), and Islas e islotes de Bahía Chamela Sanctuary (BCH), while in Carrizales (CRZ) the highest abundance was represented by the phylum Tenericutes (5.19 %). Also, 23 families were shared among the six sites. A total of 181 families were recorded in PNII, 109 in PNIM, 55 in BCH, 61 in CUM, 59 in CRZ, and 43 in PB. The most abundant family at five of the six sites was Helicobacteraceae (37.5 %-64.9 %), whereas Saprospiraceae (16.1 %) predominated at PNII (Figure 2D). The families Vibrionaceae (2.9 %-15.2 %), Flavobacteriaceae (1.3 %-13.1 %), Desulfovibrionaceae (0.8%-15.2%), Campylobacteraceae (1.8 %-6.8 %), and Spiroplasmataceae (1.8 %-13.9 %) were also relevant for their high abundances.

The bacterial families associated with T. roseus showed important changes in the PB and PNII sites (Figure S2) on a temporal scale. In the first year, PB recorded relatively high values of Helicobacteraceae (73.53 %), changing in the second year by the predominance of Desulfovibrionaceae (21.55 %) and Flavobacteriaceae (28.72 %). A similar pattern was observed in PNII, where the Flammeovirgaceae (19.63 %) and Acholeplasmataceae (15.58 %) families represented the most abundant component in the first year, while Saprospiraceae (22.38 %) and Xanthobacteraceae (8.64 %) dominated in the second year.

3.2. Alpha and Beta Diversity

The PERMANOVA model of alpha diversity (N, ⁰D, ¹D, ²D, ²¹D) explained 51.7 % of the total variation; the Site factor was the only one that presented significant differences (37.3 % of explained variation) (Table 1). The a posteriori tests showed that the spatial difference is given by the PNII with respect to the other sites, except for the PNIM (Table S4); these two island sites are located in the northern limit of the study region. The PCO ordination of bacterial alpha diversity explained 98.1 % of the total variation at the site level (Figure 3A). The first principal coordinate (PCO1) was associated with the highest ⁰D, ¹D, and ²D values characteristic of the PNII. In contrast, the second principal coordinate (PCO2) was related to high N, found in the PNIM. The remaining sites, especially PB and BCH, correlated with ²¹D.

The PERMANOVA model of bacterial family composition and abundance explained 53.5 % of the total variation (Table 1). The individual factors Year and Site presented significant differences, and spatial change was the most important (19.5 % of the variation explained). During the second sampling year, the interannual variation was generated mostly by the increase in the abundances of the families Helicobacteraceae, Vibrionaceae, Desulfovibrionaceae, Saprospiraceae, Xanthobacteraceae, Planococcaceae, and Spiroplasmataceae. In addition, the abundance of the families Flammeovirgaceae and Bartonellaceae decreased. The results of the a posteriori tests of the Site factor were similar to those of the alpha diversity, showing that the bacterial assemblage of the PNII was different from the other sites, except for PB (Table S4). According to SIMPER analysis, the families with the highest contribution to average dissimilarity both in years and sites were Helicobacteraceae, Desulfovibrionaceae, Burkholderiaceae, Flavobacteriaceae, Spiroplasmataceae, and Desulfonatronumaceae (currently Desulfonatronaceae) (Table S5). The families Desulfovibrionaceae and Spiroplasmataceae were only recorded in the second year (2018-2019), when Helicobacteraceae, Flavobacteriaceae, and Desulfonatronaceae showed higher relative abundance (Figure 2C). Similarly, the families Helicobacteraceae, Desulfovibrionaceae, and Desulfonatronaceae were representative at all sites except in PNII, whereas, Burkholderiaceae and Spiroplasmataceae were found in low abundances at PNII and PB sites. Flavobacteriaceae was present at most sites (except CUM), recording its highest relative abundance at PB (Figure 2D). The spatial patterns of the bacterial assemblage were observed on a PCO ordination that explained 75.7 % of the total variation (Figure 3B). High abundance values of Desulfonatronaceae associated mainly with PB and BCH were observed at the first principal coordinate (PCO1). In contrast, Helicobacteraceae had higher representation at PNIM, CRZ, and CUM sites in the second principal coordinate (PCO2).

The PERMANOVA model of the environmental variables explained 52.3 % of the total variation (Table 1). The Site factor was the only one that presented significant differences, with 34.4 % of the explained variation. However, subsequent tests did not allow to determine significant differences (P>0.05; Table S4). According to SIMPER analysis, the environmental variables with the greatest contribution to spatiotemporal differentiation were α-chlorophyll, dissolved oxygen, salinity, ammonium, fecal coliforms, and the cover of live coral, sponge, rock, and sand (Table S6). The nMDS ordination showed a stress value of 0.15 (Figure 3C) and projected spatial distinctiveness of the PNII; this site positively correlated with a high sponge cover. Temporal differences in PB were attributed to fecal coliforms (Figure 3C).

3.3. Relation of Bacterial Assemblage and Environmental Variables

Spatial RDA ordination was shown to be a significant model (P-value=0.0437) and R²_adj=62.0 % (Figure 4A). The best combination of variables predictive of spatial change in the bacterial assemblage consisted of depth, and live coral, sand, and sponge covers (Figure 4A, Table S1). Sponge cover was the only significant individual environmental variable (Table S8). The PNII was characterized by the highest sponge cover and average depth, which were positively related to higher abundances of Desulfobacteraceae and Flammeovirgaceae families and negatively related to the Helicobacteraceae family (Figure 4A). The highest live coral cover correlated with CUM and CRZ, followed by PNIM. These sites were mainly associated with the families Spiroplasmataceae and Thermoactinomycetaceae. In contrast, sites BCH and PB presented the highest sand cover. BCH was associated with high abundances of Campylobacteraceae and Desulfonatronaceae. PB correlated positively with Enterobacteriaceae, Desulfovibrionaceae, and Flavobacteriaceae; however, it correlated negatively with Burkholderiaceae and Clostridiales Family XI Incertae Sedis.

The spatiotemporal RDA resulted in a significant model (P-value = 0.0002) with an adjustment of R²_adj=67.73 % (Figure 4B). The best combination of predictor variables for this model consisted of environmental variables, such as sea surface temperature, salinity, dissolved oxygen, fecal coliforms, and light extinction coefficient, and benthic conditions, such as the presence of sponges, macroalgae, and sandy substrate (Figure 4B, Table S1). Sponge cover, salinity, and dissolved oxygen were significant on their own (Table S8).

The spatiotemporal model complemented the results of the spatial RDA. In particular, these models highlight the spatial differentiation of PNII and the temporal differentiation of PB and BCH in the second year. PNII was characterized by the highest average sponge and macroalgae cover. In the first year, PNII correlated with the families Methylobacteriacaee, Desulfobacteraceae, and Flammeovirgaceae, among others. However, in the second year, it was positively associated with Defluviitaleaceae and Saprospiraceae (Figure 4B). In addition, the second year of sampling recorded higher average fecal coliform values in PB and BCH; this positively related to an increase in the Desulfovibrionaceae, Desulfonatronaceae, and Enterobacteriaceae families. In contrast, for both years, the other sites correlated positively with a high abundance of Helicobacteraceae and Campylobacteraceae, but negatively with high values of dissolved oxygen, temperature, fecal coliforms, and sponge coverages.

4. Discussion

The results of this study evidenced that the structure and composition of the T. roseus bacteriome in the Mexican Central Pacific (MCP) showed certain spatial and temporal “stability” in most sites (PNIM, BCH, CUM, CRZ, and PB). These results agree with those reported by Hernández-Zulueta et al. [43], who found no spatial and temporal variation of the bacterial microbiota of Pocillopora damicornis and P. verrucosa in the MCP. These results suggest that bacterial assemblages are not influenced by physicochemical variables such as sea surface temperature, considering that the region presents daily variations of ±3°C [44] and an annual temperature range of ±10°C [45]. Littman et al. [46] and Carlos et al. [47] also observed this pattern in corals from the Great Barrier Reef in Australia and Brazil, respectively. These results are relevant because bacterial microbiota has been reported to contribute to the adaptation, resilience, and resistance of marine invertebrates to thermal stress generated largely by vents such as El Niño Southern Oscillation (ENSO) [48,49]. Therefore, the spatial stability of bacterial assemblages could be due to urchin-microbiota specificity. This specificity might reflect the microbiota’s functions to maintain the host’s health, such as digestive health and innate immunity [12].

The only site that showed significant spatial and temporal change in the bacterial assemblage of T. roseus was Isla Isabel National Park (PNII). The PNII is an insular site located at the boundary of two biogeographic provinces [26], which promotes the richness and abundance of echinoderms representative of temperate and tropical zones [10,50]. Finally, the fact that it is an island site 28 km from the coast [43] could contribute to the differences in the T. roseus bacterial assemblages between the PNII and the rest of the sites sampled in the MCP.

Regarding the temporal difference of the T. roseus bacteriome in the PNII, a positive correlation was observed in the Flammeovirgaceae family with the macroalgal coverage of the PNII in the first year analyzed (2017-2018). Increases in this family of opportunistic bacteria have been associated with prolonged thermic stress [51] and disease in asteroids [52]. However, this pattern changed in the following year (2018-2019), with the predominance of the Saprospiraceae and Xanthobacteraceae families. The Saprospiraceae family has been associated with bald disease in sea urchins [13], while an abundance of Xanthobacteraceae has been found in bleached algae [53]. Therefore, our results suggest that the high values of Flammeovirgaceae and Xanthobacteracea families in PNII might be remnants of the 2014/2016 ENSO event [54].

With the exception of the PNII, the T. roseus bacteriome maintained a similar pattern across the years studied (2017-2018 and 2018-2019), with no differences associated with local conditions [55,56]. These findings support the hypothesis of spatial and temporal adaptability of species to their environment [57]. This adaptability might be the result of a coevolutionary process [58] or the reflection of the host’s ability to structure its bacterial microbiota, modifying the abundance of specific taxa obtained from the environment, regardless of site-linked variations [55].

The bacteriome was primarily characterized by the prevalence of Proteobacteria and Bacteroidetes. These phyla are the most abundant members of the bacterial microbiota of echinoids [59,60,61]. In the present study, they represented the major component of the bacterial assemblage of T. roseus, collectively contributing to 43.6 %-68 % of the bacterial relative abundance of the sea urchin. Particularly, these families dominate in the coelomic fluid, feces, and intestinal tissue of some echinoid species [12,14,62], and they indicate the mutualistic relationship between the sea urchin and its resident microbiota [11] to carry out digestion and nutrient uptake processes [12].

Furthermore, echinoderm assemblage in the MCP is strongly related to local conditions, in particular, to structural elements of the benthic habitat [63]. Therefore, variables representative of habitat structure explained the spatial change of the bacterial assemblage of T. roseus. Live coral cover (LCC) was an important predictor variable for T. roseus bacterial assemblage. CRZ had the highest average LCC (67.6 %) of all sites, while BCH, CUM, and IMNP presented values between 16 %-21 %. Hernández-Zulueta et al. [64] observed that the environmental variables that explain the variation in bacterial assemblages of mucus and tissue of P. damicornis and P. verrucosa, were the coverages of fleshy macroalgae, live coral, and sponges. These results suggest that the coral bacterial assemblage is interconnected with the bacterial microbiota of other structural elements of the benthos.

Due to the close relationship of sea urchins with the sediment, the sandy substrate represented a key component of the bacterial assemblage of echinoids. The bacterial microbiota of marine sediment is highly diverse and different from other organisms [59,60]. Thus, the sediment acts as a reservoir of microorganisms that the holobiont can obtain and filter, selecting specific taxa according to the internal conditions of the sea urchin [55]

Punto B, the site with the least coral development (LCC<5%) and the greatest anthropogenic impact, showed the highest abundance of the Enterobacteriaceae family, which has been associated with fecal coliforms [65,66]. In addition, this site is in front of a large urban development and receives constant wastewater discharges [43], leading to changes in bacterial assemblages [65,66]. Punto B also showed dominance of Prolixibacteraceae and Sphingobacteriaceae during the second year of sampling. Prolixibacteraceae has been reported as the second most abundant family in wetlands of contaminated sites due to their proximity to coastal urban developments [67]. While Sphingobacteriaceae includes members resistant to physical disturbance and heavy metals [68,69,70]. In addition, Flavobacteriaceae and Desulfovibrionaceae presented high relative values in PB and BCH from 2018 to 2019. Although Flavobacteriaceae is found in low amounts in urchins [11,13], this family predominates in sea urchins inhabiting near urban developments [62]. Similarly, other marine invertebrates, e.g., corals, with a disease or some level of environmental stress show an overrepresentation of some bacterial taxa, including Flavobacteriales and other members of Desulfovibrionaceae [71,72].

This study presents the first spatiotemporal characterization of the bacterial assemblage associated with the sea urchin T. roseus in the MCP. Our results indicate that the changes in environmental conditions derived from seasonality do not represent a significant factor influencing the T. roseus bacteriome at the regional level. The structure, composition, and abundance of the T. roseus bacterial assemblage at most sites (PNIM, BCH, CUM, CRZ, and PB) suggested the plasticity of the holobiont in the face of spatial and temporal environmental variability in the MCP. In addition, the dominance of some families of the order Campylobacterales (i.e., Helicobacteraceae and Campylobacteraceae) is critical because of their possible relationship with urchin health at the physiological level. However, the environmental differences of the PNII site suggest that the bacterial microbiota of the pink sea urchin T. roseus may be dynamic and present important changes in its composition and abundance.

Echinoid diseases have not yet been reported in the MCP. However, the particular prevalence of some bacterial taxa implicated in diseased or stressed marine organisms (i.e., Flavobacteriaceae, Sphingobacteriaceae, Saprospiraceae, and Flammeovirgaceae, among others), suggests that some populations of T. roseus are under stress, mainly in the PB and PNII sites. In addition, the low representation of the Helicobacteraceae and Campylobacteraceae families in the PB and PNII sites could indicate unhealthy hosts. The dominance of Enterobacteriaceae suggests that PB and BCH coral ecosystems (2018-2019) were exposed to possible anthropogenic contamination.

Finally, this work showed that reef habitat structure significantly influences the local bacterial assemblage of T. roseus. The relationship we found between high live coral cover (LCC) and sites with similar assemblages emphasizes the importance of continuing efforts to protect and conserve the main coral ecosystems in the MCP.