Transcriptomic analysis on the formation of the viable putative non-culturable state of beer-spoilage Lactobacillus acetotolerans

Lactic acid bacteria (LAB) are the most common beer-spoilage bacteria regardless of beer type, and thus pose significant problems for the brewery industry. The aim of this study was to investigate the genetic mechanisms involved in the ability of the hard-to-culture beer-spoilage bacterium Lactobacillus acetotolerans to enter into the viable putative non-culturable (VPNC) state. A genome-wide transcriptional analysis of beer-spoilage L. acetotolerans strains BM-LA14526, BM-LA14527, and BM-LA14528 under normal, mid-term and VPNC states were performed using RNA-sequencing (RNA-seq) and further bioinformatics analyses. GO function, COG category, and KEGG pathway enrichment analysis were conducted to investigate functional and related metabolic pathways of the differentially expressed genes. Functional and pathway enrichment analysis indicated that heightened stress response and reduction in genes associated with transport, metabolic process, and enzyme activity might play important roles in the formation of the VPNC state. This is the first transcriptomic analysis on the formation of the VPNC state of beer spoilage L. acetotolerans.


Results
Formation of the VPNC state. Cellular viability and culturability of L. acetotolerans strains were evaluated once a week, and the formation of VPNC state by L. acetotolerans was obtained and verified after 16 ± 0.8 subcultures (112 ± 5.6 days).
Overview of the transcriptional analysis. cDNA libraries of nine L. acetotolerans samples were constructed, sequenced and generated with a total of 15,525,486-24,368,138 reads (Table 1). At least 89% of the total reads were mapped to the reference genome of L. acetotolerans strain BM-LA14527 (GenBank accession number: LTDX00000000). More than 43% of the total reads were mapped in proper pairs and the duplication rates range from 34% to 52%. Gene expression levels determined by the average values of Reads Per Kilobases per Million reads (RPKM) demonstrated that approximately 1800 total predicted genes were expressed in the normal, mid-term and VPNC states of the 3 L. acetotolerans strains (Table S1). Adjusted P value < 0.05 and │ log 2 (fold change)│ > 1 were used to identify DEGs. Diversity in identified DEGs was obtained, including 66/56 (normal state versus mid-term state for BM-LA14526), 35

