1 INTRODUCTION

Vibrio harveyi is a facultatively anaerobic, bioluminescent, and Gram-negative rod-shaped bacterium (Johnson and Shunk, 1936). Bacterial infection causes serious vibriosis in marine fish and invertebrates (Austin and Zhang, 2006), particularly in juvenile populations (Pujalte et al., 2003), leading to massive deaths and major economic losses in aquaculture. V. harveyi has been isolated from multiple sources including common snook (Centropomus undecimalis) in USA (Kraxberger-Beatty et al., 1990), silvery black porgy (Acanthopagrus cuvieri) and cultured brown spotted grouper (Epinepelus tauvina) in Kuwait (Saeed, 1995), and sunfish (Mola mola) in Spain (Hispano et al., 1997). V. harveyi-associated diseases are major constraint on the production of marine invertebrates, particularly in South America and Asia. V. harveyi infections have also been reported in common dentex (Dentex dentex L.) (Company et al., 1999), salmonids (Salmo salar L.) (Zhang and Austin, 2000), seahorse (Hippocampus kuda) (Alcaide et al., 2010), farmed sole (Solea senegalensis) (Zorrilla et al., 2003), summer flounder (Paralichthys dentatus) (Gauger et al., 2006), and gilthead sea bream (Sparus aurata) (Haldar et al., 2010). In addition, V. harveyi has been isolated from several crustacean larvae, including black tiger prawn (Penaeus monodon) (Lavilla-Pitogo et al., 1990), kuruma prawns (Penaeus japonicas) (Liu et al., 1996), rock lobster (Jasus verreauxi) (Diggles et al., 2000), and oysters (Ortigosa et al., 1994). Infected animals became inappetent and developed inflammation and necrotic subdermal cysts in various species of fishes, such as crustaceans and mollusks (Austin, 2010).

Bacteria could enter viable but non-culturable (VBNC) or normal dormant state by a variety of environmental stresses (Wolf and Oliver, 1992; Biosca et al., 1996; Jiang and Chai, 1996). Moreover, nutrient deficiency is the most common environmental stress in natural ecosystems for microorganisms. It is reported that some species could survive for a long time in response to starvation condition such as V. harveyi, V. cholerae, V. fischeri, V. anguillarum, V. campbellii, V. mimicus, V. proteolyticus, V. vulnificus (Jiang and Chai, 1996; McDougald et al., 1998). Like other bacterial species, V. harveyi may be chronically exposed to nutrient-poor waters and its ability to survive in water for extended periods of time is an adaptive strategy against stress. Imaging showed that starved cells alter their morphology from rod to spherule cells (Novitsky and Morita, 1976; Sun et al., 2016). Starved cells will adjust their gene expression to adapt to environmental stress, with potential changes to transcriptional regulators, heat shock proteins, virulence regulators, and genes related to metabolism. However, little is known about gene expression profiles under starvation stress.

Hence, we obtained whole transcriptomic profiles of V. harveyi DY1 strain under starvation stress by high-throughput sequencing and determined the gene expression changes compared to non-starved cells. This is the first comprehensive transcriptomic analysis of starvation effects in pathogenic V. harveyi, and the results greatly strengthen our understanding of the molecular pathogenic and transmission mechanisms of V. harveyi strain DY1.

2 MATERIAL AND METHOD 2.1 Bacterial strain

The V. harveyi DY1 strain was first isolated from cultured megalopa of swimming crab (Portunus trituberculatus) in 2011 by Zhang et al. (2014) in Jiangsu, China. DY1 strain was recovered from -80℃ as glycerol stocks and incubated with LB broth (Hopebio Qingdao, China) at 28℃ overnight with shaking.

2.2 Starvation stress and cell enumeration

