Applied Microbiology and Biotechnology https://doi.org/10.1007/s00253-018-8891-y

APPLIED GENETICS AND MOLECULAR BIOTECHNOLOGY

Integrating molecular and ecological approaches to identify potential polymicrobial pathogens over a shrimp disease progression Wenfang Dai 1,2 & Weina Yu 1,2 & Lixia Xuan 1 & Zhen Tao 1 & Jinbo Xiong 1,2 Received: 2 October 2017 / Revised: 16 February 2018 / Accepted: 25 February 2018 # Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract It is now recognized that some gut diseases attribute to polymicrobial pathogens infections. Thus, traditional isolation of single pathogen from disease subjects could bias the identification of causal agents. To fill this gap, using Illumina sequencing of the bacterial 16S rRNA gene, we explored the dynamics of gut bacterial communities over a shrimp disease progression. The results showed significant differences in the gut bacterial communities between healthy and diseased shrimp. Potential pathogens were inferred by a local pathogens database, of which two OTUs (affiliated with Vibrio tubiashii and Vibrio harveyi) exhibited significantly higher abundances in diseased shrimp as compared to healthy subjects. The two OTUs cumulatively contributed 64.5% dissimilarity in the gut microbiotas between shrimp health status. Notably, the random Forest model depicted that profiles of the two OTUs contributed 78.5% predicted accuracy of shrimp health status. Removal of the two OTUs from co-occurrence networks led to network fragmentation, suggesting their gatekeeper features. For these evidences, the two OTUs were inferred as candidate pathogens. Three virulence genes (bca, tlpA, and fdeC) that were coded by the two candidate pathogens were inferred by a virulence factor database, which were enriched significantly (P < 0.05 in the three cases, as validated by qPCR) in diseased shrimp as compared to healthy ones. The two candidate pathogens were repressed by Flavobacteriaceae, Garvieae, and Photobacrerium species in healthy shrimp, while these interactions shifted into synergy in disease cohorts. Collectively, our findings offer a frame to identify potential polymicrobial pathogen infections from an ecological perspective. Keywords Shrimp gut microbiota . Health status . Random Forest model . Virulence gene . Co-occurrence network

Introduction Gut microbiota contributes fundamental roles in improving host fitness, such as regulating host important metabolic and physiological functions, and barrier against pathogen invasion (Dai et al. 2017; Kamada et al. 2013; Krajmalnikbrown et al. 2012; Xiong et al. 2017a). As a result, dysbiosis in the gut microbiota generally causes host metabolic syndrome and disease (Kasubuchi et al. 2015; Zhu et al. 2016). To suit the Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00253-018-8891-y) contains supplementary material, which is available to authorized users. * Jinbo Xiong [email protected] 1

School of Marine Sciences, Ningbo University, Ningbo 315211, China

2

Collaborative Innovation Center for Zhejiang Marine High-Efficiency and Healthy Aquaculture, Ningbo 315211, China

remedy to the case, it is of foremost importance to identify the causal agents. It has been proposed that gut bacteria missing in diseased individuals are candidate probiotics (Alou et al. 2017). Following this idea, it is likely that the enriched detrimental gut bacteria are potential pathogens. Traditionally, a single pathogen was isolated from host disease tissues (Carraro et al. 2011). However, this approach is time-consuming and fails to reproduce ecological niches and symbiotic relationships in situ, which may bias the identification of causal agents (Carraro et al. 2011; Ndoye et al. 2011). In addition, there is increasing evidence that some disease are caused by polymicrobial pathogens (Peters et al. 2012). The invasion of a pathogen is an ecological process, which is influenced by resident species, such as competition with gut commensals (Mallon et al. 2015). For these reasons, distinguishing taxa that contribute the divergences in gut microbiotas between host health status may provide a new avenue for identifying potential pathogens. Indeed, studies have repeatedly shown that host disease severity and incidence are closely associated with the degree of gut microbiota

Appl Microbiol Biotechnol

dysbiosis (Boursier et al. 2016; Carding et al. 2015). During health, a consortium of microbes provides a protective effect, while a disruption (i.e., loss of gatekeepers) in these interactions can begin a self-perpetuating cycle of dysbiosis (Byrd and Segre 2016; Widder et al. 2014). For instance, using an ecological network, multiple taxa that associated with host health status are identified, and removal of these gatekeepers causes collapse of microbial community (Wagner et al. 2017). Similarly, Buffie et al. (2015) found that Clostridium scindens conferred colonization resistance to Clostridium difficile infection. This probiotics effect was further validated by adoptive transfer. Thus, screening the resident taxa that antagonize pathogens facilitates the identification of probiotics. By distinguishing between age- and disease-discriminatory taxa, our recent work has screened gut bacterial biomarkers for quantitative diagnose shrimp health status with 91.5% accuracy (Xiong et al. 2017b). In addition, a bacterial pathogen 16S rRNA gene database based on the taxonomic list derived from the virulence factors database using 99% cutoff has been established, which allows researchers to screen potential pathogens (Chen et al. 2016). This method has been widely applied to identify potential pathogens in soil (Sun et al. 2017) and domestic poultry feces (Zhuang et al. 2017). Collectively, these available findings provide clues to identify candidate pathogens and resident taxa that are a barrier against pathogens, which in turn facilitates rational design of microbiota-based diagnostics and therapeutics from an ecological perspective. To date, the identification of candidate pathogens primarily applies bacterial 16S rRNA gene as biomarker (Zhu et al. 2016). However, some pathogenic bacteria, i.e., Vibrio harveyi and Vibrio campbellii, share identical 16S sequences (Gauger and Gómezchiarri 2002; Sawabe et al. 2007). Thus, it is inadequate to accurately infer pathogens via solely relying on the bacterial 16S rRNA gene. Currently, virulence genes have been served as additional molecular markers for validating potential pathogens. For example, vhh and vch genes have been integrated to distinguish pathogenic bacteria V. harveyi (Conejero and Hedreyda 2004) and V. campbellii (Luis and Hedreyda 2006). In this regard, detecting the virulence genes would assist the identification of potential pathogens. In this study, we tracked the dynamics of shrimp gut bacterial assembly over a disease progression. Shrimp samples were collected from identically managed ponds, thereby ruling out the effects of host genetics, diet, and rearing conditions that are known determinable factors in shaping the variations in gut microbiotas (Benson et al. 2010; Ussar et al. 2015). This design allowed us: (i) to characterize the divergence of gut bacterial communities between healthy and diseased shrimp, and (ii) to identify candidate pathogens. The characteristics of candidate pathogens were integrally validated by similarity percentage (SIMPER) analysis, a random Forest model, gatekeeper feature, and abundances of virulence genes. Given the fundamental roles of resident taxa in barrier