GO functional analysis of DEGs.
To understand the function of DEGs, GO function annotation and enrichment analysis were performed on the total, up-regulated, and down-regulated DEGs of three strains (Tables  S4-S6, with significantly enriched terms marked in bold). Significantly enriched GO terms composed of three specific categories (biological process, cellular component, and molecular function) of three strains were illustrated in Fig. 3.
As to the enriched GO terms in L. acetotolrans strain BM-LA14526, "organic phosphonate transport", "outer membrane-bounded periplasmic space", and "organic phosphonate transmembrane transporter activity" were significantly enriched in the categories of biological process (BP), cellular component (CC), and molecular function (MF), respectively. For the enrichment GO terms of up-regulated DEGs, "carbohydrate metabolic process" and "1-phosphofructokinase activity" were of significant enrichment in BP and MF categories, respectively. Of the 2 significantly enriched GO terms of down-regulated DEGs, "organic phosphonate transport" and "organic phosphonate transmembrane transporter activity" were identified in BP and MF categories, respectively.
The L. acetotolrans strain BM-LA14528 had the most significantly enriched GO terms in common. Among the 8 significantly enriched terms in BP category, "regulation of transcription, DNA-dependent", "response to stress", "protein folding", and "copper ion export" were the dominant terms. Three significantly enriched terms in CC category were "external encapsulating structure", "polyphosphate kinase complex", and "tricarboxylic acid cycle enzyme complex". The GO term "copper-exporting ATPase activity" highly represented the 11 terms of significant enrichment in MF category. Of the 17 significantly enriched GO terms of the up-regulated DEGs, the dominant terms in BP and MF categories were the same with that of the total DEGs. As for the significantly enriched GO terms of the down-regulated DEGs, 4, 3, and 4 terms were identified in the categories of BP, CC, and MF, respectively, and the dominant term in CC category was the same with that of the total DEGs.  (Tables S7-S9, with significantly enriched terms marked in bold) and significantly enriched COG categories of the three strains were also analyzed ( Fig. 4).
For the COG categories acquired by L. acetotolerans strain BM-LA14526, "[P] Inorganic ion transport and metabolism" (also significantly enriched in the COG categories of down-regulated DEGs) and "[V] Defense mechanisms" were significantly enriched. However, no COG category of the up-regulated DEGs was found to be of significant enrichment. Among the three strains, L. acetotolerans strain BM-LA14527 had the most significantly enriched COG categories, with 3 ("[L] Replication, recombination and repair", "[I] Lipid transport and metabolism", and "[J] Translation, ribosomal structure and biogenesis"), 1 ("[I] Lipid transport and metabolism") and 1 ("[E] Amino acid transport and metabolism") of the total, up-regulated and down-regulated DEGs, respectively. In the COG categories enriched by L. acetotolrans strain BM-LA14528, only two categories (" [O] Posttranslational modification, protein turnover, chaperones" and "[KT] COG1983 Putative stress-responsive transcriptional regulator") of the up-regulated DEGs were determined significantly enriched. KEGG pathway analysis of DEGs. To identify the pathways of DEGs, DEGs were mapped to the KEGG database and KEGG enrichment analysis were performed on the total, up-regulated, and down-regulated DEGs of three strains (Fig. 5, Tables S10-S12, with significantly enriched terms marked in bold).
Of total and up-regulated DEGs in L. acetotolerans strain BM-LA14526, "ABC transporters" was significantly enriched. However, no significantly enriched pathways were in common with the down-regulated DEGs. As for L. acetotolerans strain BM-LA14527, 2 ("Fatty acid biosynthesis" and "Synthesis and degradation of ketone bodies"), 4 ("Fructose and mannose metabolism", "Synthesis and degradation of ketone bodies", "Valine, leucine and isoleucine degradation", and "Glycerolipid metabolism"), and 1 ("ABC transporters") were significantly enriched in the total, up-regulated, and down-regulated DEGs, respectively. The pathways of "Biotin metabolism", "Pathways in cancer", "Protein processing in endoplasmic reticulum", and "Renal cell carcinoma" were identified to be significant enrichment in L. acetotolerans strain BM-LA14528. In addition, "Biotin metabolism" and "Renal cell carcinoma" were also significantly enriched pathways of down-regulated DEGs, while "Protein processing in endoplasmic reticulum" was also significantly enriched in pathways of up-regulated DEGs.
qRT-PCR validation. Ten genes with different expression profiles were selected to validate the results of RNA-seq. The representative 10 genes were selected based on three criteria: (i) Gene function. (ii) Expression level. (iii) Gene position in the genome. The 10 genes (5 up-regulated and 5 down-regulated) were selected from different significantly enriched GO terms, COG categories and KEGG pathways, and located in different positions in the genome. According to the results, the mRNA levels of the 10 representative genes determined by qRT-PCR were consistent with those from the RNA-Seq analysis.

