Introduction

With the development of modern industrialization and urbanization, global warming has made a thoroughly grave threat to the all-rounded hierarchical biosphere, from marine ecosystems to terrestrial biosphere1,2,3. Rather than directly leading to ecosystem collapse or species extinction, warming tends to gradually reconstruct the eco-structure and bring biodiversity variations due to abundance alterations4,5. Along with the deeper evaluations in warming to different bio-hierarchies, compared to mammals, ectothermic animals are more sensitive to temperature changes as warming inevitably affects their distribution, behavior, and health6,7,8. Indeed, through chronic temperature increases and, more frequently, extreme heatwave events, global warming mediates the healthy status of ectothermic animals due to temperature effects on physiology, metabolism, and immune capability9,10,11. In particular, animal’s health is also tightly connected to the stability of intestinal microbiota (IM), not only since microbe is deem to a sensitive ecological indicator for temperature changes12,13, but also because intestinal microbiota plays essential roles in the health of their hosts due to extensive metabolic and immune interactions14,15,16.

The host-IM relationship is considered as one of the most intricate interactive relationship17. Considerable variation in IM has been described among individuals, therefore, identifying compositional patterns across large populations and geographies potentiates microbiota-based diagnostics or prevention of disease18. Hence, defining a large number of microbial compositions as only 2–4 types with distinctive compositions and functions that are dominated by different keystone genera, which is also called “enterotypes”, is an effective strategy to form a unified landscape of IM to a specific host18,19,20,21. Based on network analysis, different keystone genera show potential competing relationship19,21, which brought convenience to evaluate microbial variations among enterotypes. More broadly, enterotypes revealed a ratio of Bacteroides and Prevotella as a healthy marker aligned with inflammation and dietary in human18,22, yet direct transitions were rare between Bacteroides-dominated and Prevotella-dominated communities, suggesting the presence of a transition barrier between two enterotypes22. This relative stability is mainly attributed to the well-regulated host immunity and long-term dietary interventions23,24. However, these phenomena are usually reported in mammals rather than in ectothermic animals, such as aquaculture animals like fish and shrimp25,26. The IMs of ectotherm tend to be much more dynamic and less diverse than those of endotherm due to the lack of enough strong deterministic factors in aquaculture animals, such as habitat, environmental parameter, and geographic origin14,27,28,29.

Aquaculture products, also called “Blue Food”30,31, are conducive to tackle micronutrient and protein deficiencies in the world, which are therefore considered to potentially satisfy global food and nutrition demands32. Pacific white shrimp, Litopenaeus vannamei, is one of the most successful models of worldwide utilizing marine bioresources in more than 38 countries and areas, with global production of 5.6 million tons and valued at over $28 billion in 202233. At present, shrimp industry is currently being threatened by several bacterial diseases, such as acute hepatopancreatic necrosis disease (AHPND), hepatopancreas necrosis syndrome (HPNS) and white feces syndrome (WFS)28,34,35,36. The outbreaks of WFS and HPNS are proved to be attributed to shrimp IM dysbiosis and multiple bacterial markers were further identified28,34,35,37. Notably, the probability of WFS occurrence is significantly increasing during high temperatures, reaching 33–34 °C in pond water34,38. Indeed, WFS, once a tough problem in many low-latitude Asian countries, has spread to high-latitude areas in recent years through chronic temperature increases34,39,40. In parallel, shrimp mobilizes its immunity to pathogens during temperature fluctuation, leading to intestinal microbiota alterations41,42. As observed in aquaculture practices, both WFS and HPNS can be observed in shrimp ponds in 1–3 days after dramatical temperature change. In this sense, shrimp health is facing disease challenges along with the more frequent ambient warming observed in aquaculture practices.

Aquaculture is responsible for the continuing impressive growth in the supply of protein consumption. As a model of aquaculture ectothermic animals, shrimp is expected to be particularly vulnerable to temperature changes. However, it remains poorly understood how ambient warming mediates shrimp health. Herein, two hypotheses are proposed: (i) shrimp IM can be stratified into different enterotypes, which are relevant to shrimp health status; (ii) ambient warming shapes shrimp microbial structure and enterotype distribution. To address these hypotheses, a total of 1,369 shrimp IM data was obtained from major representative shrimp culture regions in nine countries. Based on the multi-omics analysis, shrimp IM was stratified into three enterotypes that were closely relevant to health statuses. Moreover, temperature was the most shaping factor for microbial structure and led to the variation of probability of disease, which was further validated through ambient warming and antimicrobial peptide gene knockdown experiments. Collectively, these findings provided a globally comprehensive catalog for shrimp intestinal microbiome, and revealed that warming-driven enterotypes mediated host health, which helped to improve shrimp health culture management from a microecological perspective and alerted the inevitable challenge of global warming in aquaculture field.

Results

Catalog of the global shrimp intestinal microbiome

To investigate the intestinal microbial characteristics of Pacific white shrimp as well as the relationship between IM and ambient warming, a total of 15 cohorts from nine countries (Brazil, China, Ecuador, Indonesia, Malaysia, Mexico, Thailand, USA and Vietnam) which accounted for 79.8% global production (Supplementary Table 1)35,43,44,45,46,47,48,49,50,51, comprising 1369 IM data were enrolled in this study (Table 1). Water physicochemical parameters (temperature, salinity, dissolved oxygen (DO) and pH) of 1150 samples among these 1369 data collected by our team were determined in field (Table 1). This large dataset covered the major shrimp aquaculture areas in the world, including different ages (15–91 day-post-hatching), aquaculture patterns (earthen ponds, indoor farming system and Aquamimicry mode) (Supplementary Table 2), which were fully representative for the global shrimp aquaculture conditions.

