Microbial signature profiles of Penaeus vannamei larvae in low-survival hatchery tanks affected by vibriosis

Ethics statement

No ethics statement or permissions were required for this study. However, in order to ensure the ethical treatment of animals, a two-step anesthesia-euthanasia protocol of animal submersion in 60% and 95% alcohol was used to sacrifice them.

Sources of samples and microbiome data

An observational study was conducted in a commercial shrimp hatchery (64 tanks, 35 tons, South America) to investigate the microbiome of P. vannamei larvae naturally affected with vibriosis (Fig. S1). Two investigations were generated from the samples collected in the observational study (Fig. S1). In the first investigation, already published by the authors of this article (Reyes et al., 2022), the microbiome of larvae collected from high (22 samples from seven tanks) and low (25 samples from other seven tanks randomly selected from 51 tanks) survival tanks affected by AHPND was compared for future development of probiotics to mitigate AHPND (Fig. S1). In the second research (this article), we re-used the microbiome data obtained from the low survival AHPND-affected tanks of our first study (accession number PRJNA800805). This data was then compared with the microbiome data from zoea 2 syndrome-affected larvae (20 samples from five different tanks, microbiome data presented for the first time in this article) with the new aim of comparing the microbial profiles of larvae affected by two types of vibriosis: AHPND and zoea 2 syndrome (Fig. S1).

Sample collection

Samples of larvae of the following stages were collected from each one of the 12 tanks: Mysis 3 (M3), Postlarvae 4 (PL4), Postlarvae 7 (PL7), and Postlarvae 10 (PL10), except in two of the seven tanks affected by AHPND, where only samples for M3, PL4, and PL7 were collected in each of the two tanks because the populations died at the PL9 stage (Fig. S1). In addition, a sample (M3) from a third tank affected by AHPND failed DNA quality control (Fig. S1). In conclusion, a total of 45 samples of P. vannamei larvae from the 12 commercial hatchery tanks (25 and 20 samples from tanks affected by AHPND and zoea 2 syndrome) were used for the microbiome analysis of this article (Fig. S1). Survival [mean ± standard error of the mean (SEM)] at harvest from the 12 tanks was 24.4 ± 5.2%. In addition, observations for clinical signs of disease were made at each sampling (5 to 7 min) at the larval stages PL7 and PL10. Ranges of water quality parameter were pH = 7.9–8.1, salinity = 35 g/L, dissolved oxygen concentration = 5.0 mg/L, and temperature = 32.5–33 °C. Larvae were supplemented with four commercial diets (Skretting PL, Frippak Fresh, Zeigler Z Plus, and Advance Brine Shrimp).

Sample processing and high-throughput sequencing

The samples were processed for three kinds of analysis: (1) Polymerase chain reaction (PCR) analysis to detect bacterial infections caused by AHPND (detection of PirAB toxin genes) and discard infections with common shrimp pathogens (white spot syndrome virus–WSSV, infectious hypodermal and hematopoietic necrosis virus–IHHNV, Enterocytozoon hepatopenaei–EHP, and shrimp hemocyte iridescent virus–SHIV), (2) histological analysis to confirm the presence of lesions compatible with bacterial infections, and discard lesions caused by other common shrimp pathogens, and (3) microbiome analysis through 16S rRNA amplicon high-throughput sequencing (HTS).

Approximately, 1 g of each sample (M3, PL4, PL7, and PL10) was washed with a 2% NaCl sterile solution to remove external impurities, macerated, and divided into two aliquots. The first aliquots of macerates (500 µl) for each sample were processed for PCR analysis. The genomic DNA (gDNA) extraction was carried out using the phenol-chloroform protocol. The PirAB toxin genes were amplified with the AP4 F1/R1 - AP4 F2/R2 primers (Dangtip et al., 2015). WSSV, IHHNV, EHP, and SHIV infections were diagnosed with the IK1/2-IK3/4, IHHNV309F/R, SWP 1F/1R-SWP 2F/2R, and SHIV-F1/R1-SHIV-F2/R2 primers, respectively (Lo et al., 1996; Tang, Navarro & Lightner, 2007; Jaroenlak et al., 2016; Qiu et al., 2017).

Twenty larvae of each sample collected at the PL7 and PL10 stages were preserved in Davidson’s fixative solution for histological analysis. Shrimp tissues were processed according to a previously published protocol where the sections were cut at 5 um and stained with Mayer-Bennet haematoxylin and eosin (Bell, 1988).