Overnight cultures were diluted (1:100) with fresh LB broth and grown to the mid-log phase at 28℃. Cultures were collected by centrifugation at 12 000×g for 10 min. The pellets were washed twice with sterile artificial seawater to rinse nutrients off and suspended in the starvation regimes. Then the cells in triplicate were placed in Erlenmeyer flasks containing 100 mL sterile artificial seawater to give a final concentration of 10⁸ CFU/mL. Cells were incubated in the dark at 25℃ and monitored for a period of four weeks according to the previous method. Cells in starvation regimes were sampled at 1, 7, 14, 21, and 28 days and counted by plate count method (Sun et al., 2016). The cell number was counted and converted to base-10 logarithms to fit normal distribution model.

2.3 Morphology analysis

The morphology changes of V. harveyi DY1 strain after four weeks of starvation was observed by SEM as previously described (Arias et al., 2012). The nonstarved and starved cells were centrifuged at 10 000×g for 10 min at 4℃. The pellets were washed with PBS, and fixed with glutaraldehyde (0.25 g/L). A graded ethanol series was used to dehydrate the samples, then coated with gold palladium alloy in an Electron Microscopy Science (EMS 550X), and examined by Zeiss EVO 50 (Zeiss, Germany).

2.4 Library construction and Illumina sequencing

Cells starved for four weeks and the non-starved cells were subjected to total RNA isolation and cDNA preparations. Three independent trials were performed for sequencing. Total RNA was extracted from cells using the EasyPure RNA kit (TransGen Biotech, Beijing, China). rRNA was removed after total RNA was collected from prokaryocyte. Fragmentation buffer was added for splitting mRNA to short fragments. Taking these short fragments as templates, random hexamer-primer were used to synthesize the first-strand cDNA and the second-strand cDNA was synthesized using buffer, RNase H, dNTPs, DNA polymerase I. Short fragments were purified with QiaQuick PCR extraction kit and resolved with elution buffer for end reparation and adding poly (A). After that, the short fragments were connected with sequencing adapters and the UNG enzyme was used to degrade the second-strand cDNA, and the product was purified by MiniElute PCR purification kit before PCR amplification. At last, the library was sequenced using Illumina HiSeq2000. Images generated by sequencers were converted by base calling into nucleotide sequences, which were called raw data or raw reads and are stored in FASTQ format. Raw reads produced from sequencing machines contain dirty reads that contain adapters, unknown or low quality bases. After removing reads with adaptors and low quality reads, clean reads were obtained. The clean reads were mapped to reference genome and genes sequences respectively using SOAP2 (Li et al., 2009).

2.5 Gene annotation and difference analysis

After data filtering, we compared the clean reads with the reference genome and sequence which has been deposited at DDBJ/EMBL/GenBank under the accession numbers CP009467.2 and CP009468.1. Gene coverage is the percentage of a gene covered by reads. This value equals to ratio of the number of bases in a gene covered by unique mapping reads to number of total bases in that gene. The threshold of the P-value in multiple tests was determined by setting FDR (false discovery rate) at 0.000 1. The FPKM (fragments per kb per Million reads) method was applied to calculate gene expression (Mortazavi et al., 2008). To screen the DEGs in DY1 cells potentially related to starvation stress, statistical and analysis of gene expression was applied using fold change > 2 and P≤0.001 as standards.

2.6 Functional annotation and GO/KEGG enrichment analysis

For gene ontology analysis, the alignment results were parsed for assigning GO terms by using Blast2GO software. Gene sequences of V. harveyi DY1 strain were aligned against Gene Ontology database. The calculated P-value went through Bonferroni Correction, taking corrected P≤0.05 as a threshold. GO terms fulfilling this condition are defined as significantly enriched GO terms in DEGs.

KEGG pathway enrichment analysis identifies significantly enriched metabolic pathways or signal transduction pathways in DEGs comparing with the whole genome background. The calculating formula is the same as that in GO analysis.

2.7 Quantitative real-time PCR analysis