Table 1 Summary of the shrimp intestinal 16S rRNA amplicon sequencing data used in the present study

After performing 16S rRNA gene amplicon sequencing and quality control, a total of 61 million clean reads were obtained, resulting in 13,204 amplicon sequence variants (ASVs). Rarefaction analysis showed that the ASV number reached a plateau as the sample number increased, and Shannon index and Chao1 index showed a similar trend (Supplementary Fig. 1a), suggesting that the whole intestinal composition at a global scale was well represented in this study. The ASVs were further identified into 63 phyla (Supplementary Fig. 1b), which showed Proteobacteria, Tenericutes, and Bacteroidetes were the most abundant phyla at global scale (Fig. 1a, Supplementary Fig. 1c). The most relative abundant genera were Vibrio (30.8%), Candidatus Bacilloplasma (16.5%), Shewanella (8.6%) and Photobacterium (7.0%) (Fig. 1b). The alpha diversity of shrimp IM was further investigated to assess the complexity of microbial communities by Chao1 index and Shannon index. In global scale, the observed ASVs number for one sample ranged for 304 to 2148, and the Shannon index ranged from 0.121 to 9.814. The lowest alpha diversity of shrimp intestinal microbiota was observed in China, which the highest one was observed in Thailand, which indicated that the alpha diversity varied among countries (Fig. 1c). Apparent latitudinal pattern was detected in the observed ASV number, Chao1 index and Shannon index, as the richness and diversity peaked at low latitude (Fig. 1d).

Fig. 1: Catalog of global shrimp intestinal microbiota.
figure 1

a A total of 1369 shrimp intestinal data were collected from nine countries to profile the microbial composition of pacific white shrimp. The relative abundance of bacterial communities at phylum level. b The relative abundance of bacterial communities at genus level. c The alpha diversity of shrimp IM at different countries indicating by Chao1 index and Shannon index. d Latitude distribution of microbial diversity, indicating by observed ASV number, Chao1 index and Shannon index with second order polynomial regression lines, as diversity peaked at low latitude. e A total of 19 core ASVs were obtained. The total relative abundance of core ASVs were present in 52.4% at shrimp IM. f Fifteen of the 19 ASVs belonged to Proteobacteria, three ASVs were Tenericutes and one ASV was Firmicutes. The heatmap showed the relative abundance (normalized by Z score) of core ASVs at different countries.

Core shrimp intestine microbes at global scale

Following the previous study, the core ASVs for shrimp IM were determined based on the abundance (>0.01%) and occurrence frequency in all samples (>80%) according to the definition of the core microbes52, and a total of 19 ASVs were identified as the core ASVs for the global shrimp IM (Supplementary Fig. 1d). The total relative abundance of the core ASVs was present at 52.4% in shrimp IM (Fig. 1e). Fifteen of the 19 ASVs belonged to the phylum Proteobacteria, three ASVs were Tenericutes, and one was belonging to Bacteroidetes (Fig. 1f). These core ASVs are reported to play important roles in shrimp health, for example, being relevant to WFS or AHPND outbreaks (Table 2). yet these core ASVs showed a great abundance variation at different countries (Fig. 1f).

Table 2 Advances of the core ASVs in shrimp

Shrimp IM was stratified into three enterotypes

To evaluate the enterotype in Pacific white shrimp, a Dirichlet multinomial mixtures (DMM)-based enterotype analysis was conducted to investigate the presence and characteristics of shrimp enterotypes based on the underlying compositional structure. Eventually, three enterotypes were identified based on the Dirichlet scores (Supplementary Fig. 2), which were dominated by a relatively high abundance of Vibrio (ET V, n = 496), Shewanella (ET S, n = 455), and Candidatus Bacilloplasma (ET CB, n = 418), respectively (Fig. 2a, b). It was evident to find that the relative abundance of Vibrio was positively correlated with PCoA dimension 1, while Candidatus Bacilloplasma was negatively linked to dimension 1, and Shewanella showed a high abundance at the case that both Vibrio and Candidatus Bacilloplasma were with a relatively low abundance, suggesting that the abundance of three keystone genera is negatively correlated with each other. In parallel, the stratification of enterotype was strongly supported by their functional differences. A metagenomic sequencing of a subset of 131 IM samples (Supplementary Table 3) from six cohorts, containing 1291.65 GB data, was further conducted to evaluate the microbial functional differences among enterotypes. The metagenomic sequencing data was filtered with shrimp host genome to obtain the effective microbial data (45.6–80.3%, Supplementary Table 3), and was further used to conduct the functional annotation, which shown that the microbial functions of three enterotypes were distinctive from each other (MRPP, P = 0.001, Fig. 2c).

Fig. 2: Enterotype stratification for shrimp.
figure 2

a A DMM-based enterotype analysis was conducted to identify three enterotypes. PCoA showed that the microbial structure of three enterotypes was separated. Relative abundance of keystone genera, Vibrio, Shewanella and Candidatus Bacilloplasma were related to dimension 1 and 2. b Abundance distribution of the three keystone genera at different enterotypes. c A metagenomic sequencing of a subset of 131 IM samples was further conducted. PCoA showed that the potential functions of three enterotypes were distinctive from each other (PerMANOVA, P = 0.005**). d The alpha diversity indicated by Chao1 index and Shannon index among three enterotypes were significantly different (ANOVA, P < 0.05*). e Assembly patterns of different ecological processes for three enterotypes. f Spearmans correlation to determine the relationship among keystone genera and other microbes. g The distribution of enterotypes among different countries.

