The experiment was designed to examine the effects of normal- and low-salinity, and alternating salinity stress on the gut microbiota in razor clam (S. constricta) in a fully crossed two-factorial design for duration of 96 h. Adult razor clams with average shell length of 5.12 ± 0.31 cm and average body weight of 6.68 ± 0.56 g were sampled from a shellfish aquaculture farm located at Ningbo city, Zhejiang Province, China. After transportation to the lab, the clams were respectively acclimatized under two salinity concentrations (i.e., 20 and 5 psu) for 7 days in filtered flowing seawater (temperature at 22.35 ± 0.21°C and dissolved oxygen at 5.60 ± 0.34 mg/L by constant aeration) with 50% daily water exchange, and fed once daily with microalgae Chaeroeeros moelleri containing rich nutritional composition such as carbohydrate, fat and protein. This microalgae was gained from the Marine Biotechnology Laboratory of Ningbo University. The clams were kept singly in aerated glass jars after being randomly divided into experimental treatments. The clams were either left at their acclimation salinity or stressed to the opposite salinity, resulting in four salinity treatments: normal salinity-acclimated (20 psu, selected as control group, abbreviated as NN), low salinity-acclimated (5 psu, selected as control group, abbreviated as LL), low salinity-stressed (transferred from normal salinity to low salinity, abbreviated as NL) and normal salinity-stressed (transferred from low salinity to normal salinity, abbreviated as LN) ( Figure 1A ). Each group was set with three replicate tanks. In order to monitor the responses of gut microbiota to four salinity treatments, four clams from each tank from each group were sampled at 0, 24, 48, 72 and 96 hours post stress (hps), respectively ( Figure 1A ). Correspondingly, the whole guts from every two clams were dissected aseptically and pooled to form one biological sample due to that inadequate DNA could be obtained from a single clam in trial runs. Ultimately, 120 samples (5 samplings × 6 replicates × 4 treatments) were collected for gut microbial analysis.

The total gDNA (genomic DNA) of clam gut microbiota was extracted using the FastDNA™ SPIN Kit (MO BIO Laboratories) following the instruction manual. The concentrations of gDNA extracts were quantified, and the purity of gDNA was examined using the NanoDrop ND-2000 spectrophotometer. We amplified the V3-V4 regions of gut bacterial 16S rRNA genes with bacterium-specific primer pair: 338F and 806R. Each sample was amplified in 50 μL reaction system in triplicate to minimize PCR-induced bias. For each sample, the PCR cycle condition was carried out based on the description in our previous study (Dai et al., 2023). Further, PCR products were visualized in a 1.5% agarose gel to validate the expected size. Then, we purified and quantified the PCR product from each sample with the commercial PCR fragment purification kit and Quant-iT PicoGreen dsDNA assay, respectively. Each sample was pooled in equimolar amount, and the final amplicons were sent to sequence on Illumina MiSeq platform.

After sequencing, the raw data was processed with an integrated Dix-seq pipeline (Wei et al., 2020). In detail, the pipeline performed the Trimmomatic to filter the barcode and low-quality sequences. After filtering, the atlas-utilis was employed to remove chimeric reads. Paired-end reads were merged with default parameters by performing USEARCH (Edgar, 2010). Finally, the bacterial phylotypes were discerned with UCLUST and binned into operational taxonomic unit (OTU) table at 97% similarity threshold (Edgar, 2010). After taxonomical assignment, the OTUs attached to archaea, eukaryota, chloroplasts, erroneous sequences and those unassigned at the bacteria domain were ruled out from the dataset. We identified the representative sequences and taxonomy classification for each OTU via RDP Classifier in the SILVA database, and constructed a phylogenetic tree of 16S rRNA genes from the filtered alignment with FastTree and IQ-TREE (Price et al., 2010; Lan et al., 2012). Finally, the α-diversity indices of clam gut microbial samples were generated, and the bacterial β-diversity distance between clam gut microbial samples was evaluated based on the OTU clustering results.

We constructed a Sloan’s neutral community model (NCM) to estimate the contribution of stochasticity governing the clam gut microbial community structure (Sloan et al., 2006). The nonlinear least-squares between the OTUs frequency and their abundances were assessed via the Östman’s approach, in which R² determines the overall fit to the neutral community model (Chen et al., 2019). A R² value is close to 1, indicating that the community assembly is fully stochastic.

To assess the effect of alternating salinity stress on the interspecies interactions within clam gut microbiota, bacterial network analyses were performed with MetaMIS (Metagenomic Microbial Interaction Simulator) software by introducing time-series dataset of gut bacterial compositions at OTU level for NN, NL, LN and LL group samples, respectively. MetaMIS is a tool based on the generalized Lotka Volterra equations, and is designed to infer underlying microbial interactions based on the metagenomic abundance profiles (Shaw et al., 2016). Positive or negative interactions between two nodes (i.e., OTUs) in each microbial interaction network were inferred. This approach systematically detects the interaction patterns of each node with other nodes, such as mutualism (+/+), competition (−/−), parasitism or predation (+/−), commensalism (+/0), amensalism (−/0), and no effect (0/0) in four salinity stress groups, respectively.

