Microbiota of the Digestive Glands and Extrapallial Fluids of Clams Evolve Differently Over Time Depending on the Intertidal Position

The Manila clam (Ruditapes philippinarum) is the second most exploited bivalve in the world but remains threatened by diseases and global changes. Their associated microbiota play a key role in their fitness and acclimation capacities. This study aimed at better understanding the behavior of clam digestive glands and extrapallial fluids microbiota at small, but contrasting spatial and temporal scales. Results showed that environmental variations impacted clam microbiota differently according to the considered tissue. Each clam tissue presented its own microbiota and showed different dynamics according to the intertidal position and sampling period. Extrapallial fluids microbiota was modified more rapidly than digestive glands microbiota, for clams placed on the upper and lower intertidal position, respectively. Clam tissues could be considered as different microhabitats for bacteria as they presented different responses to small-scale temporal and spatial variabilities in natural conditions. These differences underlined a more stringent environmental filter capacity of the digestive glands.


Introduction
Host-associated microbiota are organized into trophic groups of microorganisms that interact with each other and the organ they colonize. These microorganisms come from the parents (vertical transmission) or from the immediate environment (horizontal transmission). The colonized tissues offer particular physicochemical conditions and act as microhabitats. Both the host and its associated microbiota contribute to the metabolic specificity of these organs. As a result, each tissue is colonized by specific microbiota communities that are both complex and dynamic [1]. Associated microorganisms contribute to varying degrees to many metabolic processes and gene expressions that can have major impacts on host nutrition, detoxification, growth, immune system development, and resistance to infections [2][3][4]. They probably play a key role in hosts acclimation to different ecological conditions in fluctuating environments as coastal zones.
It has been shown that marine host-associated microbiota may be influenced by different biotic and abiotic factors (season, location, host age, health, and genotype) [5][6][7][8]. Being filter feeders, bivalves interact with numerous microorganisms coming from the environment. They maintain tissuespecific microbiota that differ from that of their surrounding environment [9], indicating that both organisms and tissues select their bacteria [10,11]. Nonetheless, these host-associated microbiota may vary spatially and temporally reflecting seasonal changes and geographic locations [12][13][14][15][16]. Intrinsic factors, like life stage, genetic, host tissue, or host taxonomy may also induce different variations of microbiota composition and dynamic [17][18][19]. However, the mechanisms underlying the interactions in bivalve microbiota and how these later are likely to modify their host acclimation to a changing environment are still largely unknown.
The Manila clam (Ruditapes philippinarum) is a sediment burrower that mainly ingests benthic microalgae and sedimented phytoplankton [20]. Its microbiota is influenced by environmental factors such as season [21], contaminants [22], and intertidal position [11]. They also differ among the heart, gonad, mantle, hemolymph, gills, or gut [21,23], indicating an organ-specific microbiota. This suggests that environmental factors may have different effects on microbiota dynamics depending on the considered tissues and so differently affect their functions.
Almost nothing is known about clam's extrapallial fluids (EF) microbiota, although this compartment plays an important role in immunity and calcification processes [24], which could be influenced by microbial changes. Indeed, it has been previously shown that microbiota from extrapallial fluids of eastern oyster (Crassostrea virginica) is influenced by location and season [10].
Microbiota is important in adaptation to fluctuating environment like the intertidal zone, and bivalves are good models for that. To highlight specific responses of tissue-associated microbiota in the face of environmental small-scale variations, we investigated the dynamics of clam microbiota during the spring revival, a period that is characterized by marked environmental changes impacting the host species [25]. Animals were deployed at two tidal levels and were sampled in February, March, and April. Two tissues that are differently exposed to the environment were compared (with and without depuration), the EF, which are at the interface between the animal and the environment, and the digestive glands (DG), an internal organ that is exposed to what is ingested. This study will provide a better understanding of the dynamic and scales of variation of microbiota associated with the DG and EF of healthy clams in natural conditions.