Discussion
In nature, bacteria are unavoidably exposed to variety of stress conditions, thus strategies have evolved to enhance survival under a broad array of conditions. Despite the prevalence of the VPNC state as one of the survival strategies adopted by non-sporulating bacteria, few genetic mechanisms involved in the VPNC state have been elucidated. Thorough understanding of the global gene expression profiling of the VPNC state may provide valuable Our previous study demonstrated that L. acetotolerans strains are hard to culture in routine medium, capable of entering into and resuscitating from the VPNC state, and maintain beer spoilage capacity during VPNC conditions 6 . In this study, RNA-seq and bioinformatics analysis were performed to map global gene expression associated with VPNC L. acetotolerans strains, with the goal of identifying the genetic mechanisms contributing to the VPNC state. After sequencing and mapping to the reference genome, gene expression levels were measured by RPKM and DEGs identified by the adjusted P value < 0.05 and │ log 2 (fold change)│ > 1. According to the expression variation in different phase, DEGs identified in the three L. acetotolerans strains were divided into 3 types, designated type I (the initial phase, from normal to mid-term state), II (the latter phase, from mid-term to VPNC state), and III (the whole process of the VPNC formation) ( Table 2). Among the type I DEGs, 9/12, 34/35, and 43/39 were up/down regulated in L. acetotolerans strains BM-LA14526, BM-LA14527, and BM-LA14528, respectively. Also, 12/30, 89/46, and 0/7 DEGs were identified to be up/down regulated among type II DEGs. As for the 11 type III DEGs of strain BM-LA14526, 7 were up-regulated in the initial phase and down-regulated in the latter phase, 3 were down-regulated in the initial phase and up-regulated in the latter phase, and 1 was down-regulated in the whole process. For type III DEG of strain BM-LA14527, 4/1 were up/down regulated in the whole process, 1 was up-regulated in the initial phase and down-regulated in the latter phase, and 2 were down-regulated in the initial phase and up-regulated in the latter phase. As none of common DEG was found in all of the three strains, the genetic mechanisms of the VPNC entry might vary widely from strain to strain. Notably in strain BM-LA14528, 82 type I DEGs were identified, with only one type II DEG and no type III DEG obtained, which suggested the initial phase of VPNC state formation may play more important role than the latter phase during the complex process of VPNC state formation. Since the horA gene expression levels were diminishing from normal state via mid-term state to the VPNC state, the growth ability of the three L. acetotolerans strains in beer might decrease during the formation of the VPNC state. However, the horA gene was not found in DEGs from all comparison groups, indicating insignificant effect of VPNC state formation on the growth of L. acetotolerans strains in beer.
The functions and pathways of DEGs associated with the formation of the three VPNC L. acetotolerans strains were compared by enrichment analysis. Analysis of GO functional enrichment demonstrated that significantly enriched terms of up-regulated DEGs in L. acetotolerans strains were related to "carbohydrate metabolic process", "pyruvate metabolism", "1-phosphofructokinase activity", "metabolic process", "coenzyme A metabolic process", "hydroxymethylglutaryl-CoA reductase (NADPH) activity", "hydroxymethylglutaryl-CoA synthase activity", "coenzyme binding", "regulation of transcription, DNA-dependent", "response to stress", "protein folding", "copper ion export", and "copper-exporting ATPase activity". It was reported that lactobacilli modify the transport and metabolism of carbohydrates to adapt to the carbon source available in different media and to withstand a panel of stresses [22][23][24][25] . During the formation of the VPNC state, L. acetotolerans strains may up-regulate the carbohydrate metabolic process to withstand stress or adapt to depletion of available nutrients in the culture medium. Activity of NADPH-generating systems was proposed to be critical for maintenance of the viability of VPNC cells 26 . Su et al. induced the VPNC state of Rhodococcus and analyzed the transcriptome of VPNC state determining that "coenzyme binding" was the dominant up-regulated GO term in the molecular function category 27,28 . Coupled with our results, these findings indicate that "coenzyme binding" may somehow play a role in the formation of the VPNC state in more than one species. Environmental stress conditions affect the abundance of proteins involved in transcriptional regulation 29 . Regulation of bacterial signaling and gene transcription is often mediated by two-component regulatory systems, consisting of a membrane-bound histidine protein kinase and a cytoplasmic response regulator 30 . Several up-regulated genes in the category of "stress resistance", "regulation of transcription, DNA-dependent" and "response to stress" potentially indicate involvement of their protein products in the VPNC state 26 .
As for COG category enrichment analysis, significantly enriched up-regulated categories were associated with "[V] Defense mechanisms", "[I] Lipid transport and metabolism", "[O] Posttranslational modification, protein turnover, chaperones", and "[KT] COG1983 Putative stress-responsive transcriptional regulator". "[P] Inorganic ion transport and metabolism" and "[E] Amino acid transport and metabolism" were significantly enriched down-regulated categories. The control of amino acid transport and metabolism might be necessary for maintaining the carbon-nitrogen balance in VPNC cells 26 , with amino acid metabolism-related genes regulated in response to environmental signals, such as an abundant nutrient supply or nutrient deprivation 30 . The up-regulated categories enriched in "defense mechanisms" and "stress-responsive transcriptional regulator" suggests that L. acetotolerans engages defense and stress response systems to finally enter into the VPNC state in order to survive under stress conditions, including brewed beverages. On the other hand, reducing amino acid transport and metabolism may be a strategy to decrease energy requirements and enhance the function of survival systems.
The enrichment analysis of KEGG pathway indicated that, "ABC transporters", "Pyruvate metabolism", "Synthesis and degradation of ketone bodies", "Valine, leucine and isoleucine degradation", "Fructose and mannose metabolism", "Protein processing in endoplasmic reticulum" were dominant up-regulated pathways. The "ABC transporters" pathway is important for the import of essential nutrients and export of toxic molecules in bacteria 31 . These transporters may provide prerequisite nutrients and metabolites for VPNC cell formation or continually rid cells of toxins in adverse envrionments 26 . For lactobacilli, modification of pyruvate metabolism via enzymatic increase had been verified for bacterial survival under acid stress or starvation conditions [32][33][34][35][36][37] . In the current study, the up-regulation of DEGs involved in "pyruvate metabolism" pathway and GO term suggested that L. acetotolerans up-regulated related enzymes to survive in beer (acid and starved environment). The two up-regulated DEGs designated gene 0267 and gene 0269, which were both involved in "Synthesis and degradation of ketone bodies" and "Valine, leucine and isoleucine degradation" pathways of L. acetotolerans strain BM-LA14527, encoded hydroxymethylglutaryl-CoA synthase and acetyl-CoA acetyltransferase, respectively. It was reported that L. buchneri increased the amount of enzymes related to the metabolism of arginine, proline, glycine, serine, and threonine under ethanol stress conditions 38 , while L. brevis increased the amount of cysteine sulfinate desulfinase/cysteine desulfurase-related enzyme and pyridoxal phosphate-dependent decarboxylase under hop-stressed conditions 39 . Similar to L. buchneri and L. brevis, L. acetotolerans might increase these genes in response to alcohol or hops. Formation of the VPNC state was also affected by other up-regulated genes with unknown functions. A total of 5 up-regulated DEGs associated with the "Fructose and mannose metabolism" pathway in L. acetotolerans strain BM-LA14527 encoded fructokinase, tagatose-6-phosphate kinase, PTS fructose transporter subunit IIC, mannitol-1-phosphate 5-dehydrogenase, and glycosyl transferase. The single up-regulated DEG (gene 1112) involved in "Protein processing in endoplasmic reticulum" pathway L. acetotolerans strain BM-LA14528 encoded molecular chaperone.
Dominant down-regulated pathways in L. acetotolerans strains included "Fatty acid biosynthesis", "Biotin metabolism", "Cyanoamino acid metabolism". The presence of such compounds increase stress on lactobacilli causing decreased synthesis of enzymes involved in fatty acid biosynthesis 36,[40][41][42][43][44] . However, increases in fatty acid biosynthetic enzymes in lactobacilli have been previously shown to protect the cellular membrane against environmental stresses 45 . The down-regulation of genes involved in "Fatty acid biosynthesis" pathway in the three L. acetotolerans strains might be a stress response and survival part of the genetic mechanism in the formation of the VPNC state. Furthermore, two other down-regulated genes encoded biotin-(acetyl-CoA-carboxylase) ligase (involved in "Biotin metabolism" pathway) and asparagine synthase (involved in "Cyanoamino acid metabolism" pathway), might somehow affect the formation of VPNC L. acetotolerans strain BM-LA14528.
It is noteworthy that the significantly enriched GO terms, COG categories, and KEGG pathways identified in the VPNC state between the three L. acetotlerans used in this study differed considerably from one another and findings that were previously published using other bacterial species 26,27,38 . Thus, the genetic mechanisms contributing to VPNC entry and maintenance are still loosely defined and will require validation by gene knockout studies to confirm their role in this biological process.
This study clarified the genome-wide transcriptional variation in the formation of the VPNC L. acetotolrans strains. Gene expression of the VPNC state indicated strong involvement of stress response pathways and reductions in key transport, metabolic, and enzymatic processes may play an important role in its formation. Moreover, these results determined that the molecular response underlying VPNC entry is complicated and differs between strains and species. It is unclear whether these responses are conserved in all beer types or all conditions that promote VPNC formation but are important questions for further investigation. Additionally, further studies Scientific RepoRts | 6:36753 | DOI: 10.1038/srep36753 should be conducted to determine whether other beer-spoilage bacteria can enter into the VPNC state and if they express the similar transcriptomic changes as L. acetotolerans. Elucidating how beer-spoilage bacteria regulate VPNC entry will afford opportunities to inhibit this process and alleviate beverage contamination and food safety concerns in the brewing industry.

