Composition, diversity and function of intestinal microbiota in pacific white shrimp (Litopenaeus vannamei) at different culture stages

Sample collection

From July to October 2015, 22 intestine samples were collected from 5 shrimp ponds in a commercial shrimp farm, Maoming, Guangdong, China (21.68°N, 110.88°E). Healthy shrimp were collected every 15 days from the larval stage (stage1, 15 dph) to adult stage (stage 5, 75 dph) (Table S1).

Each pond was approximately 2,600 m² and the average depth was 1.5 m. Shrimp larvae with average length of 0.7 cm were cultured at a stocking density of 200,000 shrimps each pond. The water temperature was relatively stable at approximately 32 °C. The pH value ranged from 7.5 to 8.61. The concentration of NH₃-N, NO₂-N, NO₃-N, PO4 3 − and SO4 2 − were in range of 0.0089∼1.1095 mg L⁻¹, 0.0022∼0.9869 mg L⁻¹, 0.0323∼3.3007 mg L⁻¹, 0.0171∼0.3131 mg L⁻¹ and 0.0012∼0.3777 mg L⁻¹. There was no antibiotic application during the culture period. Some probiotics, including Lactobacillus and Bacillus from Guangdong Zhongtai Biology Co., Ltd. (Guangdong, China), have been mixed with feed and applied to ponds once a week.

Sampling was according to the previously reported methods (Oxley et al., 2002; Rungrassamee et al., 2014). The shrimp’s surface was sterilized with 70% ethanol and the intestine was aseptically dissected. The intestine was put into a 2 mL centrifuge tube which contained sterile glass beads and 1.5 mL PBS buffer. The tube was thoroughly vortexed for 3 min and centrifuged at 10,000 g for 1 min. Samples were immediately stored at −80 °C before DNA extraction.

DNA extraction and sequencing

Total DNA was extracted by the PowerFecal DNA Isolation Kit (MoBio, Palo Alto, CA, USA) following the manufacturer’s directions. The concentration and purity of total DNA were determined by NanoVuePlus Spectrophotometer (GE Healthcare, USA) and 1% agarose gels. The primer pair 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used to amplify the V4 hypervariable region of 16S rRNA gene, which was modified with a barcode tag with a random 6-base oligos (Bates et al., 2011). Sequencing libraries were generated via using TruSeq DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, USA). In addition, the library quantity was assessed on Qubit 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA). The libraries were sent for sequencing by Illumina Hiseq2500 platform (Illumina, San Diego, CA, USA), which was conducted by Novogene Bioinformatics Technology Co.,Ltd. (Beijing, China). Raw data generated from Hiseq2500 platform were paired-end reads.

Data analysis