Based on the enterotype stratification, twelve genera, such as Vibrio, Pseudomonas, Exiguobacterium, and Rhodopirellula were abundant in ET V, which was enriched with functions related to linolenic lipid metabolism, secondary bile acid biosynthesis and carbohydrate digestion and absorption (ANOVA, P < 0.05, Supplementary Fig. 3). Eight genera, such as Shewanella, Stenotrophomonas and Gemmobacter, as well as functions including mineral absorption, N-Glycan biosynthesis, and polycyclic aromatic hydrocarbon degradation in ET S. Sixteen genera, such as Candidatus Bacilloplasma, Photobacterium and Aeromonas, were overrepresented in ET CB, with high abundance functions of fatty acid biosynthesis, fatty acid metabolism, and galactose metabolism.

The difference of intestinal microbial diversity and structure among the three enterotypes was further investigated. The alpha diversity among the three enterotypes was significantly different (ANOVA, Chao1 index: P = 0.005; Shannon index: P = 0.001 Fig. 2d). The relative abundance of Candidatus Bacilloplasma and Shewanella were positively corelated to Shannon index and Chao1 index (P = 0.001, R2 = 0.269 and P = 0.021, R2 = 0.184, respectively), yet Vibrio showed an opposite trend (P = 0.002, R2 = 0.311, Supplementary Fig. 4). The standardized effect size of the mean nearest taxon distance measure was used to further determine the processes governed the assembly of microbiota in shrimp intestine, which showed that the deterministic process weakened in ET CB compared to ET V and ET S (Fig. 2e). The relationship between keystone genera and other microbes were further determined by Spearman’s correlation analysis. Vibrio was positively correlated to 24 genera but negatively correlated to 131 genera in ET V (Fig. 2f), indicating a competitive relationship with other microbes. ET S contained 101 positive correlations and 50 negative correlations to Shewanella, while the number was 132 positive and 40 negative correlations to Candidatus Bacilloplasma in ET CB (Fig. 2f). Furthermore, the regional distribution of enterotype among countries showed a great variation, while ET CB seemed to be more distributed in Thailand and ET S hold a large distribution in Vietnam (Fig. 2g).

Enterotype was relevant to shrimp health status

To evaluate the relationship between shrimp enterotype and host health, a subset including 165 diseased shrimp (98 WFS shrimp and 67 HPNS shrimp with clinical and histological symptom, Supplementary Fig. 5) and 167 control shrimp (collected from the same pond and without any diseased syndrome) were chosen (Supplementary Table 4). For the 98 WFS shrimp, 31 (32%), 4 (4%), and 63 (64%) shrimp was belonging to ET V, ET S and ET CB (Fig. 3a). For the 67 HPNS shrimp, 40 (60%),10 (15%) and 17 (25%) shrimp was belonging to ET V, ET S and ET CB. Meanwhile, for the shrimp without any disease signs, 53 (32%), 45 (27%), and 69 (41%) shrimp was belonged to ET V, ET S, and ET CB, which showed an obvious enterotype distribution bias for shrimp diseases. Furthermore, a random forest model discriminated healthy and diseased shrimp was conducted to predict the probability of disease index (POD) in shrimp. Based on a ten-fold cross-validation error curve, a total of 159 and 96 top-ranking disease-discriminatory ASV were identified for WFS and HPNS (Supplementary Fig. 6a and d). The random forest model was applied to all shrimp IM data, which showed that the POD values for WFS in ET CB were the highest, while POD values for HPNS were significantly higher in ET V than those in ET S and ET CB (ANOVA, P = 0.001, Fig. 3b). Interestingly, most of these diseased-biomarker ASVs were also significantly varied among three enterotypes (Supplementary Fig. 7), which also indicated that the diseased-related metacommunity was associated to shrimp enterotypes.

Fig. 3: Enterotypes were relevant to shrimp health.
figure 3

a A subset including 165 diseased shrimp (98 WFS shrimp and 67 HPNS shrimp with clinical and histological symptom) was chosen to evaluate the relationship between enterotype and host health. b A random forests model was conducted to predict the probability of disease for WFS or HPNS. POD values were significantly different among three enterotypes (ANOVA, P < 0.01**). c The molecular ecology networks of different enterotypes, which showed that ET V and ET S hold higher average connectivity and lower average path lengths than ET CB. d Relationship between POD and abundance of the three keystone genera. POD value was significantly correlated with ratio of Vibrio/Candidatus Bacilloplasma, while POD value decreased as Shewanella increased. e The ratio Vibrio/Candidatus Bacilloplasma could effectively discriminate shrimp enterotype at dimension 1 of PCoA. f The ratio Vibrio/Candidatus Bacilloplasma was negatively correlated to Chao1 index and Shannon index.