Function potentials of clam gut bacterial communities were inferred using PICRUSt2 v2.4.1 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) according to the KEGG (Kyoto Encyclopedia of Genes and Genomes) database for phenotype prediction (Douglas et al., 2019). Significantly differential functional pathways between NN96 and NN0, between NL96 and NL0, between LN96 and LN0 and between LL96 and LL0 were selected using response ratio analysis at a confidence interval (CI) of 95% (Dai et al., 2022a).

One-way analysis of variance (ANOVA) was used to detect the significant differences in the clam gut bacterial compositions at family level among four salinity stress groups at same stress time (Churchill, 2004). Based on Bray-Curtis distance, a principal component analysis (PCA) and analysis of similarity (ANOSIM) were performed to assess the differences in the β-diversity of clam gut bacterial communities among four salinity stress groups with the increase of stress time (Mishra et al., 2017). A permutational multivariate analysis of variance (perMANOVA) was carried out to qualify the contributions of salinity stress, hps, and their interaction on the variation in clam gut bacterial communities based on “vegan” package in R (R Development Core Team, 2013). Venn diagram was applied to count the number of shared and unique OTUs among four salinity stress groups at same time point and at each salinity stress group over stress time.

The survival rate of S. constricta was unaffected by the 96-hour salinity stress in NN and LL groups ( Figure 1B ). In contrast, the mortality was firstly observed on 24 hps and showed a continuously decreased trend in NL and LN groups. Notably, the cumulative survival rate in NL group was 78% on 96 hps, which was significantly lower than that in LN group (93%) ( Figure 1B ).

A total of 17,689,879 high-quality sequences were gained with an average of 147,416 ± 55,299 (mean ± standard deviation) sequences per sample, which generated 1784 OTUs across our enrolled 120 samples corresponding to NN, NL, LN and LL groups. Among these OTUs, there were 35, 42, 52, 34 and 56 OTUs shared among four salinity stress groups at 0, 24, 48, 72 and 96 hps, respectively, while unique OTU number varied accordingly in each salinity stress group ( Supplementary Figure S1A ). In addition, the NN and LL groups displayed more shared OTUs than the NL and LN groups during 96-hour salinity stress period ( Supplementary Figure S1B ). The trend of unique OTU numbers was ordered as NN48 (125) > NN96 (88) > NN24 (87) > NN0 (76) > NN72 (66), NL24 (178) > NL48 (155) > NL96 (102) > NL0 (85) > NL72 (77), LN72 (155) > LN96 (136) >LN24 (98) > LN0 (95) > LN48 (32), LL72 (165) > LL24 (132) > LL48 (114) > LL0 (75) > LL96 (65) within 96 hps ( Supplementary Figure S1B ).

Alternating salinity stress affected the α-diversity of clam gut bacterial communities. The α-diversity index such as observed species in NN and LN groups increased at 24 hps and then dropped continuously with lowest level at 96 hps ( Figure 1C ). However, an opposite trend was observed in NL and LL groups at 24 hps. Notably, the observed species was lower in NL and LL groups at 96 hps than that at 0 hps ( Figure 1C ). The Shannon diversity in NL and LL groups kept an increased trend from 0 to 48 hps and then decreased at 72 and 96 hps ( Figure 1D ). A completely opposite pattern was detected in LN group within 96 hps. For NN group, the Shannon diversity firstly raised at 24 hps but declined at 48 hps, then increased at 72 hps and finally decreased at 96 hps ( Figure 1D ).

The bacterial phyla were predominant by Tenericutes (20.4% ± 12.7%, mean ± standard deviation) and Gammaproteobacteria (16.1% ± 7.0%), followed by Alphaproteobacteria (14.1% ± 13.2%), Campilobacterota (12.5% ± 5.8%), Bacteroidetes (5.3% ± 4.9%), Actinobacteria (4.6% ± 3.3%), and Firmicutes (4.1% ± 5.4%) across the samples. Additionally, Betaproteobacteria (0.5% ± 0.3%), Fusobacteria (0.3% ± 0.7%), Planctomycetes (0.3% ± 0.4%), and Verrucomicrobia (0.1% ± 0.1%) were identified across the samples at low abundances ( Supplementary Figure S2A ). These large standard deviations might be due to the highly temporal succession of gut bacterial compositions ( Figure 2 ). The relative abundance of Tenericutes decreased significantly (P < 0.05) in four salinity stress groups at 24, 72 and 96 hps compared with corresponding 0 hps. In contrast, the relative abundances of Alphaproteobacteria and Gammaproteobacteria drastically increased (P < 0.05) at 24 hps but decreased at 96 hps in the four salinity stress groups. Notably, the relative abundances of Firmicutes, Actinobacteria, Fusobacteria, Planctomycetes, Betaproteobacteria and Verrucomicrobia exhibited consistent declined trend in four salinity stress groups at 96 hps ( Supplementary Figure S2A ). At finer family level, the dominant families were Arcobacteraceae (12.4% ± 5.8%), Vibrionaceae (5.5% ± 4.3%), Pseudomonadaceae (5.2% ± 3.8%), Sphingomonadaceae (4.7% ± 6.8%), Methylobacteriaceae (4.5% ± 6.7%), Shewanellaceae (2.3% ± 3.4%), Flavobacteriaceae (2.3% ± 1.9%) and Rhodobacteraceae (2.2% ± 1.8%) ( Supplementary Figure S2B ). However, alternating salinity stress led to an apparent divergence in dominant taxa distribution in clam gut within 96 hps ( Figure 3 ). It is worth noting that the relative abundances of potentially pathogenic Pseudomonadaceae and Shewanellaceae increased in LN and LL groups at 96 hps relative to that at 72 hps. An opposite trend for this phenomenon was observed in NN and NL groups ( Figure 3 ).