Experimental design
The Manilla clams Ruditapes philippinarum were hatched in April 2016 and were provided by the SATMAR hatchery (France). In October 2017, clams were placed at the site IFREMER Br08 located in the Bay of Brest at Pointe du Château (48° 20′ 06.19″ N, 4° 19′ 06.37″ W, Brittany, France). Clams were arranged in two duplicated mesh bags at two different intertidal levels corresponding to 20% (L20) and 56% (L56) of time spent emerged at low tide (Fig. 1).
Healthy clams (i.e., without signs of brown ring disease) were collected from the two intertidal levels over three consecutive days (to avoid sampling bias) in T1 (February 2018, 19.6 ± 3.4 mm), T2 (March 2018, 18.1 ± 2.4 mm), and T3 (April 2018, 19.5 ± 2.8 mm). For each sampling day, sediments were sampled in triplicate at each level, next to the clams, and seawater was collected in triplicate 2 h before low tide. The seawater samples (1L) were successively passed through 8 (SW8) and 0.22 μm (SW022) polycarbonate filters (Whatman, USA) that were kept at − 80 °C until DNA extraction and sequencing. For each day, clams (n = 10) were collected by level for bacterial microbiota analyses, and half of them (n = 5) were directly dissected for tissue sampling, while others (n = 5) were transferred to the laboratory and grouped by day and intertidal level for depuration process Fig.1 Representation of experimental design detailing the deployment of clams on the two intertidal levels (L20 and L56 corresponding to 20% and 56% of time passed in emerged sediment, respectively) in order to empty the DG and reduce environmental microorganisms [26]. After 14 days in tanks (30 L) fed with flowthrough (3 L/min) filtered seawater (10 and 5 μm sand filters and UV treatment before two 1 µm filters and a second UV treatment), tissues were dissected and collected before DNA extraction. Tanks were cleaned every day to avoid biofilm formation and bacterial proliferation. No feed was added.

Bacterial Genomic DNA Sequencing
Microbiota analyses were performed from digestive glands (DG) and extrapallial fluids (EF) of both depurated and nondepurated clams, as detailed by Offret et al. [11]. Briefly, DG and EF were aseptically dissected from clams directly after sampling or after depuration. Bacterial genomic DNA from both clam tissues and environmental samples was extracted using the PowerLyzer PowerSoil DNA Isolation Kit (Qiagen, USA). For each sample, the DNA was amplified by PCR targeting the variable V3V4 region coding for the bacterial 16S rRNA gene using the 341F (5′-CCT ACG GGNGGC WGC AG-3′) and 805R (5′-GAC TAC HVGGG TAT CTA ATC C-3′) primers [27]. Finally, paired-end sequencing with a 300-pb read length was performed at McGill University (Génome Quebec Innovation Centre, Montréal, QC, Canada) on a MiSeq system (Illumina).

Bacterial Microbiota Analysis
Sequences were processed with the Galaxy-supported FROGS pipeline ("Find, Rapidly, OTUs with Galaxy Solution"). This pipeline was developed in order to analyze large sets of amplicon sequences and produce abundance tables of operational taxonomic units (OTUs) with their taxonomic affiliation [28].
Firstly, bacterial 16S rRNA gene sequences of samples were preprocess. Paired-end reads were merged using Flash [29] with 10% mismatch tolerated. Primers and adapters were removed using cutadapt [30]. Clustering was achieved using the Swarm method with an aggregation distance (d = 3) [31]. Afterwards, chimeras were detected and removed with VSEARCH (de novo UCHIME method) before applying a cross-sample validation step [28,32,33]. Next, a relative abundance filter was applied on OTUs with a threshold of 0.005% (of global relative abundance), as recommended by Bokulich et al. [34]. Taxonomic affiliation was done on this new data by blastn + [35], using the Silva 132 16S database [36] formatted to include taxonomic levels up to species [28]. For analysis, cyanobacterial and eukaryotic sequences (chloroplast and mitochondrial 16S rRNA gene from algal cells) were removed.

α-Diversity
Microbiota α-diversity was assessed at the OTUs level, and the rarefied richness and Shannon indices were computed. The period, level, tissue, and depuration effects (and their interactions) were tested on clam's α-diversity indices. For the environmental samples, the effects of period and level (and their interaction) were tested on the sediment and the effects period and filter (and their interaction) on the seawater. For each model, the effects were tested using a linear model on Box-Cox transformed data and a generalized linear model (GLM) on non-transformed data. The linear model was either a two-way analysis of variance (ANOVA) or an analysis of covariance (ANCOVA) with the number of sequences as a co-variable. Finally, the models with the greater adjusted R 2 were kept.

β-Diversity
The relative abundance of the whole bacterial microbiota communities was assessed at the order level at each period [8]. A Kruskal-Wallis test was used to highlight significative modifications (P < 0.05) of relative abundances of bacterial orders between sample groups. In relative abundances compared with microbial communities, β-diversity was also assessed at the OTUs level. On the one hand, the OTUs relative abundance tables were transformed using the Box-Cox-chord transformation with exponent 0 which corresponds to a log-chord transformation [37]. Then, the variation in microbiota composition and structure between all samples was visualized with a PCA of these transformed OTUs relative abundances. The effects of period, level, and their interactions were tested using PERMANOVA [38] when group dispersions were homogeneous. Multivariate homogeneity of group dispersions was assessed for the DG, EF, sediments, and seawater samples separately. All analyses were carried out using R3.6.2 [39], with functions from the phyloseq [40], vegan [41], emmeans [42], RAM [43], and VIM [44] packages.