Moreover, considering that the stability and complexity of the intestinal microecosystem were reported to be linked to shrimp health and disease occurrence35,53, a molecular ecological network analysis (MENA) was constructed to understand microbial interaction in the whole intestinal community. Compared to ET V and ET S, ET CB showed the lowest average connectivity and highest average path lengths (Fig. 3c), indicating that species interactions in ET CB were less cooperative and less complex than ET V and ET S. Furthermore, the POD value was significantly correlated with ratio of Vibrio and Candidatus Bacilloplasma (WFS: P = 0.001, R2 = 0.213; HPNS: P = 0.001, R2 = 0.175), while POD value decreased as Shewanella increased (WFS: P = 0.001, R2 = 0.186; HPNS: P = 0.020, R2 = 0.124, Fig. 3d), which indicated that the increases of Vibrio or Candidatus Bacilloplasma contributed to probability of shrimp disease occurrence.

Overall, in consideration of the abundance of the Candidatus Bacilloplasma, Vibrio and Shewanella, as well as their governed ecological processes, a brief scheme was purposed to distinguish shrimp enterotype: (I) Vibrio/Candidatus Bacilloplasma > 1: ET V; (II) high abundance of Shewanella and Vibrio/Candidatus Bacilloplasma ≈ 1: ET S; (III) Vibrio/Candidatus Bacilloplasma < 1: ET CB (Fig. 3e). This scheme highlighted the significance of the ratio of Vibrio and Candidatus Bacilloplasma in distinguishing shrimp microbial diversity (Fig. 3f), enterotype and health statuses (Figs. 3a, 3b and 3d, Supplementary Data 1 and 2). These findings suggested that shrimp enterotype was tightly relevant to host health.

Temperature was the most shaping factor for shrimp IM

To evaluate the relationship between IM and environmental factors, a variance partitioning analysis (VPA) was performed to quantify the relative contributions made by environmental parameters and geographic distance to the microbial structure of IM. The water parameters explained 26.9% of the observed variation, while geographic distance explained 5.8%, and the interaction effect was 8.5% (Fig. 4a). Mantel test revealed that the microbial structure exhibited a strong correlation with environmental water parameters (r2 = 0.486, Table 2). Further partial Mantel tests were performed to assess the associations between water variables and microbial structure. Temperature held the highest correlation to microbial structure variation (Table 2). Particularly, the correlation with temperature (rM = 0.444) was more than three times as strong as with salinity (rM = 0.141). To further validate this result, a total of 4 PCoA plots were constructed to observe the influence of temperature, pH, salinity, and dissolved oxygen, followed by PerMANOVA to test the significance difference among groups, which showed that different levels of temperature, salinity and pH exerts a profound influence on microbial composition (Supplementary Fig. 8). Similarly, aggregated boosted tree (ABT) models demonstrated that temperature was the primary factor affecting both phylogenetic diversity and phylotypes, accounting for nearly 27.8% of the relative influence for phylogenetic diversity and 29.6% for phylotype (Supplementary Fig. 9a). The structural equation modeling (SEM) further partitioned the direct and indirect effects of multiple factors on microbial diversity and community structures, which also suggested that temperature exerted the largest influence to microbial community structure (Fig. 4b).

Fig. 4: Temperature was the most important factor shaping intestinal microbial structure.
figure 4

a Variance partitioning analysis was used to quantify the relative contributions made by environmental parameters and geographic distance to microbial structure of IM. b Structural equation modeling was conducted to partition the direct and indirect effects of multiple factors on microbial diversity and community structures. c Enterotype distribution at different ambient temperature. d Relationship between temperature and alpha diversity. As temperature rising, the Chao1 index and Shannon index increased. e Spearman’s correlation analysis of temperature and the three keystone genera.

To deeper understand the temperature effect on shrimp enterotypes, the distribution of enterotypes at different temperatures was evaluated, which showed that ET V and ET S distributed at lower temperatures than ET CB (Fig. 4c). Indeed, temperature exerts the highest effect on microbial alpha diversity, demonstrated by alpha diversity increased as the temperature raised (Chao1 index: P = 0.011, R2 = 0.391; Shannon index: P = 0.039, R2 = 0.199, Fig. 4d). Additionally, Spearman’s correlation analysis demonstrated that 14 of the 19 core ASVs were significantly correlated to temperature variation (P < 0.05, Supplementary Fig. 9b). The relative abundance of Vibrio was negatively correlated to temperature (P = 0.023, R2 = 0.155, Fig. 4e), while the relative abundance of Shewanella and Candidatus Bacilloplasma increased under relatively higher temperature (Shewanella: P = 0.040, R2 = 0.103; Candidatus Bacilloplasma: P = 0.008, R2 = 0.258). Overall, the above results indicated that temperature played a primary role in shaping shrimp microbial composition.

Experimental ambient warming altered shrimp enterotypes and antimicrobial peptide expressions

To confirm and assess the effect of temperature on shrimp IM composition, shrimp were cultured in 20 °C, 24 °C, 28 °C, 32 °C and 36 °C ambient for one week (named groups A, B, C, D, and E, respectively, n = 60 for each group) (Fig. 5a). However, in group E (36 °C), shrimp’s food intake efficiency significantly decreased, and worse still, over 50% shrimp died within 4 days. As temperature decreased, the relative abundance of Vibrio tended to decrease (P = 0.015, R2 = 0.318), while Shewanella and Candidatus Bacilloplasma showed an opposite trend (Shewanella: P = 0.031, R2 = 0.179; Candidatus Bacilloplasma, P = 0.003, R2 = 0.243, Fig. 5b), which led to a dropping ratio of Vibrio/Candidatus Bacilloplasma (P = 0.001, R2 = 0.396, Fig. 5c). More broadly, the alpha diversity gradually increased in a warming ambient from 20°C to 36 °C, as supported by Chao1 index and Shannon index (Chao1: P = 0.001, R2 = 0.356; Shannon: P = 0.017, R2 = 0.227, Fig. 5d). Additionally, MENA revealed that the 32 °C and 36 °C showed lower average connectivity than those in other groups (Fig. 5e). A PCoA revealed that the microbial structures among different ambient temperatures were markedly distinct from others (PerMANOVA, P = 0.008, Fig. 5f). Notably, the POD values for WFS in 28 °C, 32 °C and 36 °C were significantly higher than those in 20 °C and 24 °C, while the POD values for HPNS showed an opposite trend (Fig. 5g). These above results were consistent with the results from the global data analysis, which confirmed that temperature exerted a profound influence on the microbial composition and structure.