The 45 samples were allocated into two groups according to the disease condition: affected by AHPND: samples collected from tanks where AHPND was detected by PCR (at any larval stage) and histological analysis (only for PL7 and PL10 larval stage), and affected by zoea 2 syndrome: samples from tanks where AHPND was negative by PCR (at any larval stage) and also by histological analysis (only for PL7 and PL10 larval stage), but histological lesions were compatible with zoea 2 syndrome.

The second aliquots of macerates (500 µl) obtained from each sample were distributed into 2-mL cryovial tubes for a snapshot freezing with liquid nitrogen (−196 °C) and immediately stored at −80 °C for further microbiome analysis and its functional prediction. These macerates remained at −80 °C until all samples were collected at the end of the production cycle of all 12 tanks. Then, all macerates were processed for gDNA extraction using the ZymoBiomics DNA Microprep Isolation Kit (Zymo Research, Irvine, CA, USA), following the manufacturer’s instructions. A 20 µl aliquot of gDNA was diluted in DNase-free ultrapure water, and DNA quality and concentration (A260/280 = 1.8−2.0) were examined through a NanoDrop One Microvolume UV-Vis Scanning Spectral Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The gDNA of all 45 samples was then submitted to Novogene Incorporation (Beijing, China) for amplification of the V3–V4 region of the 16S rRNA gene using the 341F/806R primers (Takahashi et al., 2014), amplicon library construction, and paired-end sequencing with the Illumina NovaSeq 6000 P250 platform. The DNA sequences of the 20 samples from the five tanks affected by zoea 2 syndrome were deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession number PRJNA870476. The sequences of the 25 samples from the seven tanks affected by AHPND are available under accession number PRJNA800805 (Reyes et al., 2022).

Sequence analysis

Sequencing reads were processed with the Divisive Amplicon Denoising Algorithm (DADA2) version 1.16 (Callahan et al., 2016). Sequences were filtered with the filterAndTrim function to obtain reads higher than Q30. Sequences with more than two expected errors were removed [maxEE = c (2,2)], and the last 10 nt of the forward reads and the last 20 nt of the reverse reads were trimmed [truncLen = c (240,230)]. The filtered reads were dereplicated, denoised, and merged with the forward and reverse amplicon sequences. Chimeric sequences were then removed using the removeBimeraDenovo function. The table of amplicon sequence variants (ASVs) was generated, and the taxonomy was assigned using the public database SILVA version 138.1 (McLaren Michael & Callahan Benjamin, 2021).

Data analysis

Larvae survival at harvest from tanks affected by AHPND and zoea 2 syndrome were compared by the non-parametric Wilcoxon-Mann–Whitney test using the stats package in R.

The variability of the alpha diversity indices and the relative abundance of the microbiome were studied for the samples collected at the four larval stages of tanks affected by AHPND and zoea 2 syndrome. Shannon and ACE indices were used to estimate the alpha diversity (Shannon & Weaver, 1949; Hughes et al., 2001). Differences in Shannon and ACE indices, and relative bacterial abundance at family and genus levels between the disease conditions were evaluated through the non-parametric Wilcoxon-Mann–Whitney test.

The microbiome structure, and the Vibrio communities among larval samples (beta diversity) collected from tanks affected by the two diseases condition was estimated by the multivariate analysis of the Bray–Curtis dissimilarity index and visualized by principal coordinate analysis (PCoA). The null hypothesis of equality of microbiomes and Vibrio communities between samples collected from tanks affected by AHPND and zoea 2 syndrome was evaluated by a similarity analysis (ANOSIM). All previous statistical analyses were considered significant with P < 0.05.

Co-occurrence network analyses at the phylum and genus levels were performed to evaluate the interactions among bacterial taxa at each disease condition. The co-occurrence network was built through a Spearman correlation matrix of the ASV abundances using the trans_network$new function. Co-occurrence was considered robust when positive or negative correlation coefficients were greater than 0.6 and P < 0.01 (Barberán et al., 2012). The cal_sum_links function was used to calculate for each ASV (edge) the number of other ASVs with significant correlations (P < 0.01) to the correspondent edge and plotted with the plot_sum_links function.