Sequences Analysis Information
The total number of sequences before processing was 13,406,675. After processing, 7,369,012 sequences were kept (i.e., 54.97%). After removing sequences affiliated to chloroplasts and samples with a number of sequences inferior to 3000, 432 samples were kept: 160 DG's samples,

α-Diversity Comparison for Clam Tissues and Environmental Samples
The effect of the number of sequences was significant for observed richness (F = 7.6, P < 0.01) justifying the use of a rarefied richness for describing α-diversity patterns. Models applied to clam's tissues microbiota α-diversity indexes (rarefied richness and Shannon indices) highlighted two main factors: depuration and tissue (including the seawater fractions), followed by the period (Fig. 2). The depuration induced a significant decrease of rarefied richness (F = 4.7, P < 0.01) and Shannon indices (F = 62.2, P < 0.001), while tissue differences resulted in an average EF α-diversity that was always significantly higher than those for the GD (rarefied richness: F = 29.4, P < 0.001; Shannon: F = 14.6, P < 0.001). The two seawater fractions exhibited two clearly different α-diversity, where the bacterial communities from SW8 fraction were always greater than those from the SW022 one (rarefied richness: F = 211.9, P < 0.001). Nonetheless, on average α-diversity of microbiota associated with the two clam tissues was always significantly lower than those of bacterial communities from sediment and SW8 (post hoc Tukey test, P < 0.001).
The period factor was significant for all models, except for the Shannon indices of sediment (F = 1.4, P > 0.05). Overall, regarding both indices, average α-diversity of bacterial communities from non-depurated DG and seawater samples significantly decreased during the studied period, notably between T1 and T3 (post hoc Tukey test, P < 0.001). In the other hand, non-depurated EF and sediment α-diversity significantly increased between T1 and T3.

Order Level Bacterial Diversity from Clam Tissues and Environmental Samples
Based on the 30 most abundant bacterial orders, DG, EF, sediment, and seawater fractions showed clearly distinct bacterial communities in terms of both relative abundance and taxa composition (Fig. 3). The two fractions of seawater (SW8 and SW022) were dominated by the orders Rhodobacterales (55-70% in SW022 and 20-30% in SW8) and Flavobacteriales (12-15% in SW022 and 10-22% in SW8). However, their composition differed for less abundant orders (less than 6% of total relative abundance). Indeed, SW022's  , depuration, period). Asterisks indicate significant pairwise comparisons between the three periods (*, P < 0.05; **, P < 0.01; ***, P < 0.001, post hoc Tukey test) bacterial community was mostly composed of Oceanospirillales and Cellvibrionales, whereas SW8's bacterial community was more diverse with Desulfobacterales, Campylobacterales, and Myxococcales, notably present in low relative abundance (inferior to 5%). Sediment bacterial communities were mostly composed of numerous few abundant orders, accounting for more than 25% of the total relative abundance, while the most abundant orders were mainly affiliated to Desulfobacterales (13-19%) and Actinomarinales (6-8%).
The non-depurated clams' DG microbiota was mostly composed of a decreasing order of contribution of Rickettsiales, Mycoplasmatales, and Diplorickettsiales. DG microbiota from clams sampled in T1 was characterized by a lower relative abundance of Rickettsiales and a higher relative abundance of Spirochaetales and Clostridiales, while higher relative abundances of Campylobacterales and Flavobacteriales were observed respectively in T2 and T3. One of the most abundant order found in DG microbiota, the Diplorickettsiales, drastically decreased in T3 compared to T1 and T2. In depurated clams, the DG microbiota was globally composed of the same orders as non-depurated clams but in different relative abundances. The main changes were observed for Spirochaetales, increasing from around 6 to 20% between the non-depurated and depurated clams, while Pirellulales relative abundances were divided by 7 in depurated clams.
Although it harbored some taxa (at the order level) in common, EF microbiota of non-depurated clams was clearly different from that of the DG's (F = 56.661, P = 0.001, 999 permutations). It was mainly composed of Oceanospirillales, Spirochaetales, and Rickettsiales. Both Rhodobacterales and Francisellales orders clearly contributed to the differences with DG microbiota. In T1, the EF microbiota was characterized by a higher relative abundance of Spirochaetales and Oceanospirillales and a lower relative abundance of Alteromonadales, Sneathiallales, and Rhodobacterales. A higher relative abundance of Campylobacterales and Flavobacteriales were observed in T2 and T3, respectively. In EF microbiota from depurated clams, the same orders were globally found than in non-depurated clams but in different relative abundances, notably for Alteromonadales, Francisellales, and Pseudomonadales which were present in much