Fig. 5: Experimental warming altered shrimp enterotypes distribution and regulating antimicrobial peptide expressions.
figure 5

a Schematic representation of warming experiment combined with multiomics analysis and RNAi strategy. b Relative abundance of keystone genera at different temperature. c The Ratio of Vibrio and Candidatus Bacilloplasma decreased as ambient temperature raised. d As temperature increased, the Chao1 index and Shannon index tended to be higher. e MENA demonstrated that as temperature raised, the complexity and stability of species interaction network decreased. f PCoA showed that the microbial structure of different temperature was distinctive from each other (PerMANOVA, P = 0.008). g POD values for WFS and HPNS varied significantly at different ambient temperature. The above results were consistent with the results from the global data analysis, which confirmed that temperature exerted a profound influence on the microbial composition and structure. h A transcriptome analysis for shrimp intestine in different temperature was conducted to further explore how ambient warming alter IM composition. PCoA showed that samples from different temperature groups were separated (PerMANOVA, P = 0.044). i Alteration of keystone genera after silencing the expression of Pen2, Pen3, Pen4, Alf1, Alf2, Alf3 and Alf4 with RNAi strategy. The variation of keystone genera was consistent with the antimicrobial peptide expression status under different ambient temperature.

To further explore how ambient warming alter IM composition, a transcriptome analysis for shrimp intestine in different temperature was conducted. After aligning to shrimp genome, a total of 32,499 transcripts were obtained. PCoA was performed to determine the transcript profiles across groups (Fig. 5h), which showed that samples from different temperature groups were separated (PerMANOVA, P = 0.044). Moreover, a Procrustes analysis to show the correlation of microbial composition and transcriptome, which suggested that the gene expression was linked to microbial variation (Supplementary Fig. 10). A total of 2276 differentially expressed genes (DEGs) were identified among five groups by ANOVA (Supplementary Fig. 11). Among these DEGs, there were 32 DEGs whose expression were significantly related to the relative abundance of Vibrio, Shewanella and Candidatus Bacilloplasma. Some DEGs were previously reported to be responsible for antimicrobial activity in shrimp, such as antimicrobial peptides like penaeidins (PENs, LOC113808998, LOC113808997) and antilipopolysaccharide factors (ALFs, LOC113800363, LOC113807189, LOC113810045) (Supplementary Table 5). To assess the relationship between antimicrobial peptides and microbial composition, the expression of Pen2, Pen3, Pen4, Alf1, Alf2, Alf3 and Alf4 was knockdown using RNAi strategy (Supplementary Fig. 12). Compared to the control group (dsGFP), the relative abundance of Vibrio increased after silencing of Pen4, and Alf4 (ANOVA, P = 0.033 Fig. 5i). The relative abundance of Candidatus Bacilloplasma increased in case of silencing Pen3 and Alf2, yet decreased after silencing Pen4. This phenomenon was consistent with the antimicrobial peptide expression status under different ambient temperatures, which validated that ambient warming altered shrimp enterotypes and antimicrobial peptide expressions.

Discussion

Global warming has become a well-recognized concern among the public. Since ectothermic animals occupy a majority in food web, their distributions and healthy statuses are essential for the whole ecosystem. Changes in microbiota may be caused by rising temperatures, leading to potential intestinal dysbiosis and an increase in disease prevalence54. In practical terms, more frequent disease outbreaks were observed in aquatic ectotherms during temperature fluctuation34,55,56. Herein, the present study provided solid evidence for the underlying mechanism that warming was relevant to the probability of shrimp diseases by altering enterotypes, which alerted the unavoidable challenges for aquaculture animals under the context of global warming.

Pacific white shrimp is currently facing devastating diseases relevant to IM dysbiosis. Given the complexity and variation of the microbial ecosystem, forming a fulfilled landscape for shrimp and deeper evaluating the features of IM dysbiosis may help to reveal the etiology for shrimp disease and prevent outbreaks. Moreover, shrimp IM was stratified into three enterotypes, which exhibited distinctive patterns on function, diversity and specie interaction (Fig. 2). Indeed, even mammal enterotypes are usually independent of age, gender, cultural background and geography19, some investigations still potentially indicate that enterotype is associate with different health states22,57,58. A potential explanation for this correlation is that the relative abundance of keystone genera correlated with multiple metabolic and immunologic functions can form thorough intestinal ecological status linking to host health22. An obvious bias in shrimp enterotypes and diseases was observed (Fig. 3). ET V was dominated by Vibrio, which is reported as an opportunistic pathogen and is relevant to IM dysbiosis in shrimp suffering from HPNS, AHPND28,34,35,40. Herein, Vibrio showed a large percentage of negative links to other microbes, which indicated that a high abundance of Vibrio might inhibit the other microbes and, therefore, lead to lower alpha diversity and a higher probability of HPNS. Besides, Candidatus Bacilloplasma showed more positive links to other microbes than those in ET V and ET S (Fig. 2). However, the high percent of positive microbe interaction links is detrimental to form a stable microbial network under external disturbance59. A fewer complex species interactions network and higher POD value for WFS was found in ET CB, which was consistent with previous studies that lower complexity and stability is a typical feature of IM dysbiosis for WFS37,60. The ratio between keystone genera was proposed as an indicator for not only stratification itself but also ecological state, such as ratio of Bacteroides and Prevotella in human18,22. Given the bias between enterotype and disease, ratio of Vibrio and Candidatus Bacilloplasma was proposed, which was significantly related to microbial diversity, enterotypes and diseases (Fig. 3). Hence, enterotypes were crucial indicators for shrimp health, which could be further applied in future aquaculture to identify or manage the healthy status of shrimp.