RNA were extracted using an EasyPure RNA kit (TransGen Biotech, Beijing, China) and reversetranscribed by TransScript One-Step gDNA Removal and cDNA Synthesis Supermix (TransGen Biotech, Beijing, China) according to the manufacturer's instructions. The cDNA was synthesized using anchored oligo (dT) 18 primer and incubated for 15 min at 42℃. Real-time PCR was performed with 10 μL SYBR^® SuperMix, 1 μL diluted cDNA, 0.4 μL of forward primer, 0.4 μL of reverse primer and 8.2 μL of nuclease-free water. Reactions were incubated at 95℃ for 10 min, followed by 40 cycles of 95℃ for 10 s, 60℃ for 30 s. Each sample was run in triplicate for analysis. The relative gene expression were calculated using 2^-ΔΔCt method and16S rRNA gene was chosen as an internal control for normalization. Specific primers were designed according to the corresponding sequences of V. harveyi (Table 1).

2.8 Statistics analysis

The data are presented as the mean±SD (n=3). One-Way Analysis of variance (one-way ANOVA) was used to evaluate the differential expression by SPSS (20.0).

3 RESULT 3.1 Morphology analysis

Vibrio harveyi was found decreased in size and changed in shape from a rod to a sphere after four weeks of starvation shown in SEM. There was also a significant decrease in size of cells after starvation compared to the initial population (Fig. 1).

3.2 De novo assembly and gene annotation

In this study, RNA-Sequencing was used to study the adaptation of V. harveyi under four-week starvation using Illumina HisSeq^TM 2000 system. After filtering the raw data, 14 422 512 and 8 662 506 clean reads were obtained in control and starvation stress group (Table 2). Proportions of clean reads mapped back to genome and genes can provide an overall assessment of the sequencing. The reads were aligned to V. harveyi reference genome using the SOAPaligner/soap2 software and mismatches no more than 5 bases were allowed in the alignment. There were 12 375 624 (85.81%) and 6 846 716 (79.04%) clean reads were mapped to reference genome in control and starvation stress group (Table 2). There were 5 605 402 (38.87%) and 3 477 390 (40.14%) reads and ratio of perfect matches in total mapped reads in the control and starvation stress group (Table 2). There were 5 337 genes annotated in the transcriptome analysis. V. harveyi strain DY1 revealed 4 679 and 4 661 annotated genes before and after starvation stress, respectively.

3.3 Gene expression difference analysis

A total of 666 differentially expressed genes (DEGs) were identified in the samples from the fourweek starved group compared to the samples from the untreated group, including 411 down-regulated genes and 255 up-regulated genes (Fig. 2). Thus, about 14.2% of the total number of genes showed significantly altered expression levels after starvation stress for four weeks. The up-regulated DEGs included 20 genes with more than 10-fold expression change between the starved group and the untreated group. Similarly, of the down-regulated DEGs, there were four genes with more than 10-fold expression change.

3.4 GO functional classification of DEGs

GO analyses were performed to classify gene functions of V. harveyi under starvation stress. GO analysis of DEGs in V. harveyi after starvation revealed that in the biological process catalogue, the most abundant category of the differentially expressed genes control metabolism. There were 136 upregulated genes assigned metabolic process, such as hydroxyethylthiazole kinase and peptidase. There were 170 downregulated genes assigned metabolic process including hemolysin, oxidoreductase, 6-phosphogluconate dehydrogenase, et al. In cellular component, the top three categories were cell (55 upregulated, 88 down regulated), cell part (53 upregulated, 86 down regulated) and membrane process (51 upregulated, 64 down regulated). In the molecular function catalogue, the most abundant terms of the differential expression gene were catalytic activity, followed by binding and transporter activity process, primarily including cell division protein ZapB, flagellar motor switch protein FliM, sulfur transfer protein TusE and transcriptional regulator (Fig. 3).

3.5 KEGG enrichment analysis of DEGs