The PCA plot illustrated dissimilarities in varying degree in the gut microbiota among four salinity treatments and among different stress time ( Figure 2A ). For instance, overall, the four salinity treatments markedly contributed to the divergence of the gut bacterial community at each time point ( Figures 2B–F ). At 0 hps, a significant difference in gut bacterial community was observed only between NL and LN groups. In the NN group, apparent separation was observed between each pair of groups among stress time, with an exception between NN24 and NN48 and between NN72 and NN96. For the NL group, the gut bacterial structures were dramatically different between each pair of groups within 96 hps, with an exception between NL72 and NL96. For LN group, the gut microbiota in LN0 was clustered away from that in LN24, LN72 and LN96, and each pair of the latter three groups exhibited apparent separation. For LL group, the gut microbiota were remarkably distinct between each pair of groups within 96 hps, with an exception between LL24 and LL0/LL48 ( Figure 2 ). These findings were further confirmed by a dissimilarity test ( Supplementary Table S1 ). Further, PerMANOVA revealed that salinity stress, hours post stress and their interaction dramatically contributed to 4.5%, 9.5% and 6.0% variation in the gut bacterial communities, respectively (P < 0.01 for the three variables) ( Table 1 ).

The R² value was 0.279 in NN group, which increased to 0.491 and 0.448 in LN and LL groups, respectively, but was comparable with that NL group (0.270) on 0 hps ( Figure 4 ). Therefore, low salinity stress enhanced stochasticity underlying the gut microbiota. The subsequent alternating salinity stress aggravated stochasticity. For example, the R² value in NL (R² = 0.382) group at 96 hps was higher than that at corresponding 0 hps ( Figures 4B, D ). However, the ecological processes did not change in NN group at 96 hps due to the R² value (0.270) was comparable with that at 0 hps ( Figure 4A ).

To estimate the effects of alternating salinity stress on the species-to-species interactions within gut microbiotas, bacterial network analyses based on time-series dataset of gut microbial compositions were performed with MetaMIS for the samples in NN, NL, LN and LL groups, respectively ( Figure 5 ). We exhibited the gut bacterial interaction relationships among core taxa in four salinity stress groups. Compared to that in the NN group (66.7%), the positive interactions decreased in NL (48.1%), LN (59.0%) and LL (51.9%) groups, suggesting a lower degree of cooperative activities in the gut microbiota of latter three groups. In contrast, interspecies interaction networks presented more complex and higher connected topology in NL, LN and LL groups compared to NN group ( Figure 5 ). Notably, there were more number of core OTUs that were attached to Vibrio, Pseudomonas, and Shewanella taxa (core ratio = 100%, Supplementary Table S2 ) in the gut network of LL group in relation to NN group ( Supplementary Table S3 ). In the four microbial networks, a variety of interaction types were observed between core OTUs. These included competition (-/-), parasitism/predation (+/-), commensalism (+/0), and amensalism (-/0) ( Supplementary Figure S3 ). For instance, the most frequent interactive relationship for OTU3 was mutualism (+/+) in NN group, which was shifted to parasitism or predation (+/-) in NL, LN and LL groups. For OTU5, the most frequent interaction was competition (−/−) in NL and LL groups, but was mutualism (+/+) in NN group and parasitism or predation (+/-) in LN group ( Supplementary Figure S3 ). These results indicated that alternating salinity stress affected the interspecies interaction relationships of gut microbial networks.

To evaluate how the functional potentials of gut microbiota were influenced by alternating salinity stress, the KEGG pathways were compared among four salinity stress groups ( Figure 6 ). We showed the functional pathways that prominently varied between each two pair groups in four salinity stress groups at 0 and 96 hps. Notably, the functional potentials exhibited an opposite pattern in two alternating salinity stress groups (i.e., NL and LN). For instance, the pathways related with the metabolisms of amino acid (lysine biosynthesis), carbohydrate (amino sugar and nucleotide sugar metabolism, pentose phosphate pathway), energy (carbon fixation in photosynthetic organisms, photosynthesis), lipid (glycerolipid metabolism), and nucleotide (purine metabolism, pyrimidine metabolism), replication and repair (DNA replication), and immune system (NOD-like receptor signaling pathway) dramatically decreased (P < 0.05) in NL group at 96 hps compared to corresponding 0 hps, whereas these pathways performed opposite tendency in LN group. Similar findings for this phenomenon were also observed in the two control groups (i.e., LL and NN) ( Figure 6 ).