Based on the unique barcode, sequences were assigned to samples and then removed off the barcode and primer sequence by QIIME (Version 1.7.0, http://qiime.org/index.html) (Caporaso et al., 2010). In order to merge paired-end reads when at least some of the reads overlap the read generated from the opposite end of the same DNA fragment, FLASH (Version 1.2.7, http://ccb.jhu.edu/software/FLASH/) was used to get raw tags (Magoc & Salzberg, 2011). In terms of quality control, raw tags with low quality (quality value  ≤19, homopolymers  ≥3 bases and sequence length ≤200 bp) were filtered by QIIME according to the QIIME quality filtering process in a bid to obtain the high-quality clean tags. Tags were compared with Gold database (http://drive5.com/uchime/uchime_download.html) by UCHIME algorithm (http://www.drive5.com/usearch/manual/uchime_algo.html) so as to remove off chimera sequences and then the effective tags were finally gained (Edgar et al., 2011).

Sequences with over 97% similarity were considered as the same OTUs for further annotation (Edgar, 2013). To align the sequences, the GreenGene Database (http://greengenes.lbl.gov/download) was used as a reference database (DeSantis et al., 2006). Later, the taxonomic information was annotated by RDP classifier (Version 2.2, http://sourceforge.net/projects/rdp-classifier/) with 80% confidence threshold. OTUs abundance information was normalized using a standard of sequence number corresponding to the sample with the least sequences. The Venn diagram, which was used to find out the shared OTUs among groups, was conducted by Draw Veen Diagram online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). Alpha diversity, showing the complexity of species for one sample through 5 indices, including Chao, Shannon, Simpson, ACE and Good’s coverage, was calculated by QIIME following the tutorial (http://qiime.org/scripts/alpha_diversity.html) and displayed via R software (Version 2.15.3). Beta diversity, used to evaluate differences of samples in species complexity, was calculated by QIIME (http://qiime.org/scripts/beta_diversity.html). Unweighted pair-group method with arithmetic means (UPGMA) was conducted to report the hierarchical clustering of samples by QIIME following the guidance (http://qiime.org/scripts/jackknifed_beta_diversity.html). Statistical analyses of alpha diversity were calculated by analysis of variance (ANOVA) to compare the significant differences at different culture stages by SPSS (Version 21). Multiple-response permutation procedure (MRPP) was conducted to test significant difference between any two of compared culture stages using the vegan package in R (Cai, 2006). Permutational analysis of multivariate dispersions (PERMDISP) was used to test whether the microbial community varied at different culture stages by using the vegan package in R (Anderson, 2006). Permutational multivariate analysis of variance (PerMANOVA) was conduct to compare microbial composition and function dissimilarities (Anderson, 2001). A calculated P value < 0.05 was considered to be statistically significant.

Microbial function prediction based on 16S rDNA data

The OTU table was used to generate the inferred metagenomic data by using PICRUSt (version 1.1.0) to predict the metagenomic functional capacity (Langille et al., 2013). The abundance values of each OTU were firstly normalized to its 16S rRNA copy number respectively. Predicted functional pathways were annotated by using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2012) at level 2 and level 3 KEGG orthology groups (KOs) (Langille et al., 2013). The accuracy of the predictions of the metagenomes was assessed by computing the nearest sequenced taxon index (NSTI). The associated metabolic pathways were identified by means of employing the HMP unified metabolic analysis network (HUMAnN) (Abubucker et al., 2012). Moreover, the KOs at different culture stages were further examined by PerMANOVA. The relationships among functional capacities were analyzed by principal component analysis (PCA).

Accession number

The raw data in this study have been deposited in the GenBank Sequence Read Archive database. The accession number is SRX2946975.

Composition of intestinal microbiota

Quality and chimera filtration of the raw data produced totally 1,408,105 high quality sequencing reads from 22 samples, belonging to five culture stages, with an average of 64,005 reads, ranging from 41,250 to 79,515 (Table 1). By performing the alignment at an average length of 253 bp, OTUs were clustered at 3% distances, among which each OTU represented a unique phylotype. Finally, 2,496 OTUs were obtained and the number of OTUs detected in each sample ranged from 585 to 1,239, with an average of 880 OTUs (Table 1).

OTUs were identified into 43 phyla. Sequences that could not be classified into any known groups were assigned as ‘others’. The most relative abundant phyla in all samples were Proteobacteria (63.5%), Cyanobacteria (7.0%), Tenericutes (6.5%), Fusobacteria (5.3%), Firmicutes (4.1%), Verrucomicrobia (3.6%), Bacteroidetes (3.6%), Planctomycetes (2.9%), Actinobacteria (0.8%) and Chloroflexi (0.4%) (Fig. 1). Proteobacteria was the most abundant phylum among 21 samples except sample E2, in which Cyanobacteria was the most abundant phyla (27.8% relatively abundance).

At genus level, a total of 289 taxa were identified. The top 10 genera were Candidatus_Xiphinematobacter (3.4%), Propionigenium (3.4%), Synechococcus (2.7%), Shewanella (1.3%), Cetobacterium (1.1%), Bacillus (0.9%), Robiginitalea (0.7%), Fusibacter (0.5%) and Arcobacter (0.5%) (Fig. S1). The abundance of Lactobacillus and Bdellovibrio were 0.04% and 0.002% respectively.

For further investigation of the dominant microbiota that exists in all samples, Veen diagram was constructed to identify dominant OTUs presented in intestine (Fig. 2). There were 59 OTUs shared among 22 samples, representing 83.1% of the total reads. Among the shared OTUs, 17 OTUs (28.8%) belonged to Proteobacteria, 5 OTUs (8.5%) belonged to Cyanobacteria and 5 OTUs (8.5%) belonged to Bacteroidetes.

Diversity, similarity and function analysis

The diversity and richness indices of all samples from five ponds were calculated in an effort to illustrate the complexity of each sample (Table S2). The completeness of sequencing was estimated with Good’s coverage, showing the probability of a randomly sequence already detected in the sample. The rarefaction curves approached the plateau (Fig. S2). The Good’s coverage ranged from 0.989 to 0.996, suggesting that additional 90 to 250 reads needed to be sequenced before discovering new OTUs. Shannon index and Simpson index were often used to quantify the diversity. The Shannon index ranged from to 1.936 to 6.592, while the Simpson index ranged from 0.273 to 0.968. The richness of each sample was calculated via Chao index and ACE index. Chao index ranged from 558 to 1,386, while ACE index ranged from 599 to 1,416.

The similarity and difference in different intestinal microbiota samples were further investigated. UPGMA clustering showed that almost all of the individual samples were clustered into groups according to the culture stage (Fig. 3). The OTU number, Shannon index, Simpson index, ACE index and Chao index were shown in boxplot graph (Fig. 4), and ANOVA showed that there was extremely significant difference in the OTU number, Shannon index, Simpson index, ACE index and Chao index at different culture stages (P value < 0.05) (Table 2). MRPP and PERMDISP showed the intestinal microbiota differed significantly between any two of compared stages (P value < 0.05) (Table 3). Among the top 10 phyla, PerMANOVA demonstrated the abundance of Proteobacteria, Fusobacteria, Tenericutes, Verrucomicrobia, Planctomycetes and Chloroflexi changed significantly at different culture stages (P value < 0.05) (Fig. 5).

Functional prediction of the intestinal microbiota

The changes in the presumptive functions of the intestinal microbiota of pacific white shrimp were examined by predicting the metagenomes using PICRUSt. The accuracy of the prediction was evaluated by computing the NSTI, and the mean of the samples was 0.171 ± 0.023. The metagenomic prediction showed the intestinal microbiota was enriched with functions that were related to transporters, ATP-binding cassette (ABC) transporters, DNA repair and recombination proteins, two component system, secretion system, bacterial motility proteins, purine metabolism, ribosome, pyrimidine metabolism, peptidases and transcription factors (Table 4). The relative abundance of transporters(5.25%) and ABC transporters (3.06%) were at the highest level during all culture stages. Two component system (2.49%) as well as DNA repair and recombination proteins (2.43%) were the thirdly and fourthly most abundant KOs. The range of KOs related to membrane transport, cell motility, energy metabolism, and the signal transduction was wide, while other KOs varied at small range. Moreover, PerMANOVA showed there were totally 199 KOs shifting significantly at different stages (P value < 0.05), including the KOs which belonged to amino acid metabolism, carbohydrate metabolism, energy metabolism, membrane transport and nucleotide metabolism (Fig. 6). PCA revealed that the functions of intestinal microbiota from the same culture stages were clustered closer, with the first two components explaining a total of 58.67% of the variation (Fig. 7). The results suggested the functional KOs of the intestinal microbiota varied a lot according to different culture stages.