Differential abundance analysis between the disease conditions was evaluated through a linear discriminant analysis (LDA) of the ASVs with effect size (LEfSe) (Segata et al., 2011). ASVs with an LDA cutoff of 2.0 and significant differences (P < 0.05) between both disease conditions were considered potential microbial signatures. The robustness of the microbial signatures obtained with the LEfSe analysis was evaluated through a Bayesian t-test analysis using the BayesFactor package (Morey & Rouder, 2011). The Bayes factor BF₁₀ was calculated as the ratio of the probability of H₁ (different abundance of ASVs between both disease conditions) over the probability of H₀ (equal abundance of ASVs between both disease conditions) based on the data. BF₁₀ with values: 1–3, 3–10, and >10 showed weak, moderate, and strong evidence in favor of H₁ (Morey & Rouder, 2011; van Doorn et al., 2021).

The ability of the bacterial functionality was predicted using the 16S rRNA amplicon data against the pathway reference profiles (Ref99NR) of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2012). Tax4fun2 was used to predict the functionality of the ASVs at hierarchical levels 1, 2, and 3 (Wemheuer et al., 2020). ASVs with similarities higher than 97% to the reference database were considered Orthologous KEGG groups (KOs). A differential abundance analysis between both disease conditions was performed to identify significant (P < 0.05) functional biomarkers through LEfSe analysis of the KOs, with an LDA cutoff of 2.0. All downstream analyses, except the Bayesian analysis, were performed with the microeco package (Liu et al., 2021) in the R software version 4.2.1 (R Core Team, 2022).

Sequence analysis

The sequence analysis of the 45 samples generated approximately six million clean reads, with an average of 146,311 clean reads per sample (Table S1). The error rate across the 45 samples averaged 0.03% for all the clean reads (Table S1). On average, 88.06% of these reads exhibited a Phred score Q >30 and a GC content of 53.72% (Table S1). Amplicons were clustered into 8,617 ASVs, with an average of 392 ASVs per sample (Table S1). Good’s coverage was higher than 99.99%, indicating an optimal sequencing depth (Table S1).

Microbiome composition and statistical analysis

No significant differences in the Shannon and ACE indices of the larvae microbiome were found between disease conditions (P = 0.876 and 0.100, Fig. S2). The distribution of the relative abundance of the microbiome showed the prevalence of 17 bacterial phyla, representing 29 classes, 75 orders, 105 families, and 182 genera in all larval samples.

At the phylum level, Pseudomonadota exhibited the highest and most similar abundance in tanks affected by AHPND (72%) and zoea 2 syndrome (77%) (Table S2 and Fig. S3). Bacteroidota (16%) and Bacillota (3%) were the following phyla with the highest global average abundances (Table S2 and Fig. S3).

At the family level, Rhodobacteraceae, Vibrionaceae, and Flavobacteraceae exhibited the highest abundance at both AHPND and zoea 2 syndrome-affected tanks (Fig. 2 and Table S2). Vibrionaceae exhibited almost a two-fold increase of abundance ( P = 0.015) in tanks affected by zoea 2 syndrome, especially at the M3 stage (49%), in contrast to the same larvae stage in AHPND-affected tanks (11%) (Table S2). Flavobacteraceae was most abundant at the M3 stage from AHPND-affected tanks, being three-fold more abundant than at the same larval stage from tanks affected by zoea 2 syndrome, but overall was not significantly different between the two diseases (P = 0.120) (Fig. 2 and Table S2).

At the genus level, Marinibacterium exhibited the highest abundance, especially at the PL4 and PL7 (22%) stages of tanks affected by zoea 2 syndrome (Fig. 2B and Table S2). Catenococcus and Vibrio were the second and third most abundant genera, respectively, in tanks affected by zoea 2 syndrome, with almost a two-fold increase compared to the corresponding AHPND cohort (P = 0.016, and P = 0.013, respectively) (Fig. 2B and Table S2). Gilvibacter was seven-fold more abundant (P = 0.044) in AHPND-affected tanks, especially at the M3 stage (18%). Bacillus was the least abundant genus in larvae at both disease conditions (Table S2).