against pathogens, we further explored how the resident taxa interact with the candidate pathogens to sustain shrimp health.

Materials and methods Experimental design and sample collection The investigated shrimp ponds are located in Ningbo, Eastern China. The greenhouse ponds are cement bottomed with uniform size (2000 m2) and water depth (1.2 m), which were managed identically in terms of disinfectant seawater (salinity 25‰) inputs, 5% daily water exchange, water temperature (28 °C), stocking density (360,000 individuals per pond), feed type, and schedule. Congeneric larval shrimp (Litopenaeus vannamei) were inoculated into each pond on 8 April 2016. A disease emerged at the adult stage of shrimp in the three of our monitored six ponds on 1st July (84 days after inoculation). Given the disease exacerbation, shrimp were urgently harvested on 10th July. We collected healthy and diseased samples from each pond on 1st, 4th, and 10th July, respectively. The diseased shrimp exhibited typical symptoms of white fecal syndrome, such as inactivity, lack of appetite, red hepatopancreas, white guts, and white fecal strings (Sriurairatana et al. 2014), which is distinct of these infected by white spot syndrome virus (WSSV, the most common viral pathogen in shrimp, disease signs were loose cuticle, reddish discoloration, and white spots on the inside surface of the carapace) (Reddy et al. 2013). On 10th July, WSSV loads were comparable (unpaired t test, P = 0.452) between shrimp health status, with 5.64 ± 0.86 × 103 copies per milligram tissue in healthy and 5.96 ± 1.23 × 103 in diseased shrimp. These low loads (< 105 copies per mg tissue) are insufficient to cause L. vannamei disease (Qiao et al. 2015). Under these premises, WSSV was not the primarily causal agent here. To improve statistical power, two to three technical replicates for a given pond were employed. In total, we analyzed 56 biological samples, with 26 healthy and 30 diseased cohorts (Table S1). Shrimp were separately stored and aerated in tanks with water from the corresponding ponds during transportation.

Gut sample collection and genomic DNA (gDNA) extraction On the sampling day, because an adequate amount of DNA could not be obtained from the intestine of a single shrimp in the initial trial runs, three shrimp intestines from each health status for each pond were dissected on ice and pooled to compose one biological sample. Genomic DNA (gDNA) was extracted using the FAST DNA Spin kit (MO BIO Laboratories, Carlsbad, CA, USA) according to the manufacturer’s instructions. The gDNA extracts were quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies,

Appl Microbiol Biotechnol

Wilmington, USA) and then stored at − 80 °C prior to amplification.

Gut bacterial 16S rRNA gene amplification and sequencing The bacteria-specific primers 341F (5’-CCTAYGGG RBGCASCAG-3′) and 806R (5’-GGACTACNNGGGTA TCTAAT-3′) were used to amplify the V3-V4 regions of the 16S rRNA gene. To minimize PCR-induced bias, each sample was amplified in triplicate (50 μL reaction system) with the following conditions: 25 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 45 s, with a final extension at 72 °C for 10 min. To verify the expected band size, PCR products were visualized in 1.5% agarose gel. For each sample, triplicate amplicons were combined and purified using a PCR fragment purification kit. The concentrations of purified amplicons were measured using a PicoGreen-iT dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). Identical amounts of amplicons from each sample were pooled in a single tube, and were sequenced using a MiSeq platform (Illumina, San Diego, CA, USA), producing 2 × 300 bp paired-end reads.

Processing of Illumina sequencing data The paired-end reads were joined with FLASH (Magoč and Salzberg 2011) and were processed with the Quantitative Insights Into Microbial Ecology pipeline (QIIME v1.9.0) (Caporaso et al. 2010b). Briefly, the sequences with truncated or ambiguous bases at any site of more than three consecutive bases receiving a Phred quality score (Q) < 20 were deleted, as the truncated reads that had < 75% of their original length. Chimeric sequences were identified using the UCHIME algorithm (Edgar et al. 2011) and removed from the subsequent analysis. Bacteria phylotypes were identified using UCLUST (Edgar 2010) and were classified into operational taxonomic units (OTUs, 97% similarity level). The most abundant sequence of each OTU was selected as the representative sequence and then taxonomically assigned in the Greengenes database (release 13.8) (DeSantis et al. 2006) using PyNAST (Caporaso et al. 2010a). After taxonomies had been assigned, OTUs that were affiliated with Archaea, Eukaryota, Chloroplasts, and those unassigned to the bacteria domain were excluded from the dataset. The filtered alignments were then used to generate a maximum-likelihood tree using FastTree (Price et al. 2009) for phylogenetic analysis. To correct for unequal sequencing depth, we used a 20× randomly rarefied subset of 21,300 sequences per sample to calculate the distances between samples. Raw sequence data were deposited in Sequence Read Archive at DDBJ under the accession number DRA005997.