Understanding the forces that underlie IM’s stability and resilience is the premise for engineering IM to improve host health. Herein, warming increased shrimp intestinal diversity but decreased the complexity and stability of species interaction network, which was consistent with previous studies in other ectothermic animal that climate warming induces alterations in diversity and function, leading possible negative consequences for common lizard’s survival54. Temperature rising was the primary factor to intestinal microbial structure, and exerted a profound influence on alpha diversity (Chao1 index and Shannon index, Fig. 4d). Indeed, the temperature ranges of Thailand (28.0–35.7) and China (21.0–31.9) were broader than those in other countries (Table 1), which probably led to the large range of alpha diversity in this global investigation. Brazil is the only country exhibiting only one of the three possible enterotypes. In this study, only 4 shrimp intestinal data in Brazil were included (PRJNA473127). The temperature ranges from 28.0 °C–30.0 °C, which showed no preference to three enterotype (Fig. 4c). However, the water salinity of this samples ranges from 34‰ to 35‰, which is at a relatively high level. As we found in this study, salinity, after temperature, was the second shaping factor to shrimp intestinal microbiota (Partial Mantel test, temperature: rM = 0.444; salinity: rM = 0.141, Table 3). Previous studies showed that the abundance of Vibrio significantly increases under high salinity conditions43, which might explain why these four samples were identified as ET V in this study. Notably, WFS-bias ET CB distributed in higher temperature ranges than ET V and ET S (Fig. 4), reflecting the observed phenomenon that WFS usually breakouts in high temperature34. More broadly, warming altered shrimp diversity and structure, as well as the transition among different enterotypes (Fig. 5). Many studies report that once the intestinal system settles into an ecological state, interactions of microbe-microbe and microbe-host interactions can prevent transition into alternate state unless they first traverse other transitional states22. This transition barrier was mainly attributed to long-term dietary and host immunity in mammal’s intestine ecosystem27. However, this strong deterministic force is seldom found in aquaculture animals. For instance, a positive correlation between temperature and alpha diversity was observed in this study (Fig. 4d). Indeed, warming regulates the expressions of immune genes including TLR, IMD and Casp3 which were especially related to the resistance to Vibrio, leading to a dramatic alteration in shrimp intestinal diversity41. In higher temperature range, the increase of Candidatus Bacilloplasma and decrease of Vibrio might lead the rising of alpha diversity, which was also consistent with the result in Fig. 3f. Commonly, aquaculture animals have much more dynamic gut microbiomes than terrestrial vertebrates due to the lower host selectivity26. Our results might potentially explain the epidemiological features for WFS, which implied that warming would mediate shrimp health by altering its enterotypes. In consideration of the phenomenon that WFS has spread from low-latitude to high-latitude areas along with temperature rising in recent years34,39, immediate actions are urgently needed to prevent farming losses when facing climate changes.

Table 3 Mantel test and partial Mantel test for environmental parameters to shrimp intestinal microbiota

Collectively, main findings of this study were synthesized into a conceptual model as depicted in Fig. 6. Shrimp IM was stratified into three enterotypes that were linked to shrimp diseases outbreaks at a global scale, and further experimentally demonstrate warming as the most sharping factor to shrimp enterotype distribution and antimicrobial peptides expression. All these results could be a newly potential ecological indicator for global warming in ectothermic animal field, and called for heightened vigilance to temperature changes in farming practices.

Fig. 6: A graphic summary for the main findings in the present study.
figure 6

Integrated with multi-omics analysis, global shrimp intestinal microbiota was stratified into three enterotypes and temperature was the most shaping factor to shrimp microbial structure. Firstly, as ambient warming, the relative expression of different antimicrobial peptides (Alfs and Pens) significantly altered. Subsequently, the relative abundance of Candidatus Bacilloplasma increased at a relative higher temperature, yet Vibrio showed an opposite trend. Differences of the dominated keystone genera led to distinctive enterotype statuses, and a briefly scheme was purposed to distinguish shrimp enterotype. For ET V, the probability of HPNS increased due to high abundance of opportunistic pathogen Vibrio and its competitive relation to other microbes. For ET CB, the probability of WFS increased for the fragile network and IM dysbiosis, while ET S showed a more stable and complex network, which thereby led to an obvious disease bias in shrimp. Taken together, warming-driving enterotypes mediated shrimp health.

Methods

Global IM data collection for pacific white shrimp