Methods
Bacterial strains and growth conditions. Three L. acetotolerans strains, BM-LA14526, BM-LA14527, and BM-LA14528, were isolated from turbid finished beers (identified to be spoilage) in Guangzhou (China) in 2014 and grown anaerobically at 26 °C in MRS broth (Oxoid, UK).
Induction of the VPNC state. Induction of the VPNC state in strains BM-LA14526, BM-LA14527, and BM-LA14528 was conducted by beer subculture 7 . Firstly, the assay of prolonged adaptation to beer described by Suzuki et al. was carried out 8 . Since L. acetotolerans strain was found to be unable to grow in high-bitterness beers, Chinese beers were used for induction of the VPNC state. Approximately 10 7 cells of L. acetotolerans BM-LA14526, BM-LA14527, and BM-LA14528 were inoculated and anaerobically subcultured at 26 °C in 10 mL aliquots of the degassed commercial beer. The interval of each subculture was 7 days. Afterwards, the bacterial strain was induced to entry into VPNC state by continuous passage in beers. 1 mL exponentially growing bacterial strain which was considered to be the 0th generation was inoculated to 10 mL degassed beer and anaerobically cultured at 26 °C for 7 days to be the 1st generation. Then 1 mL 1st generation culture was inoculated to 10 mL degassed beer and anaerobically cultured at 26 °C for 7 days to be the 2 nd generation, the amount of 2 nd generation culture were measured by plate colony counting method. All the studies related to entry into VPNC state were performed in independent triplicate biological experiments. The number of culturable cells, total cells, and cell viability were accessed by a conventional plate culture procedure, AODC method, and Live/Dead BacLight bacterial viability kit (Molecular Probes, USA), respectively 7 .
RNA extraction, library construction and sequencing. Total RNA of L. acetotolerans samples were extracted using the Bacterial Total RNA Extraction kit (Sigma-Aldrich, USA) according to the manufacturer's instruction. The quality of total RNA was determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). RNase-free DNase I (Ambion Inc., Austin, TX, USA) and MICROBExpress TM kit (Ambion Inc., Austin, TX, USA) were used to remove trace DNA and ribosomal RNAs, respectively. The mRNA was fragmented ultrasonically and then reverse transcribed to cDNA. RNA-seq libraries were constructed using the Illumina Paired End Sample Prep kit and all samples were sequenced using the Illumina Hiseq 2500 platform.
RNA-seq data analysis. Raw reads were generated from image data and stored in FASTQ format. All the raw data have been deposited in the NCBI Short Read Archive (SRA) database under accession no. SRX1687039. Clean reads were aligned to L. acetotolerans BM-LA14527 (GenBank accession no. LTDX00000000) using TopHat 46 . Only mismatches and read gaps with no more than 2 bases were allowed in the alignment. Gene coverage was calculated by the percentage of genes covered by reads. The nine samples were combined to 9 comparison groups (BM-LA14526 normal state versus BM-LA14526 mid-term state, BM-LA14526 normal state versus BM-LA14526 VPNC state, BM-LA14526 mid-term state versus BM-LA14526 VPNC state, BM-LA14527 normal state versus BM-LA14527 mid-term state, BM-LA14527 normal state versus BM-LA14527 VPNC state, BM-LA14527 mid-term state versus BM-LA14527 VPNC state, BM-LA14528 normal state versus BM-LA14528 mid-term state, BM-LA14528 normal state versus BM-LA14528 VPNC state, and BM-LA14528 mid-term state versus BM-LA14528 VPNC state) to determine differentially expressed genes (DEGs) under each culture condition. DEGs were identified using DEGseq panormalage 47 and the method described by Audic S. 48 . The gene expression level was measured using RPKM method 49 and P values adjusted using the edgeR panormalage 50 . Genes with an adjusted P value < 0.05 and │ log 2 (fold change)│ > 1 were identified as DEGs. Gene Ontology (GO) 51 , Clusters of Orthologous Groups (COG) 52 and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway 53 analysis were performed to determine the putative function of DEGs.
Quantitative real-time RT-PCR (qRT-PCR) validation. In order to validate the RNA-seq data, qRT-PCR was performed to quantify the mRNA transcripts of 10 selected DEGs using the CFX96 qRT-PCR System (Bio-Rad, USA) according to the manufacturer's instructions. The RNA samples were extracted using Bacterial Total RNA Extraction kit (Sigma-Aldrich, USA) and treated with RNase-free DNase I to remove DNA contamination. Each qRT-PCR reaction was conducted in a final volume of 25 μ L. The thermal cycling profile was as follows: 42 °C for 60 min and 72 °C for 10 min; 45 cycles of 95 °C for 10 s, 60 °C for 30 s, 70 °C for 1 min and a final extension of 68 °C for 7 min. Negative control samples containing sterile water were also included. The cycle threshold values (C T ) were determined and the relative fold differences were calculated by the 2 −ΔΔC T method 54 using 16S rRNA as the reference gene. Three independent experiments were run in triplicate.