Factors Influencing Beta Diversity
From this point on, results are those of analyses conducted at the OTU's level. PCA (Fig. 4) allow clearly discriminating sample matrix (F = 46.511, P = 0.001, 999 permutations) (clams vs sediment and seawater on PC1), seawater fractions (SW8 vs SW022 on PC2), and clam tissues (EF vs DG on PC2). The SW8 bacterial community was more similar to that of sediments, whereas the one from SW022 was more similar to that of the EF.
Homogeneity of variances was previously tested for each sample grouping, revealing that multivariate dispersion was not homogenous (Fig. 5) for both DG (F = 17.333, P < 0.001, 999 permutations) and EF tissues (F = 6.893, P = 0.001, 999 permutations). Using the betadisper function, the DG microbiota had a higher heterogeneity of group dispersions than EF. The heterogeneity of dispersions occurred between the depurated samples of L56 in T3 and all the groups in T1 (Fig. 5A). It also occurred between the non-depurated samples of L20 in T2 and all the groups in T2. The dispersions were also heterogeneous between T1 and T3. In the case of the EF's microbiota (Fig. 5B), the heterogeneity of dispersions implied three groups of clams: the depurated at L20 in T1, the non-depurated at L56 in T2, and all the depurated in T3. The L20 levels of T2 and T3 seemed more similar to T1 (L20 and L56) than to L56. For the depurated clams, all clusters were overlapping, indicating a more similar community over the three periods. However, L56 and L20 of the three periods tended to pull away from each other and the centroid, indicating differences between levels.
Regarding the DG bacterial microbiota (Fig. 5a), a high heterogeneity of dispersion, which increased between T1 and T3, was observed. However, it did not prevent us from seeing the dynamic of the microbiota. A clear distinction between depurated and non-depurated clam's microbiota was observed along the first axis of the PCA for both levels (L20 and L56) in T1 and T2, which was much less the case in T3. For the non-depurated clams, the DG bacterial microbiota in T1 and T2 were very close regardless of the levels; no period or level effects were observed. It is only in T3 that the DG bacterial microbiota detached (especially for L20) from T1 and T2, and that a slight level effect appeared. For the depurated clams, the DG microbiota were more similar between the three periods. Moreover, in that case, the level effect was more pronounced in T2 than in T3.
For the EF bacterial microbiota (Fig. 5b), compared to the DG microbiota, the distinction between depurated and non-depurated was clearer along the first axis, as clusters did not overlap in the PCA regardless of levels and periods. The EF microbiota behaved differently from the DG ones. In T1, no level effect was visible (as for the DG), but in T2 and T3, the level effect becomes more important than the period effect, with a clear distinction between levels regardless of the time period. The L56 levels were clearly distinct from the other ones in both T2 and T3.
The PCA representing the sediments communities (Fig. 5c) showed a level effect but no period effect. Only the L56 samples in T1 appeared to have a distinct microbiota from other L56 samples. Multivariate analysis of variance on the Box-Cox transformed and standardized environmental data reveals a period effect (F = 1.77, P = 0.025, R 2 adj = 0.069), the level (F = 9.25, P = 0.001, R 2 adj = 0.18), and their interaction (F = 1.59, P = 0.041, R 2 adj = 0.062). The seawater's bacterial community (Fig. 5d) was clearly different between the three periods and the two filters. The heterogeneity of dispersion between the groups did not allow a multivariate analysis of variance (F = 3.6108, P = 0.037, 999 permutations).

Discussion
In this study, we investigated the diversity and the temporal dynamic of the DG and EF microbiota from Manilla clams placed on two intertidal levels in the foreshore during spring revival. Results brought more information regarding specificity of clam tissue microbiota, while both small spatial and temporal scales significantly contributed to the shaping of the DG and EF microbiota. The observed differences were attributed in a decreasing order to the depuration process,