A total of 15 cohorts from nine countries, comprising 1369 IM data of pacific white shrimp were included in this study, with 1150 sites collecting IM samples by ourselves and 219 IM samples data downloaded from National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) (Brazil (4, PRJNA473127), Ecuador (55, PRJNA800805), Indonesia (40, PRJEB37200), Mexico (84, PRJNA629597, PRJNA575880 and PRJNA792615) and USA (36, PRJNA551222)). The other 699 shrimp intestinal microbiota data collected in China was previously uploaded in NCBI SRA under PRJNA391456, PRJNA542015, PRJNA486615, PRJNA89338435,49,50,51. Another 451 samples were collected in this study in Malaysia (97), Thailand (264), and Vietnam (90), and the sequencing data was uploaded in SRA under PRJNA625419 and PRJNA893384. The shrimp intestine collection was conducted with the following steps: shrimp’s surface was sterilized with 70% ethanol, and the intestine was aseptically dissected, which was then put into a 2 mL centrifuge tube containing 1 mL 75% ethanol with 2% 0.5 mol·L−1 EDTA solution. Afterwards, samples were immediately stored at −20 °C ice before DNA extraction.

The other 219 shrimp IM data from Brazil (4), Ecuador (55 data), Indonesia (40), Mexico (84) and USA (36) were downloaded from SRA under PRJNA473127, PRJNA800805, PRJEB37200, PRJNA629597, PRJNA575880, PRJNA792615 and PRJNA551222. Water physicochemical parameters of shrimp culture environment, including temperature, salinity, pH and dissolved oxygen were determined by an YSI pro plus detector (YSI, Ohio, USA).

High throughput sequencing of 16S rRNA gene and bioinformatic analysis

For the shrimp intestine sample we collected, the total DNA of shrimp intestine was extracted by QIAamp PowerFecal DNA Kit (Qiagen, Dusseldorf, Germany) following the manufacturer’s directions. The concentration and purity of total DNA were determined by NanoVuePlus Spectrophotometer (GE Healthcare, Massachusetts, USA). Primer pair 515 F (5ʹ-GTGCCAGCMGCCGCGGTAA-3’) and 806 R (5’-GGACTACHVGGGTWTCTAAT-3’) were used to amplify the V4 region of 16S rRNA gene, which was modified with a barcode tag containing a random 6-base oligos61. Sequencing libraries were generated using TruSeq DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, USA) and the library quantity was determined by Qubit 2.0 Fluorometer (Thermo Scientific, Massachusetts, USA). The 16S rDNA amplicon libraries were sequenced by Illumina Hiseq2500 and HiSeq4000 platforms. The raw data were paired-end reads and their average length was 250 bp.

Raw data quality was filtered, and the ASV profile information was obtained using dada2 method in QIIME2 software (Version 2022.9)62,63. Taxonomic assignment was conducted using a SILVA 132 reference data. The relative abundance of each ASV was normalized using the standard of sequence number corresponding to the sample with the least number of sequences. Since all the 1369 16S amplicon sequencing data contained V3, V4, and V3 + V4 regions, the rarefaction curve was conducted based on 1150 V4 region data. The following alpha diversity and beta diversity index were calculated on genus level with all 1369 shrimp intestinal microbiota data.

Alpha diversity, which shows the complexity of species in just one sample through the Shannon index and Chao1 index, was calculated using vegan package in R software (Version 4.10). The beta diversity, based on the Bray-Curtis and Canberra distance was calculated using vegan package in R, which was further visualized to evaluate the species complexity differences of samples by principal coordinate analysis (PCoA). Random forests regression was used to identify the biomarkers of healthy and diseased shrimp and predict the probability of disease using the following parameters with randomForest package in R (cv. fold = 10, step = 0.9, replication = 50, ntree = 4000).

The POD value was defined as the ratio between the number of decision trees that was voted as ‘diseased’ and the number of total sampling trees (nvotes/ntrees). Molecular ecological network was conducted using online pipeline following instruction (http://ieg4.rccc.ou.edu/mena/login.cgi)64 and was displayed using Gephi software (Version 0.9.6). The mean nearest taxon distance (MNTD) was calculated to determine which processes govern the assembly of IM65. The negative nearest taxon index (NTI) assessing ecological processes was calculated with Picante package60,66. The mean distance between each taxon and its nearest neighbor (β-MNTD) was computed reflecting the dissimilarity between communities65. Difference between the observed β-MNTD and the mean of null distribution is referred as β-NTI. The β-NTI in combination with Bray-Curtis-based Raup-Crick (RCBray) was used to determine the relative contributions of ecological processes that determine the assembly of IM. The relative influence of community turnover determined by homogeneous and heterogeneous selections was denoted by β-NTI < −2 and β-NTI > +2 fractions, respectively67,68. If |β-NTI | < 2 but RCBray > +0.95 or < −0.95, community turnover was governed by dispersal limitation or homogenizing dispersal processes; the fraction of pairwise comparisons with |β-NTI | < 2 and |RCBray | < 0.95 represented the component of compositional turnover was governed by undominated69.

Shrimp enterotype clustering based on genus profile

Genus-level enterotypes analysis was performed according to the Dirichlet multinomial mixtures (DMM) clustering protocols with Jaccard distance by “microbiome” and “DirichletMultinomial” packages in R. The optimal clustering number among the samples was calculated using Laplace approximation to the model evidence that was set from 1 to 5. The best number of clusters was assessed by the fitting effects of DMM model. To explore the component compositional difference, the Dirichlet parameter vector was used by fitting a single mixture to the data set, which arranged ASVs by assignment strength. The most important drivers at genus level were sorted by their contribution to different clustering.

Metagenomic sequencing and gene catalog construction

A subset of 131 shrimp samples were selected for metagenomic sequencing. Sequencing libraries were generated using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, Massachusetts, USA). The PE150 libraries were generated for metagenomic sequencing, which were analyzed for size distribution by Agilent2100 Bioanalyzer (Agilent, California, USA), and then were sequenced by an Illumina HiSeq2500 platform.