A notable difference was observed when bacterial communities (beta diversity) were compared between the two disease conditions. The ANOSIM analysis of the beta diversity showed a remarkable separation of the larval microbiome affected with both disease conditions (P = 0.005, Fig. 3A). Moreover, significant differences in the shrimp microbiome were observed between disease conditions at M3 (P = 0.005), PL7 (P = 0.039), and PL10 (P = 0.040) stages, except at PL4 (P = 0.130). Closer analysis of the Vibrio genus revealed a significant separation (P = 0.001) between the two diseases (Fig. 3B).

The co-occurrence bacterial network at the phylum level for the samples collected from tanks affected by AHPND showed a higher complexity and connectivity than in samples collected from tanks affected by zoea 2 syndrome (Fig. 4, Table S3). This higher complexity of the microbial communities of samples collected from tanks affected by AHPND was reflected by higher values (between 20% and 2.2-fold higher) of the number of edges, network centralization, density, clustering coefficient, heterogeneity, and average degree, and a lower average path length (Table S3). Such complexity was observed despite the bacterial network for both disease conditions showed the same network diameter, and only a 0.3% higher number of nodes for the bacterial network of samples collected from tanks affected by zoea 2 syndrome (Fig. 4, Table S3).

Pseudomonadota was directly associated with 251 edges at the phylum level in samples affected by AHPND (Fig. 4). A total of 112 of these edges were connected to ASVs of the same Pseudomonadota phylum, 105 edges were connected to Bacteroidota ASVs, and the remaining 34 edges were connected to ASVs of the following phylum: Desulfobacterota, Bacillota, Myxococcota, Campylobacterota, Actinomycetota, Bdellovibrionota, Chloroflexi and Cyanobacteria (Fig. 4). In samples affected by zoea 2 syndrome, Pseudomonadota was directly associated with 220 edges, of which 146 were connected to ASVs of the same Pseudomonadota phylum, 41 edges were connected with Bacteroidota ASVs, and the remaining 33 edges were connected with ASVs of the following phylum: Actinomycetota, Desulfobacterota, Bdellovibrionota, Bacillota, Myxococcota, Chloroflexi and Patescibacteria (Fig. 4).

The network analysis at the genus level showed a high number of ASVs without taxonomic identification at both disease conditions (Figs. S4, S5). In the co-occurrence network of the AHPND-affected larvae, Vibrio edges were mostly associated with Catenococcus (Fig. S4), whereas 23 of the 38 edges in the co-occurrence network of the samples affected by zoea 2 syndrome were associated with ASVs belonging to the same Vibrio genus, being this number of edges 2.5-fold higher in larvae affected by zoea 2 syndrome compared to larvae affected by AHPND (Figs. S4, S5).

The LEfSe analysis at the genus level identified 22 significant ASVs, ten of which were significantly more abundant for AHPND-affected larvae (Table S4 and Fig. 5). The genera Sneathiella, Cyclobacterium, Haliea, Catenococcus, and Lewinella were abundant in AHPND-affected larvae, with Sneathiella being the most abundant (Fig. 5). On the other hand, the genera Catenococcus, Spongiimonas, Vibrio, Meridianimaribacter, Tenacibaculum, and Thalassotalea were abundant in the tanks affected by zoea 2 syndrome (Table S4 and Fig. 5). The Bayesian analysis confirmed the evidence (strong and moderate) of the differential abundance of ASVs found through the LEfSe analysis (Table S5). However, two ASVs from tanks affected by zoea 2 syndrome (ASV 56 and ASV 114), corresponding to the Thalassotalea and Vibrio genera showed weak evidence of differential abundance between the two disease conditions (Tables S4, S5).

A total of 22 bacterial pathways were differentially represented in the comparison analysis between disease conditions using the LEfSe analysis (P < 0.05) (Fig. 6). At a higher level (KEGG level 3) the function pathway (metabolism-related) was significantly more abundant in samples affected by AHPND, while cellular processes and environmental information processing were significantly more abundant in samples affected by zoea 2 syndrome (Fig. 6). At lower depth levels (KEGG levels 1 and 2) the following functional pathways were significantly more abundant in AHPND-affected tanks: amino acid and carbon metabolisms, metabolic pathways, biofilm formation, translation, two-component system, signal transduction, global and overview maps, biosynthesis of antibiotics, amino acids, and secondary metabolites (Fig. 6). In contrast, cell motility, phosphotransferase system, quorum sensing, cellular community—Prokaryotes, ABC transporter, and membrane transport were significantly more abundant in samples affected by zoea 2 syndrome (Fig. 6).