Statistical analysis A two-way (sampling date and health status) analysis of variance (ANOVA) procedure (Duncan’s Multiple Range test) was used to analyze differences in relative abundances of phyla and α-diversity of the gut microbiota over a shrimp disease progression (Churchill 2004). Principal coordinates analysis (PCoA) and analysis of similarity (ANOSIM) were conducted to evaluate the overall differences in the shrimp gut bacterial communities based on Bray-Curtis distances.

Identification of candidate pathogens and validation of their virulence genes To screen potential pathogens, gut bacterial sequences were binned into OTUs as described above, then were taxonomically classified against a local bacterial pathogens 16S rRNA gene database with 100% similarity cutoff (Chen et al. 2016). Second, SIMPER analysis was further employed to identify taxa that contribute the differences in gut communities between shrimp health status. As a result, two OTUs (V. tubiashii and V. harveyi) were searched. Third, using the profiles of the two OTUs, a Random Forest model (estimated by 10-fold cross validation) was used to evaluate the discriminative accuracy of shrimp health status. In addition, fragmentation analysis was used to test whether the two OTUs were gatekeepers in sustaining the interspecies interaction of gut microbiota. Finally, the virulence genes that were coded by the two candidate pathogens were inferred by an openaccessible virulence factors database (http://www.mgc.ac.cn/ VFs/) (Chen et al. 2005) and were validated by qPCR with corresponding specific primers (Table S3) on the 7500 Fast Real-Time PCR detection System (Applied Biosystems, Foster City, CA, USA). Specifically, 100 ng DNA was used as template in 20 μL reaction volume, in which included 10 μL SYBR Premix Ex Taq II (Takara, Dalian, China), 0. 8 μL each primer, and 0.4 μL ROX reference dye. A two-step real-time qPCR program was employed, 95 °C for 2 min, followed by 40 cycles at 95 °C for 15 s and 54 °C for 20 s. We chose three virulence genes (bca, fdeC, and tlpA genes) that are coded by the two candidate pathogens because of their important roles in facilitating pathogenic infection (Johansson and Cossart 2003). The relative abundances (virulence gene copies per 16S rRNA gene) of the three virulence genes were calculated by the comparative Ct method (Schmittgen and Livak 2008) and were compared between healthy and diseased shrimp gut microbiota, respectively.

Interspecies interaction among candidate pathogens and related resident taxa To identify gut taxa that inhibited or facilitated the establishment of the two candidate pathogens, we focused on the co-

Appl Microbiol Biotechnol

occurrence pattern among the two candidate pathogens and their directly interacted taxa. Using a stand-alone tool, MetaMIS (Metagenomic Microbial Interaction Simulator), microbial interaction network analysis was evaluated by introducing time series dataset of microbial composition at bacterial OTU-level (Shaw et al. 2016). This approach systematically examines interaction patterns, such as mutualism (+/+), competition (−/−), parasitism or predation (+/−), commensalism (+/0), amensalism (−/0), and no effect (0/0) in healthy and diseased shrimp, respectively. The topology of the resulting microbial network was visualized using Cytoscape 3.1.1 (Shannon et al. 2003).

Results Comparison of shrimp body lengths and weights between health status The average body lengths were comparable (P > 0.15 in all cases, Table S2) between healthy and diseased shrimp over the disease progression, with 8.32 ± 0.53 cm vs. 8.23 ± 0.36 cm (mean ± standard error, P = 0.239), 8.59 ± 0.37 cm vs. 8.41 ± 0.52 cm (P = 0.394), and 8.66 ± 0.34 cm vs. 8.51 ± 0.35 cm (P = 0.156) on 1st, 4th, and 10th, July, and corresponding body weights were 8.16 ± 0.35 g vs. 8.00 ± 0.41 g (P = 0.580), 8.76 ± 0.37 g vs. 8.37 ± 0.50 g (P = 0.370), and 9.20 ± 0.36 g vs. 8.90 ± 0.32 g (P = 0.193), respectively (Table S2).

Distribution of taxa in bacterial communities of healthy and diseased shrimp After quality control, a total of 1,953,763 high-quality reads and 24,966–41,845 sequences per sample (mean = 34,889 ± 4698) were obtained across the samples (N = 56, Table S2), which were clustered into 10,995 OTUs. The dominant gut bacterial phyla were γ-Proteobacteria (mean ± standard error, 54.2% ± 3.4%), followed by Actinobacteria (21.4% ± 2.9%), Bacteroidetes (17.1% ± 1.4%), Verrucomicrobia (6.3% ± 1.0%), and α-Proteobacteria (4.8% ± 0.5%) in both healthy and diseased shrimp, but significantly varied in their relative abundances (Table S4). For example, relative abundances of Verrucomicrobia and Actinobacteria dramatically decreased in diseased shrimp in relation to those in healthy ones; whereas that of γ-Proteobacteria displayed an opposite trend (Table S4).

diseased shrimp were significantly lower than the corresponding healthy shrimp over disease exacerbation (on 4th and 10th July, P < 0.01 in both cases, unpaired t test) (Fig. S1). Based on the OTUs detected across the samples, a PCoA biplot depicted obvious clusters of the gut bacterial communities between shrimp health status and over the disease progression. Notably, no clear separation was detected at the disease emergence stage, whereas the differences were apparent over diseased exacerbation (Fig. 1). These patterns were confirmed by analysis of similarity (ANOSIM), revealing that the structures of shrimp gut bacterial communities were significantly distinct (P < 0.05) between each pair of groups, but not between I701H and I701D, I704D and I710D (Table S5).

Identification of candidate pathogens Potential pathogens were firstly screened by blast against a local pathogens database (see details in the BMaterials and Methods^ section), of which V. tubiashii OTU27681 and V. harveyi OTU14490 are significantly enriched (P < 0.05) in diseased shrimp as compared with healthy ones on each sampling date (Fig. S2). SIMPER analysis showed that the two OTUs (OTU27681 and OTU14490) cumulatively contributed 64.5% dissimilarity of the gut microbiotas between healthy and diseased shrimp. Furthermore, a Random Forest model was applied to diagnose the health status of the 56 samples based on the profiles of the two OTUs, which contributed an overall 78.5% accuracy (estimated by 10-fold cross validation) (Table S6). In the healthy subjects, 20 subjects (76.9% accuracy at the 50% cutoff probability for classification) were accurately diagnosed as health. In the diseased cohorts, 24 samples (80% of diseased subjects) were correctly predicted as diseased (Table S6 and Fig. 2). Putting the pieces

Temporal variation of gut microbiota between healthy and diseased shrimp The Shannon diversity index and observed species richness among healthy shrimp were relatively stable, but linearly decreased in the diseased subjects over the disease progression (Fig. S1). In addition, bacterial richness and diversity in

Fig. 1 Principal coordinates analysis (PCoA) biplot visualized the dissimilarities of shrimp gut bacterial community among groups based on Bray-Curtis distance. Numbers represent sampling date; H healthy, D diseased

Appl Microbiol Biotechnol

Fig. 2 The predicted probabilities of shrimp disease based on profiles of the two candidate pathogens. The probability > 50% was stratified as diseased, while < 50% was stratified as healthy shrimp. Inconsistency between predicted and observed status was termed as false diagnose

together, we reasoned that the two OTUs, V. tubiashii OTU27681 and V. harveyi OTU14490, were candidate pathogens.

Fragmentation analysis of microbial co-occurrence networks To assess the degree of network fragmentation, OTU27681 and OTU14490 were removed from the gut bacterial community data. We calculated the ratio of the number of disconnected subgraphs to the overall number of nodes as proposed elsewhere (Widder et al. 2014). To ensure that the results were not subjected to a random threshold, fragmentation was calculated across a range of correlation cutoff thresholds. Based on the correlation cutoff thresholds, the default level of fragmentation of microbial networks in diseased shrimp was significantly higher than those in healthy ones (Wilcoxon rank sum test P = 0.001) (Fig. 3), indicating that removal of the two candidate pathogens contributed network fragmentation, especially in diseased shrimp. In other words, the presences of the two candidate pathogens are important for stabilizing the assembly of gut bacterial community in diseased shrimp. Thus, OTU27681 (V. tubiashii) and OTU14490 (V. harveyi) were the central phylotypes and gatekeepers in the gut microbiota of diseased shrimp.

Fig. 3 Healthy versus diseased fragmentation graph: healthy and diseased bacterial community network levels of fragmentation versus all possible constructed network thresholds. Dotted lines indicate mean fragmentation levels for Health and Disease. *Indicates a significant difference between the means of group fragmentation, Wilcoxon rank sum test W = 1426, P = 0.001

bca, fdeC, and tlpA genes that were coded by the two candidate pathogens (V. tubiashii and V. harveyi). Real-time qPCR results showed that the relative abundances of bca, fdeC, and tlpA genes in diseased shrimp increased to 2.2-, 2.9-, and 3.4folds as those in healthy cohorts, respectively (Fig. 4).

Exploring interspecies interaction of the two candidate pathogens To explore which resident taxa activate or repress the pathogens, we focused on modules of the two candidate pathogens (Fig. 5). As expected, the two candidate

Quantifying the virulence genes coded by the two candidate pathogens The virulence genes coded by the two candidate pathogens were firstly inferred using an open-accessible virulence factors database (http://www.mgc.ac.cn/VFs/) (Chen et al. 2005). To further validate the virulence, we quantified the abundances of

Fig. 4 Abundances of the three virulence genes that coded by the two candidate pathogens (V. tubiashii OTU27681 and V. harveyi OTU14490) in healthy and diseased shrimp on 10th July. Mean ± SEM and were compared using unpaired t test. *P < 0.05; **P < 0.01

Appl Microbiol Biotechnol

pathogens served as the major hubs (core ratio = 100%) in the co-occurrence network of the shrimp gut bacterial communities (Table S7). MetaMIS analysis identified 10 and 12 OTUs that interacted directly with the two candidate pathogens (i.e., OTU27681 and OTU14490) in healthy and diseased shrimp, respectively (Fig. 5). The interaction types among core OTUs exhibited distinct distribution between healthy and diseased shrimp (Fig. S3). For example, the most frequent interactive relations for OTU27681 and OTU14490 were competition (−/−) in healthy shrimp (Fig. S3A), whereas these were mutualism (+/+) in diseased subjects (Fig. S3B). Additionally, interactive strength was higher in diseased shrimp than in healthy ones (Fig. S4). In general, the interactions among the two OTUs and their resident taxa in healthy shrimp were negative, while those in diseased shrimp were positive (Fig. 5). For instance, Flavobacteriaceae OTU12579 and OTU47322 antagonized with the two candidate pathogens in healthy shrimp, which shifted into cooperative associations in diseased subjects, with concurrent increase in their abundances (Fig. 5).

Discussion

Fig. 5 Microbial interaction networks of candidate pathogens (V. tubiashii OTU27681 and V. harveyi OTU14490) in healthy a and diseased b shrimp. The solid (or dashed) arrow represents facilitative (or repressive) interaction between the two individual taxa. The sizes of

the circles were proportional to the abundance of that taxon. The two candidate pathogen OTUs were underlined. OTUs were colored at the bacterial class level

Dysbiosis in the gut microbial community is frequently coincident with host diseases (Andersen et al. 2015; Xiong et al. 2015). Given that the interspecies interaction among gut microbiota contributes pivotal roles in barrier against pathogen invasion, it is mandatory to explore how the co-occurrence pattern of potential pathogens is altered by host disease. To achieve this, we tracked the dynamics of shrimp gut bacterial communities over a disease progression. Candidate pathogens were identified and validated by their discrimination between host health status, high diagnose accuracy, gatekeeper features, and coded virulence genes. In addition, we explored how the resident taxa interact with the candidate pathogens in sustaining shrimp health. Shrimp gut microbiota was dominated by γ-Proteobacteria and Bacteroidetes species (Table S4), which is similar to the gut lineages of shrimp in previous studies (Huang et al. 2016; Xiong et al. 2015). However, this composition is distinct from that of vertebrates, such as zerbrafish (Stephens et al. 2016) and ayu (Nie et al. 2017). This divergence agrees with the host

Appl Microbiol Biotechnol

specific signatures in gut microbiota (Xiong et al. 2016). Notably, relative abundances of the dominant phyla varied significantly between shrimp health status (Table S4). It has been reported that gut Actinobacteria species can produce antibiotics to host (Protasov et al. 2017). Thus, its enrichment would benefit the host to resistant against pathogens establishment. Consistently, we found that the relative abundance of Actinobacteria was significantly higher in healthy shrimp than in diseased subjects. On the contrary, some members of γProteobacteria are known shrimp pathogens, such as Vibrio species (Austin and Zhang 2006; Vaseeharan and Ramasamy 2003), which is concordant with their enrichment in diseased individuals (Table S4). At the finer taxonomic level, the gut bacterial communities significantly differed between healthy and diseased shrimp (Fig. 1 and Table S5). Notably, alpha diversity of the gut bacterial communities was relatively stable in healthy shrimp. In contrast, the diseased shrimp exhibited a significantly lower bacterial diversity in relation to corresponding healthy subjects (Fig. S1). Consistently, it has been proposed that a reduction in diversity could lead to decreased functional stability of a bacterial community, thereby increasing the risk of developing diseases (Jones and Lennon 2010). The pathogenic bacteria 16S rRNA gene database allows us to infer potential pathogens (Chen et al. 2016). Two potential pathogens (OTU27681 V. tubiashii and OTU14490 V. harveyi) were dominant and were significantly enriched in diseased shrimp as compared with healthy ones (Fig. S2). In addition, the two OTUs cumulatively contributed 64.5% dissimilarity of the gut bacterial communities between healthy and diseased shrimp. Intriguingly, profiles of the two OTUs contributed to an overall 78.5% predicted accuracy of shrimp health status (Table S6 and Fig. 2). Indeed, OTU27681 V. tubiashii and OTU14490 V. harveyi have long been recognized as opportunistic pathogens in shrimp aquaculture (Austin and Zhang 2006; Phumkhachorn and Rattanachaikunsopon 2010). Putting these pieces together, we infer that V. tubiashii and V. harveyi are the candidate pathogens in this shrimp disease. Gatekeepers contribute important roles in sustaining the stabilization of a bacterial community, such as the fundamental organization of bacterial community network (Wagner et al. 2017; Widder et al. 2014). Therefore, a community with tightly connected species in an ecological network is more fragile by the removal of gatekeepers (Montoya et al. 2006). Theoretically, a higher level of network fragmentation means weaker co-occurrence patterns and decreased biotic interactions. Contrarily, if a network is less fragmented, then less disconnected compartments are observed (Widder et al. 2014). Consistent with this assertion, removal of V. tubiashii OTU27681 and V. harveyi OTU14490 from the co-occurrence networks resulted in a higher co-occurrence network fragmentation in diseased shrimp than in healthy ones, which supported the gatekeeper features of the two OTUs (Fig. 3). Consistently, it has repeatedly shown that selective removal

of gatekeeper species generates a high fragility of the networks (Wagner et al. 2017; Widder et al. 2014). Given that gatekeepers contribute disproportional importance in sustaining the integrity and function of a bacterial community (Widder et al. 2014), we infer that the two candidate pathogens, V. tubiashii and V. harveyi, exhibit a dominant effect on stabilizing the assembly of gut bacterial community in diseased shrimp. This feature, again, supports the causal role of the two candidate pathogens in shrimp disease. It is known that pathogenic bacteria code a variety of virulence factors to interact with host immunity and physiology, thereby causing morbidity (Wu et al. 2008). However, close relatives between commensals and pathogens share similar and even identical 16S rRNA genes (Větrovský and Baldrian 2013). Thus, we further verified the pathogenic role by quantifying the abundances of three virulence genes (bca, fedC, and tlpA) that were coded by the two candidate pathogens. As expected, the abundances of the virulence genes in diseased shrimp were significantly higher (P < 0.05 in the three cases) than in healthy ones (Fig. 4). The three detected virulence genes contribute to pathogenic infection to host. For example, the fdeC gene is involved in adhesion capability of pathogens and direct interaction with the host tissues (Nesta et al. 2012; Wu et al. 2008), while the expression of the tlpA gene can suppress the ability of TIRcontaining proteins to induce the transactivation and DNA-binding activities of NF-kB, thereby subverting normal host cell processes (Newman et al. 2006). It should be noted that we did not detect the expression level of virulence genes due to insufficient amount of mRNA can be extracted from shrimp guts. However, the expression of virulence genes is energy cost, thus pathogens only express virulence factors until it is mandatory to maximize their fitness during infection (Kitamoto et al. 2015). For this reason, a higher abundance could indicate a higher expression level. Accordingly, the increased abundances of virulence genes facilitate the ability of candidate pathogens to infection, thereby resulting in shrimp disease. Interspecies interaction plays a key role in sustaining host health (Ramayo-Caldas et al. 2016; Weng et al. 2017). Thus, we explored how resident taxa interacted with the candidate pathogens. The interactions among resident taxa and the two candidate pathogens (i.e., V. tubiashii OTU27681 and V. harveyi OTU14490) were generally negative in healthy shrimp, which shifted into positive interactions in the diseased subjects (Fig. 5). Seven out of 17 OTUs were affiliated with Vibrio genus, whose relative abundances in diseased shrimp were higher than in healthy shrimp (Fig. 5). The two candidate pathogens connected most interaction types (mutualism (+/+), competition (−/−), parasitism or predation (+/−) and amensalism (−/0)) to the species in the microbial interaction networks (Fig. S3), which is concurrent with their underlying

Appl Microbiol Biotechnol

biological connections in the consensus network (Fig. S4). Because amensalism (−/0) interaction was displayed, thus the two candidate pathogens were inhibited by the resident OTUs in healthy individuals. It has been proposed that a decrease in abundance or loss of keystone species can destabilize the community (Coyte et al. 2015; Wagner et al. 2017). Consistent with this assertion, Rhodobacteraceae OTU12580 was abundant in healthy shrimp, while it was lost or undetectable in diseased ones (Fig. 5). It is known that Rhodobacteraceae species are dominant members in healthy shrimp gut, which contribute important roles in carbon biogeochemical cycling. In practical, Rhodobacteraceae strains have been widely applied in aquaculture with a high probiotics effects (Newaj-Fyzul et al. 2014; Xiong et al. 2014). In this regard, the loss of Rhodobacteraceae OTU12580 may unbalance the microbial loop, thereby causing shrimp disease. Notably, Flavobacteriaceae species (OTU12579, OTU26422, and OTU47322) and Photobacterium OTU6879 (affiliated to Vibrionaceae) were negatively associated with two V. tubiashii OTU27681 and V. harveyi OTU14490 in healthy shrimp, whereas this pattern shifted into corporative interactions in disease cohorts, with increased relative abundances of these OTUs (Fig. 5). Commensal Flavobacteriaceae species occupy unique nutritional and ecological niches, thereby preventing the growth of pathogens in healthy individuals. However, pathogens can use its virulence factors to localize near/at the intestinal epithelium, which in turn benefit from the bioavailable substrates produced by Flavobacteriaceae (Kitamoto et al. 2015). Consistent with this assertion, abundances of the three virulence genes significantly increased in the diseased shrimp (Fig. 4). In addition, pathogen infection generally induces host inflammatory responses, which in turn provides a growth advantage to pathogens and their closely phylogenetic relatives (share similar niche preference) (Chen et al. 2017; Stecher et al. 2010). Thus, it is expected there were increased abundances of resident Vibrio species and Photobacterium (Fig. 5). In summary, two candidate pathogens were identified from an ecological perspective, which was validated by their enriched abundances, high discriminative accuracy between shrimp health status, increased abundances of their coded virulence genes, and gatekeeper roles in maintaining interspecies interaction in diseased shrimp. Going a step further, we identified taxa that confer resistance or facilitation to the establishment of candidate pathogens, which in turn, could guide a rational design of microbiota-based diagnose and therapy. Collectively, our study offers a new frame to identify potential polymicrobial pathogens, and provides novel insights on the pathogenic mechanism from an ecological perspective. Acknowledgements This work was supported by the Project of Science and Technology Department of Ningbo (2017C10044), the Zhejiang

Province Public Welfare Technology Application Research Project (2016C32063), and the K.C. Wong Magna Fund in Ningbo University.

Compliance with ethical standards This article does not contain any studies with human participants performed by any of the authors. All applicable international, national, and/or institutional guidelines for the care and use of animals were followed the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Conflict of interest The authors declare that they have no conflict of interest.

References Alou MT, Million M, Traore SI, Mouelhi D, Khelaifia S, Bachar D, Caputo A, Delerce J, Brah S, Alhousseini D, Sokhna C, Robert C, Diallo BA, Diallo A, Parola P, Golden M, Lagier JC, Raoult D (2017) Gut bacteria missing in severe acute malnutrition, can we identify potential probiotics by culturomics? Front Microbiol 8:899 Andersen SB, Marvig RL, Molin S, Johansen HK, Griffin AS (2015) Long-term social dynamics drive loss of function in pathogenic bacteria. Proc Natl Acad Sci U S A 112:10756–10761 Austin B, Zhang XH (2006) Vibrio harveyi: a significant pathogen of marine vertebrates and invertebrates. Lett Appl Microbiol 43:119– 124 Benson AK, Kelly SA, Legge R, Ma F, Low SJ, Kim J, Zhang M, Oh PL, Nehrenberg D, Hua K, kachman SD, Moriyama EN, Walter J, Peterson DA, Pomp D (2010) Individuality in gut microbiota composition is a complex polygenic trait shaped by multiple environmental and host genetic factors. Proc Natl Acad Sci U S A 107: 18933–18938 Boursier J, Mueller O, Barret M, Machado M, Fizanne L, Araujo-Perez F, Guy CD, Seed PC, Rawls JF, David LA, Hunault G, Oberti F, Calès P, Diehl AM (2016) The severity of nonalcoholic fatty liver disease is associated with gut dysbiosis and shift in the metabolic function of the gut microbiota. Hepatology 63:764–775 Buffie CG, Bucci V, Stein RR, McKenney PT, Ling L, Gobourne A, No D, Liu H, Kinnebrew M, Viale A, Littmann E, van den Brink MRM, Jenq RR, Taur Y, Sander C, Cross J, Toussaint NC, Xavier JB, Pamer EG (2015) Precision microbiome restoration of bile acidmediated resistance to Clostridium difficile. Nature 517:205–208 Byrd AL, Segre JA (2016) Adapting Koch's postulates. Science 351:224– 226 Caporaso JG, Bittinger K, Bushman FD, Desantis TZ, Andersen GL, Knight R (2010a) PyNAST: a flexible tool for aligning sequences to a template alignment. Bioinformatics 26:266–267 Caporaso JG, Kuczynski J, Stombaugh J (2010b) QIIME allows analysis of high-throughput community sequencing data. Nat Methods 7: 335–336 Carding S, Verbeke K, Vipond DT, Corfe BM, Owen LJ (2015) Dysbiosis of the gut microbiota in disease. Microb Ecol Health Dis 26:26191 Carraro L, Maifreni M, Bartolomeoli I, Martino ME, Novelli E, Frigo F, Marino M, Cardazzo B (2011) Comparison of culture-dependent and -independent methods for bacterial community monitoring during Montasio cheese manufacturing. Res Microbiol 162:231–239 Chen L, Yang J, Yu J, Yao Z, Sun L, Shen Y, Jin Q (2005) VFDB: a reference database for bacterial virulence factors. Nucleic Acids Res 33:325–328 Chen Q, An X, Li H, Su J, Ma Y, Zhu YG (2016) Long-term field application of sewage sludge increases the abundance of antibiotic resistance genes in soil. Environ Int 92:1–10

Appl Microbiol Biotechnol Chen WY, Ng TH, Wu JH, Chen JW, Wang HC (2017) Microbiome dynamics in a shrimp grow-out pond with possible outbreak of acute hepatopancreatic necrosis disease. Sci Rep 7:9395 Churchill GA (2004) Using ANOVA to analyze microarray data. BioTechniques 37:173–175 Conejero MJ, Hedreyda CT (2004) PCR detection of hemolysin (vhh) gene in Vibrio harveyi. J Gen Appl Microbiol 50:137–142 Coyte KZ, Schluter J, Foster KR (2015) The ecology of the microbiome: networks, competition, and stability. Science 350:663–666 Dai W, Yu W, Zhang J, Zhu J, Tao Z, Xiong J (2017) The gut eukaryotic microbiota influences the growth performance among cohabitating shrimp. Appl Microbiol Biotechnol 101(16):6447–6457 Desantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, Huber T, Dalevi D, Hu P, Andersen GL (2006) Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 72:5069–5072 Edgar RC (2010) Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26:2460–2461 Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R (2011) UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27:2194–2200 Gauger EJ, Gómezchiarri M (2002) 16S ribosomal DNA sequencing confirms the synonymy of Vibrio harveyi and V. carchariae. Dis Aquat Org 52:39–46 Huang Z, Li X, Wang L, Shao Z (2016) Changes in the intestinal bacterial community during the growth of white shrimp, Litopenaeus vannamei. Aquac Res 47:1737–1746 Johansson J, Cossart P (2003) RNA-mediated control of virulence gene expression in bacterial pathogens. Trends Microbiol 11:280–285 Jones SE, Lennon JT (2010) Dormancy contributes to the maintenance of microbial diversity. Proc Natl Acad Sci U S A 107:5881–5886 Kamada N, Chen GY, Inohara N, Núñez G (2013) Control of pathogens and pathobionts by the gut microbiota. Nat Immunol 14:685–690 Kasubuchi M, Hasegawa S, Hiramatsu T, Ichimura A, Kimura I (2015) Dietary gut microbial metabolites, short-chain fatty acids, and host metabolic regulation. Nutrients 7:2839–2849 Kitamoto S, Nagaokitamoto H, Kuffa P, Kamada N (2015) Regulation of virulence: the rise and fall of gastrointestinal pathogens. J Gastroenterol 51:195–205 Krajmalnikbrown R, Ilhan ZE, Kang DW, Dibaise JK (2012) Effects of gut microbes on nutrient absorption and energy regulation. Nutr Clin Pract 27:201–214 Luis BBS, Hedreyda CT (2006) Analysis of a gene (vch) encoding hemolysin isolated and sequenced from Vibrio campbellii. J Gen Appl Microbiol 52:303–313 Magoč T, Salzberg SL (2011) FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27:2957–2963 Mallon CA, Elsas JD, Salles JF (2015) Microbial invasions: the process, patterns, and mechanisms. Trends Microbiol 23:719–729 Montoya JM, Pimm SL, Solé RV (2006) Ecological networks and their fragility. Nature 442:259–264 Ndoye B, Rasolofo EA, LaPointe G, Roy D (2011) A review of the molecular approaches to investigate the diversity and activity of cheese microbiota. Dairy Sci Technol 91:495–524 Nesta B, Spraggon G, Alteri C, Moriel DG, Rosini R, Veggi D, Smith S, Bertoldi I, Pastorello I, Ferlenghi I, Fontana MR, Frankel G, Mobley HLT, Rappuoli R, Pizza M, Serino L, Soriani M (2012) FdeC, a novel broadly conserved Escherichia coli adhesin eliciting protection against urinary tract infections. MBio 3:e00010–e00012 Newaj-Fyzul A, Al-Harbi A, Austin B (2014) Developments in the use of probiotics for disease control in aquaculture. Aquaculture 431:1–11 Newman RM, Salunkhe P, Godzik A, Reed JC (2006) Identification and characterization of a novel bacterial virulence factor that shares homology with mammalian Toll/Interleukin-1 receptor family proteins. Infect Immun 74:594–601

Nie L, Zhou Q-J, Yan Q, Chen J (2017) Interplay between the gut microbiota and immune responses of ayu (Plecoglossus altivelis) during Vibrio anguillarum infection. Fish Shellfish Immunol 68:479–487 Peters BM, Jabra-Rizk MA, Graeme A, Costerton JW, Shirtliff ME (2012) Polymicrobial interactions: impact on pathogenesis and human disease. Clin Microbiol Rev 25:193–213 Phumkhachorn P, Rattanachaikunsopon P (2010) Isolation and partial characterization of a bacteriophage infecting the shrimp pathogen Vibrio harveyi. Afr J Microbiol Res 4:1794–1800 Price MN, Dehal PS, Arkin AP (2009) FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol 26:1641–1650 Protasov ES, Axenov-Gribanov DV, Rebets YV, Voytsekhovskaya IV, Tokovenko BT, Shatilina ZM, Luzhetskyy AN, Timofeyev MA (2017) The diversity and antibiotic properties of Actinobacteria associated with endemic deepwater amphipods of Lake Baikal. Anton Leeuw Int J G 5:1–19 Qiao G, Xu H, Wang Z, Jang IK, Qi Z, Zhang M (2015) Comparison of immune response of Pacific white shrimp, Litopenaeus vannamei, after multiple and single infections with WSSV and Vibrio anguillarum. Fish Shellfish Immunol 44(1):257–264 Ramayo-Caldas Y, Mach N, Lepage P, Levenez F, Denis C, Lemonnier G, Leplat JJ, Billon Y, Berri M, Doré J, Rogel-Gaillard C, Estellé J (2016) Phylogenetic network analysis applied to pig gut microbiota identifies an ecosystem structure linked with growth traits. ISME J 10:2973–2977 Reddy A, Jeyasekaran G, Shakila R (2013) Morphogenesis, pathogenesis, detection and transmission risks of white spot syndrome virus in shrimps. Fish Aquacult J 3:1–13 Sawabe T, Kita-Tsukamoto K, Thompson FL (2007) Inferring the evolutionary history of vibrios by means of multilocus sequence analysis. J Bacteriol 189:7932–7936 Schmittgen TD, Livak KJ (2008) Analyzing real-time PCR data by the comparative CT method. Nat Protoc 3:1101–1108 Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13:2498–2504 Shaw GTW, Pao YY, Wang D (2016) MetaMIS: a metagenomic microbial interaction simulator based on microbial community profiles. BMC Bioinformatics 17:488 Sriurairatana S, Boonyawiwat V, Gangnonngiw W, Laosutthipong C, Hiranchan J, Flegel TW (2014) White feces syndrome of shrimp arises from transformation, sloughing and aggregation of hepatopancreatic microvilli into vermiform bodies superficially resembling gregarines. PLoS One 9(6):e99170 Stecher B, Chaffron S, Käppeli R, Hapfelmeier S, Freedrich S, Weber TC, Kirundi J, Suar M, Mccoy KD, Mering CV, Macpherson A, Hardt W (2010) Like will to like: abundances of closely related species can predict susceptibility to intestinal colonization by pathogenic and commensal bacteria. PLoS Pathog 6:e1000711 Stephens WZ, Burns AR, Stagaman K, Wong S, Rawls JF, Guillemin K, Bohannan BJ (2016) The composition of the zebrafish intestinal microbial community varies across development. ISME J 10:644– 654 Sun J, Zeng Q, Tsang DC, Zhu L, Li X (2017) Antibiotics in the agricultural soils from the Yangtze River Delta, China. Chemosphere 189: 301–308 Ussar S, Griffin NW, Bezy O, Fujisaka S, Vienberg S, Softic S, Deng L, Bry L, Gordon JI, Kahn CR (2015) Interactions between gut microbiota, host genetics and diet modulate the predisposition to obesity and metabolic syndrome. Cell Metab 22:516–530 Vaseeharan B, Ramasamy P (2003) Control of pathogenic Vibrio spp. by Bacillus subtilis BT23, a possible probiotic treatment for black tiger shrimp Penaeus monodon. Lett Appl Microbiol 36:83–87

Appl Microbiol Biotechnol Větrovský T, Baldrian P (2013) The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS One 8:e57923 Wagner MB, Waite DW, Hoggard M, Douglas RG, Taylor MW, Biswas K (2017) Bacterial community collapse: a meta-analysis of the sinonasal microbiota in chronic rhinosinusitis. Environ Microbiol 19:381–392 Weng FCH, Shaw GTW, Weng CY, Yang YJ, Wang D (2017) Inferring microbial interactions in the gut of the Hong Kong whipping frog (Polypedates megacephalus) and a validation using probiotics. Front Microbiol 8:525 Widder S, Besemer K, Singer GA, Ceola S, Bertuzzo E, Quince C, Sloan WT, Rinaldo A, Battin TJ (2014) Fluvial network organization imprints on microbial co-occurrence networks. Proc Natl Acad Sci U S A 111:12799–12804 Wu HJ, Wang AH, Jennings MP (2008) Discovery of virulence factors of pathogenic bacteria. Curr Opin Chem Biol 12:93–101 Xiong J, Zhu J, Zhang D (2014) The application of bacterial indicator phylotypes to predict shrimp health status. Appl Microbiol Biotechnol 98:8291–8299

Xiong J, Wang K, Wu J, Qiuqian L, Yang K, Qian Y, Zhang D (2015) Changes in intestinal bacterial communities are closely associated with shrimp disease severity. Appl Microbiol Biotechnol 99:6911– 6919 Xiong J, Dai W, Li C (2016) Advances, challenges, and directions in shrimp disease control: the guidelines from an ecological perspective. Appl Microbiol Biotechnol 100:6947–6954 Xiong J, Dai W, Zhu J, Liu K, Dong C, Qiu Q (2017a) The underlying ecological processes of gut microbiota among cohabitating retarded, overgrown and normal shrimp. Microb Ecol 73(4):988–999 Xiong J, Zhu J, Dai W, Dong C, Qiu Q, Li C (2017b) Integrating gut microbiota immaturity and disease-discriminatory taxa to diagnose the initiation and severity of shrimp disease. Environ Microbiol 19: 1490–1501 Zhu J, Dai W, Qiu Q, Dong C, Zhang J, Xiong J (2016) Contrasting ecological processes and functional compositions between intestinal bacterial community in healthy and diseased shrimp. Microb Ecol 72:975–985 Zhuang F, Li H, Zhou X, Zhu Y, Su J (2017) Quantitative detection of fecal contamination with domestic poultry feces in environments in China. AMB Express 7:80