Quality control was conducted by Trimmomatic. The reads aligned to the shrimp genome sequence were removed with MEGAHIT (Version 1.05)70. The remaining high-quality reads were used for further analysis. The assembly of reads was executed using MEGAHIT de novo. The clean data were mapped against scaffolds using MEGAHIT. Unused reads from each sample were assembled using the same parameters. Genes (minimum length of 100 nucleotides) were predicted on scaftigs longer than 500 bp using Prodigal (Version 2.6.3)71. Afterward, a non-redundant gene catalog was constructed with Linclust (Version 2.0) using a sequence identity cut-off of 0.9. To determine the abundance of genes, reads were realigned to the gene catalog with BBMap (Version 37.68). The abundance of genes was calculated by counting the number of reads and normalizing by gene length. All genes were aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using DIAMOND (Version 0.7.9.58) with the E-value of 1E-472. Genes were assigned to the KEGG orthology by the highest scoring annotated hits. The abundance of KEGG pathway was calculated by summing the abundance of genes annotated to the same feature.

Warming experiment to shrimp IM alteration

To determine the effect of temperature in regulating shrimp IM, shrimp were cultured under different temperature. Shrimp with average weight of 7.0 ± 0.31 g were reared in the 300 L tanks at a density of 60 shrimp per tank for a three-day adaption. Then, shrimp were divided into 5 groups (n = 60 per group) and cultured under different ambient temperature for one week: (I) 20 °C; (II) 24 °C; (III) 28 °C; (IV) 32 °C; (V) 36 °C. Shrimp was fed with commercial feed (Xiajianle #2, bought from Guangdong HAID Group Co., Ltd.). This commercial contains 39% protein, 6% crude fiber, 3.5% fat, 1.8% lysine and 1% phosphorus. No antibiotic or probiotic was supplemented in the experiment. At the end of the experiment, 18 shrimp intestines were aseptically dissected from each group for 16S rRNA amplicon sequencing analysis, and 10 of these 18 intestines was also further extracted the total RNA for transcriptome analysis.

RNA interference for antimicrobial peptides in shrimp

To deeper understand the relationship between antimicrobial peptides and microbial composition variation, double stranded RNAS (dsRNAs) that targeting green fluorescent protein (dsGFP, as control), Pen2, Pen3, Pen4, Alf1, Alf2, Alf3, Alf4 were prepared by in vitro transcription using a T7 RiboMAX Express RNAi System Kit (Promega, Wisconsin, USA). The primer information was listed in Supplementary Table 6. Shrimp with average weight of 6.7 ± 0.27 g were randomly divided into eight groups (n = 20 for each group). Shrimp were injected with 50 μL PBS containing 5 μg dsRNA. At 48-h post injection, six shrimp were randomly sampled from each group to determine the RNA interference efficiency for antimicrobial peptides in intestine. Afterward, shrimp were cultured for five days. At last, for each group, six shrimp intestine was collected in to target the microbial compositional variations.

Transcriptome sequencing and expression examination for antimicrobial peptides

A total of 10 shrimp were randomly selected in each group. Each intestine was individually used for RNA extraction using the RNeasy Mini Kit (Qiagen) as one sample, and further sent for transcriptome sequencing. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB) and index codes were added to attribute sequences to each sample. The libraries were sequenced by an Illumina NovaSeq6000 platform. Clean data were obtained by removing reads containing adapter and low-quality reads from raw data using Trimmomatic. The clean reads were then mapped to the reference genome sequence and annotated based on the reference genome (NCBI Genome Database Number: ASM378908v1) using Hisat2 (Version 2.0.4)73. Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped using StringTie software (Version 1.3.4 d)74. Differential expression analysis of two groups was performed using DESeq2 software (Version 3.15)75.

To assess the knockdown efficiency of antimicrobial peptides, real-time PCR was performed with SYBR Premix Ex Taq™ II (Takara, Shiga, Japan) on a LightCycle 480 System (Roche, Basel, Swiss). Expression was determined using the 2−ΔΔCt method after normalization to the control group that using dsGFP as knockdown control.

Analysis for relationship between IM and environmental factors

Mantel test was used to test correlations between two distance matrixes, and partial Mantel was used to determine the contributions of various factors to explain community variations76. The significance was measured by random permutations (n = 999). For environmental variables, the best combination of environmental variables was selected with ranked correlations by “BioENV” function in “vegan” package. The selected environmental variables were normalized. Partial Mantel test was used to determine the contributions of various factors to explain community variations using vegan package77. Aggregated boosted tree was conducted to evaluate the relative influence of environmental variables on microbial community diversity using vegan package78. Structural equation modeling was performed to further test direct and indirect effects of geographic and environmental factors on microbial diversity by AMOS software (Version 24.0.0). A priori model was initially established according to empirical data and literature, which was then adjusted by removing non-significant links or variables and adding new links until the lowest χ2 value was archived with P > 0.05, which means the model fits the observed data.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.