Expression of some genes plays a significant role in protection against environmental stress. KEGG enrichment analysis was used to identify the most enriched pathways in V. harveyi after starvation stress. In KEGG enrichment, the DEGs were assigned into 120 pathways. The most abundant pathways were metabolic pathways, microbial metabolism in diverse environments and Biosynthesis of secondary metabolites. Flagellar assembly pathway and Butirosin and neomycin biosynthesis pathway were the most representative pathways related to starvation stress (Fig. 4). Moreover, ABC transporters pathway, flagellar assembly pathway, fatty acid metabolism pathway, fatty acid degradation pathway, biosynthesis of antibiotics pathway and beta-Lactam resistance pathway also involved many differentially expressed genes.

3.6 Verification of the DEGs by qRT-PCR

To validate the sequencing results, eight differentially expressed genes (danK, fimT, hcp, impA, aphB, rnfA, fur, htpG) of V. harveyi were randomly selected for qRT-PCR. Among these genes, four genes (danK, fimT, impA, aphB) are up-regulated genes and the other four (hcp, rnfA, fur, htpG) are down-regulated genes. The cell samples that subjected for qRT-PCR are the same as those used for transcriptome sequencing. The expression profile of these DEGs are similar to those obtained in sequencing results (Fig. 5).

4 DISCUSSION

High mortality rates in the early stages of the marine fish and invertebrates caused by V. harveyi is one of the most important reasons for economic losses during the production period. However, the survival strategies of V. harveyi in the aquatic environment are not well understood. Microorganisms must adapt to environmental stress to grow and survive (Poindexter, 1981; Kjellerberg et al., 1983; Kunttu et al., 2009). Previous studies showed that the most widespread environmental stress for microbes in a natural environment is a lack of nutrients (Wai et al., 1999; Vatsos et al., 2003; Suzina et al., 2004; Montánchez et al., 2014; Kaberdin et al., 2015; Parada et al., 2016; Montánchez et al., 2019). To study the physiological response and the long-term survival mechanism upon starvation stress is urgent and important.

In this study, we compared the size of the starved cells with cells from the control group, and observed that the starved cells were obviously shorter (Fig. 1). This implied that V. harveyi likely triggered specific regulatory mechanisms to adjust its shape and size. This result was previously reported (Sun et al., 2016). Parada et al. (2016) found that morphological change was not directly related to the entry of cells into the VBNC state, although they observed a gradual decrease in the size of V. harveyi cells during incubation. Montánchez et al. (2014) found that incubation in cold seawater for 12 h did not cause any significant morphological changes in V. harveyi, suggesting that V. harveyi likely elicits specific adaptation mechanisms maintaining its culturability under stress conditions.

The reads were mapped to the reference genome. The gene coverage was above 70% for both groups. The RNA-seq data can be used to assess the variation in expression of virulence-related and metabolismrelated factors in V. harveyi, which suggest the strategies used by these bacteria to survive under starvation conditions. There were 255 up-regulated genes and 411 down-regulated genes under starvation stress, and these genes very likely act in the response to starvation stress (Fig. 2). The analysis shown in Fig. 3 revealed 35 GO categories relevant to starvation stress of V. harveyi. The main terms associated with starvation stress were binding and transporter activity, metabolic process, and cell and catalytic activity. Previous proteomic study revealed that the level of membrane proteins participating in cellular transport, maintenance of cell structure, and bioenergetics processes remained unchanged during starvation at low temperature, suggesting that V. harveyi might need these proteins for long-term survival or for the resuscitation process after dormancy (Parada et al., 2016).