The Environmental Contribution to the Clam Holobiont
Regarding clam microbiota, both DG and EF were clearly different from the environment, although the bacterial communities in the surrounding seawater and sediments contributed to the microbiota shaping. Indeed, EF microbiota was much closer to environmental bacterial communities than DG microbiota, especially from the 0.22 to 8 μm seawater fraction, suggesting a greater sharing between EF and seawater (free and small particles associated bacteria) than with the digestive glands. OTUs present in the clams DG microbiota, such as Rickettsiales, Mycoplasmatales, Diplorickettsiales, Spirochaetales, Oceanospirillales, and Flavobacteriales, are commonly found in clams DG [11,22], with Mycoplasmatales and Rickettsiales attributed to the core members [22]. EF microbiota harbored a high relative abundance of taxa orders such as Rhodobacterales, similar to those found in seawater but not in the DG. This is probably the result of a direct exchanges between seawater and EF, maybe through the pallial line, possibly associated with less stringent EF physical and chemical conditions for the seawater bacteria than in the DG. To the author's knowledge, EF microbiota of Manila clams has never been described so far. However, Rickettsiales and Flavobacteriales were already found in high relative abundance in EF microbiota from another bivalve species, the eastern oyster Crassostrea virginica [10]. Symbol colors and shapes in legend separating the sampling period combined to intertidal level and the depuration process or seawater fractions

Tissue-Specific of Clam Microbiota Dynamics
This study revealed a different response of the two-tissue microbiota over time. Even if both DG and EF microbiota showed a strong temporal variability in terms of α-diversity, EF's α-diversity significantly increased over time in their natural environment, whereas that of the DG's decreased. Indeed, organ microbiota are shaped by the organ's physiological role, notably digestion and shell mineralization for DG and EF, respectively (Wilbur and Saleuddin, [45]; Allam and Paillard, [46]. On the one hand, the Manila clam's different organs are known to be bacterial niches with unique features [21], and so they could be considered as different microhabitats. On the other hand, this could be linked to environmental variations, such as temperature and dissolved oxygen which have already been shown to impact EF microbial composition of Pacific oysters [10]. Regarding DG microbiota, a higher α-diversity has been described from hepatopancreas microbiota of the Manila clams sampled in winter than in summer [22]. We hypothesize that the observed decrease of bacterial α-diversity in clam DG could be caused by the presence of other water-associated microorganisms such as micro-eukaryotes. Indeed, chlorophyl concentration recorded on the intertidal zone (data not shown) increased during the studied period, suggesting that the microbiota switch from low to high specificity in order to be able to metabolized the nutrients newly available.
In addition to the temporal heterogeneity, a tidal level effect has also been observed for both DG and EF, but their dynamics were very different. The greater effect of the intertidal level on EF microbiota in T2 and T3 may be explained by a greater heterogeneity of environmental conditions in L56 since the immersion period was shorter. Microbiota from bivalve fluids seemed to be particularly sensitive to the sampling zone as it was shown for hemolymph's microbiota from the Pacific oyster [47]. These results seem to hint that EF microbiota are more impacted by environmental variations when they are placed in the upper part of the foreshore. On the contrary, DG microbiota were more impacted by the period when they were placed in the lower part of the foreshore (L20), which was under a greater influence of seawater and feeding sources. So, depending on the considered tissue, microbiota dynamics vary with intertidal position.

GD Are More Stringent Environmental Filters than EF
Analysis of β-diversity components (Supplementary Information 1) also revealed distinct variations between the two-tissue microbiota. EF bacterial microbiota differences between individuals were mostly driven by species replacement for the all three-sampling periods. It is thus hypothesized that EF reached the maximum niche occupancy. On the other hand, DG bacterial microbiota β-diversity was mostly driven by richness difference, except for clams sampled in T1. These results highlight the possibility that OTUs richness might be more constrained in EF and less in DG.
In addition to the fact that α-diversity was higher in the EF microbiota, the depuration process that affected more EF than GD may illustrate that EF microbiota was composed of more transient microorganisms coming from the surrounding seawater. The analysis of β-diversity revealed a higher dispersion in DG microbiota compared to EF microbiota, illustrating a higher within-tissue diversity in the DG. All these results revealed that DG is a potentially more stringent environmental filter than EF.
In conclusion, this study has shown that environmental variations may have a strong influence on clam microbiota, through impact on both digestive gland and extrapallial fluids associated bacteria. DG and EF presented their own microbiota specificity in terms of relative abundance, richness, and dynamics. Furthermore, clam's tissue could be considered as distinct habitats from the bacterial OTUs point of view, as they presented different responses to either small-scale temporal and spatial variabilities or to a depuration process, illustrating that their selection processes differ. These differences underlined that DG could be a more stringent microhabitat probably linked to its particular physicochemical conditions. This study helps to have a better understanding of the Manila clam's microbiota dynamic in the face of environmental variation and especially of the EF microbiota which is still very poorly known. through the postdoctoral fellowship of Clement Offret. This work was also supported by the HORIZON2020 project "Preventing and mitigating farmed bivalve disease" VIVALDI (grant number 678589).

Data Availability
The sequence data have been deposited in the NCBI database under BioProject ID PRJNA748558.