The analysis of KEGG enrichment may help us to understand the molecular mechanisms utilized by V. harveyi for long-term survival under starvation stress. The most abundant KEGG pathways included antibiotic biosynthesis, ABC transporters, butanoate metabolism pathway, starch and sucrose metabolism, and beta lactam resistance pathways, which are processes that have previously been related to starvation stress (Higgins, 2001; Kim et al., 2013; Svensson et al., 2014). Genes related to nutrient transport, such as phospholipids, oligopeptides, and simple sugars, were down-regulated after four weeks starvation. However, genes related to nutrient transfer, such as the phospholipid transfer protein gene and the monosaccharide transporter gene, were significantly up-regulated under starvation stress, representing an active metabolism pathway in V. harveyi. This change may be supplementation sparked by the compensatory response, and changes in the expression of genes controlling biosynthesis of lipids and molecular transport likely affect the composition and properties of the V. harveyi cell envelope, an obvious adaptation to stress. We suspect that V. harveyi comprehensively utilized multiple strategies and adaptation mechanisms to sustain key physiological functions under nutritional stress. At least in part, V. harveyi compensated for the reduced expression of biosynthetic genes by the upregulation of transporter genes controlling the uptake of amino acids (e.g. amino acid ABC transporters). Previous research found that limitation of nutrients leads to the significant downregulation of genes controlling central carbon metabolism, biosynthesis of lipids, amino acids, and nucleotides (Montánchez et al., 2014).

To better understand the survival strategies of V. harveyi under starvation stress, the fatty acid metabolism pathway was analyzed. The genes related with fatty acid synthesis were significantly downregulated and the expression level of fatty acid degradation related genes were significantly increased, producing a large number of acetyl coenzyme A for TCA cycle. The observation of the up-regulation of fatty acid degradation correlates to previous work by Kaberdin et al. (2015). Moreover, the fine-tuning of these metabolic pathways appears to be attained through the action of small regulatory RNAs known for their essential roles in post-transcriptional control of gene expression (Kaberdin and Bläsi, 2006).

Virulence genes played a crucial role in environmental adaptation and virulence. Previous study analyzed V. harveyi adaptation in sea water microcosms at elevated temperature and found that elevated temperature also affected regulation of V. harveyi genes controlling its virulence and ancillary mechanisms (i.e. production and secretion of virulence factors, biofilm formation and motility) (Montánchez et al., 2019). The sequencing analysis revealed downregulated expression levels of virulence regulator gene toxR and flagellin A. These virulence genes regulate a considerable number of genes involved in environmental adaptation and virulence. Other down-regulated genes include flagellinassociated genes, such as flgB, flgC, flgD, flgF, flgG, flgH, flgI, flgK, flgL, fliM, and fliO. Down-regulation of flagellin-associated genes will directly reduce flagellar synthesis of V. harveyi, resulting in reduced motility and adhesive strength, which may be critical during the initials steps of infection. Thus, the pathogenicity may decrease under starvation stress. The decreased expression of flagellin-associated genes suggests that under starvation conditions, V. harveyi reduced energy requirements by reducing the synthesis of non-essential structures.

Under starvation stress, many genes were enriched in metabolic pathways involved in the synthesis of macrolides and ketolide antibiotics (Fig. 4). The sequencing revealed up-regulated expression levels of many genes related to the synthesis of macrolides and ketolides and beta lactam resistance genes, including TolC (outer membrane protein), acrA (membrane fusion protein, multidrug efflux system), and opp (oligopeptide transport system). These resistance-related genes may significantly improve the survival rate of V. harveyi.

5 CONCLUSION

This study presented the transcriptomic response profiles of the Portunus trituberculatus pathogen V. harveyi strain DY1 under starvation conditions, to elucidate the changes in V. harvey gene expression due to starvation in seawater. V. harveyi may survive under starvation conditions by regulating the expression of virulence and metabolism-related genes. The identified genes may be important highvalue drug targets, suggesting new ways to effectively control clinical infection of V. harveyi. Our findings should facilitate future study on stress adaptation and environmental survival mechanisms of V. harveyi.

6 DATA AVAILABILITY STATEMENT

All sequence data that support the findings of this study have been deposited in the NCBI Short Read. The sequence read archive (SRA) accession number: PRJNA507871.