A number of parameters potentially affecting microbial population structure during retting were examined. It is commonly admitted by farmers that the maturity of flax plants has a direct impact on the retting time and influences the choice for the pulling (up-rooting) date. Generally, straw from younger plants (flowering/green capsule stage) rets more quickly than that of more mature plants (yellow/brown capsule stage). This is thought to be related to differences in cell wall composition (e.g., pectin/lignin modifications and/or deposition) and water content (Meijer et al., 1995; Day et al., 2005; Akin, 2013). Our results showing that there was no significant difference in microbial communities and colonization dynamics between the early vs. standard cultures would suggest that differences in retting time may indeed be related to differences in cell wall structure and not to population differences.

Compared to litter decay that normally proceeds undisturbed, dew-retting is a semi-controlled process during which the straw swaths are turned by farmers to obtain a more uniform fiber separation. Our analyses revealed that this practice had a significant effect on both bacterial and fungal community membership and structure of the flax straw microbiome confirming a real microbiological effect of swath turning that probably contributes to a more uniform retting.

Although our results indicated no significant correlation between measured climatic conditions (temperature and rainfall) and community structures during the retting period it is important to remember that our study was conducted within a single year. It is possible that significant variations in community structures may occur between different seasons and further work is necessary to clarify this point.

### Microbial Dynamics

During dew-retting the relative abundance of the Bacteroidetes phylum increases while that of the Protobacteria decreases. A similar dynamic also occurs during biodegradation of field biomass from different angiosperm species (e.g., Arundo donax, Eucalyptus camaldulensis, and Populus nigra) suggesting, as might be expected, that similarities exist between the temporary dew-retting ecosystem and degradation of lignocellulose in the field (Ventorino et al., 2015). Interestingly, the bacterial dynamics of flax dew-retting appear to be closer to that of field lignocellulose degradation than to that observed during flax water retting where Protobacteria increased during retting (Zhao et al., 2016). In this latter case, the phylum Proteobacteria was mainly represented by the genera Azotobacter that increased during retting and (to a much lesser extent) Enterobacter that remained constant. For fungal phyla we observed an increase in the relative abundance of Ascomycota at the expense of Basidiomycota in contrast to the situation generally observed during both field and forest litter decomposition (Schneider et al., 2012; Kuramae et al., 2013; Voˇríšková and Baldrian, 2013). The observed increase of Ascomycota was due to the saprophytic Altenaria species (Dang et al., 2015) that has previously been linked to later stages of dew-retting (Brown et al., 1986). In contrast, Altenaria species are more abundant during initial stages of litter decay (Snajdr et al., 2011). Our results also indicated that C. herbarum and Epicoccum nigrum contributed to the increase in Ascomycota relative abundance. During this stage less recalcitrant components of the biomass (pectins, and hemicelluloses) are progressively degraded (Dilly et al., 2001). Contrary to litter decay, dew-retting is a semi-controlled process and the challenge is to limit degradation of major quality related polymers such as crystalline cellulose. In this context, changes in the relative abundance of Ascomycota vs. Basidiomycota could represent an interesting bioindicator of retting progress.

More detailed information on population dynamics at different time points during retting was provided by analyzing the relative abundance of OTUs at different taxonomic rank (e.g., phyla, classes, or genus/species level). The most abundant bacterial OTU corresponded to Sphingomonas sp. that was present throughout most of the retting period in both early and standard cultures. Although Sphingomonas species have been previously identified during bamboo and hemp retting, as well as in forest litter microbiome, this is the first time they have been found in flax retting (Fu et al., 2011; Ribeiro et al., 2015) (Urbanová et al., 2015). These species are able to hydrolyze terminal non-reducing alpha-L-rhamnose residues in alpha-L-rhamnosides giving them the ability to degrade pectin (rhamnogalacturonan I and rhamnogalacturonan II) in the middle lamella (Hashimoto and Murata, 1998). Another Sphingomonas species, S. paucimobilis is also able to degrade lignin (Masai et al., 1999; de Gonzalo et al., 2016). The second most abundant OTU corresponded to P. rhizosphaerae, present during the early and medium retting stages but decreasing in latter stages. A number of Pseudomonas species have previously been associated with retting of different fiber plants (e.g., flax, hemp, jute, ramie) (Rosemberg, 1965; Munshi and Chattoo, 2008; Duan et al., 2012; Ribeiro et al., 2015). Pseudomonas sp. is considered as one of the most efficient lignin degradation bacterium (Shui Yang et al., 2007) and the genomes of both Pseudomonas putida and Pseudomonas aeruginosa contain genes encoding endoglucanases (Talia et al., 2012). Other abundant OTUs corresponded to Rhizobium, Pedobacter, and Flavobacterium that are known to show pectinase, cellulose, and hemicellulose activities (Mateos et al., 1992; McBride et al., 2009; López-Mondéjar et al., 2016). In addition, Pedobacter has also been identified during bamboo and hemp retting (Fu et al., 2011; Ribeiro et al., 2015) or forest litter degradation (Urbanová et al., 2015). In contrast to Sphingomonas and Pseudomonas, these organisms become more abundant toward the end of the retting period and could be associated with "over-retting" when the structural integrity of the fiber starts to be degraded.

In contrast to the more evenly distributed abundance of the bacterial OTUs, fungal OTUs were dominated by one major species—C. herbarum—that rapidly increased during early retting. This species, as well as the third most abundant OTU (E. nigrum) are known to be common dew-retting agents and are believed to degrade cellulose (Brown, 1984). Of the other fungal OTUs, all have previously been associated with dew-/water-retting except for Itersonilia perplexans. Interestingly, our results also indicated that Alternaria alternata is present at the start of retting. Traditionally, the appearance of this species is used as a signal that retting is starting to go too far and that the swaths should be collected (Brown et al., 1986).

### Hydrolytic Enzyme Potential

Prediction of hydrolytic enzymes potentially present during retting was performed by using PICRUSt (Langille et al., 2013). This software successfully predicts bacterial enzymatic activities represented in different databases (e.g., KEGG Ortholog, COGs, or CAZy). Overall, a large collection of enzyme activities targeting both the main backbones and side chains of the major polysaccharide polymers were identified. Based on OTU counts, ∼38, 43, and 19 percent of the total hydrolytic enzyme potential targeted pectins, hemicelluloses, and cellulose, respectively. Despite the clear dynamics and significant changes in the straw microbiome these values remained constant throughout the retting period. Similar software does not exist for predicting fungal enzyme potential. This represents an important hurdle for obtaining a complete overview of the dew-retting process as fungi are major producers of extracellular hydrolytic enzymes (Schneider et al., 2012). Nevertheless, FUNGuild analysis showed that pathogenic taxa, present at the beginning of retting are progressively replaced by saprophytic fungi, more able to degrade lignocellulose. This change is most likely related to the fact that flax plants are still living when up-rooted.

In conclusion, we have shown that HTS metabarcoding is a powerful technique for analyzing complex bacterial and fungal community dynamics during flax dew-retting that can be used to identify different factors affecting the microbiota and—potentially—fiber isolation and quality. However, these results were obtained on samples retted in 1 year and it will be necessary to validate these data over several seasons. The use of PICRUSt data allows a predictive study of potential bacterial hydrolytic activity but should be coupled in future studies with alternative meta-omics methods such as metatranscriptomic or metaproteomic coupled with metagenomics to facilitate the assembling with appropriate reference genomes (Schneider et al., 2012; Dai et al., 2015; Hesse et al., 2015; Kuske et al., 2015; Wu et al., 2015). Such an approach would not only allow confirmation of bacterial enzyme dynamics but would also enable identification of fungal enzymes involved in this process.

### AVAILABILITY OF DATA AND MATERIALS

The microbial DNA sequencing data sets supporting the results in this article are available at the EBI ENA with accession number PRJEB20299.

### AUTHOR CONTRIBUTIONS

Conceptualization: CD, SG, and SH; Methodology: CD, SG, and SH; Experimentation: CD; Bioinformatic and statiscal analysis: CD; Writing—original draft: CD; Writing—Review and editing: SG and SH; Funding Acquisition: SG and SH.

### FUNDING

This work was funded within the framework of the collaborative French "Future project" SINFONI. CD thanks the region of Hauts-de-France and Bpifrance for their financial support.

### ACKNOWLEDGMENTS

The authors would like to thank the following people/organizations: The flax farmer C.A.L.I.R.A (Coopérative Agricole LInière de la Région d'Abbeville) for growing and retting flax plants, in particular Vincent DELAPORTE for his valuable assistance. The Genomic and Transcriptomic platform of Genopole Occitanie-Toulouse (INRA, GeT-Plage, http://get.genotoul.fr/) where sequencing was performed and more particularly Catherine ZHANCHETTA and Olivier BOUCHEZ. The Plateau d'Ecologie moléculaire et biochimie Evo-Eco-Paléo (CNRS—UMR 8198, Evolution, Ecologie et Paléontologie) for use of the LightCycler 480 and epMotion and more particularly Cécile GODÉ and Anne-Catherine HOLL for their technical support. Antoine PORTELETTE, Julien LE

### REFERENCES

ROY, Sandrine ARRIBAT, and Brigitte CHABBERT for their valuable assistance during the sampling. Alexandrine THORE for her invaluable help in the Supplementary Figure 2 graphical design.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2017.02052/full#supplementary-material

Supplementary Figure 1 | Rainfall, daily amplitudes of temperature and moisture, and retting time points at Martainneville during the 2014 dew-retting period. Raw data from the Abbeville station available on the infoclimat website (https://www.infoclimat.fr).

Supplementary Figure 2 | Schematic representation of the experimental site at Martainneville during the 2014 dew-retting campaign.

Supplementary Figure 3 | Rarefaction curves of observed OTUs in bacterial (16S) soil (A) and stem (B) samples, defined at a 97% sequence similarity cut-off.

Supplementary Figure 4 | Rarefaction curves of observed OTUs in fungal (ITS) soil (A) and plant (B) samples, defined at a 97% sequence similarity cut-off.

Supplementary Figure 5 | Alpha diversity of bacterial (A–C) and fungal (D–F) communities in soil and plant samples during flax dew-retting. (A,D) Richness (Chao1), (B,E) Evenness (Heip's), and (C,F) Diversity (Inverse Simpson). Different colors indicate the source samples (light brown, soil early harvest; dark brown, soil standard harvest; light green, plant early harvest; dark green, plant standard harvest). Plots show the distribution of results of 20,548 bacterial and 42,436 fungal sequences subsampled from each sample 1,000 times and calculated for the average (error bars for Chao1 and Inverse Simpson correspond to lower and higher bound 95% confidence intervals and for Heip's represent standard deviation).

Supplementary Figure 6 | Before- (A,C) and after- (B,D) swath turning biomarkers for bacteria (A,B) and fungi (C,D). Biomarkers were identified by analyzing differential OTUs abundance using Metastats, Indicator and LEfSe.

Supplementary Figure 7 | Heatmaps showing the total counts of all CAZyme families in early and standard harvests obtained from the OTUs table (Greengenes Database used for consensus taxonomy) and generated by PICRUSt software.

Campilho, R. D. S. G. (2015). Natural Fiber Composites. Boca Raton, FL: CRC Press.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Djemiel, Grec and Hawkins. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Sequencing of the Cheese Microbiome and Its Relevance to Industry

Bhagya. R. Yeluri Jonnala1,2, Paul L. H. McSweeney 1 , Jeremiah J. Sheehan2 and Paul D. Cotter 2,3 \*

1 Food and Nutrition Deptartment, University College Cork, Cork, Ireland, 2 Teagasc Food Research Centre, Fermoy, Ireland, 3 APC Microbiome Ireland, Cork, Ireland

The microbiota of cheese plays a key role in determining its organoleptic and other physico-chemical properties. It is essential to understand the various contributions, positive or negative, of these microbial components in order to promote the growth of desirable taxa and, thus, characteristics. The recent application of high throughput DNA sequencing (HTS) facilitates an even more accurate identification of these microbes, and their functional properties, and has the potential to reveal those microbes, and associated pathways, responsible for favorable or unfavorable characteristics. This technology also facilitates a detailed analysis of the composition and functional potential of the microbiota of milk, curd, whey, mixed starters, processing environments, and how these contribute to the final cheese microbiota, and associated characteristics. Ultimately, this information can be harnessed by producers to optimize the quality, safety, and commercial value of their products. In this review we highlight a number of key studies in which HTS was employed to study the cheese microbiota, and pay particular attention to those of greatest relevance to industry.

Keywords: high throughput sequencing, cheese, microbiota, sensory characteristics, metatranscriptomics, industry

### INTRODUCTION

Cheese has a diverse microbial community, which indeed can vary within the cheese from the core to the surface that is greatly influenced by manufacturing including ripening conditions. Understanding the composition of this community (microbiota), and its impact on the quality and safety of cheese products, is of critical importance. In addition to, in the majority of cases, consciously added starter and adjunct bacteria (which are added as a supplement), cheese contains a heterogeneous variety of other, non-starter, microorganisms. These various microorganisms can play vital roles in the development of the organoleptic properties of cheese (Fox et al., 2000), nutrient composition, shelf-life, and safety.

Historically, culture-based microbiology techniques were used to gain an understanding of the microbial component of cheese. However, it has become increasingly clear that this approach can be limited in its ability to detect "difficult-to-culture" or sub-dominant microorganisms, thereby potentially providing misleading results. As a result, culture-independent approaches have become

#### Edited by:

Florence Abram, National University of Ireland Galway, Ireland

#### Reviewed by:

Baltasar Mayo, Consejo Superior de Investigaciones Científicas (CSIC), Spain Francesco Grieco, Istituto di Scienze delle Produzioni Alimentari (ISPA), Italy Rodrigo Bibiloni, Arla Foods, Denmark Fernanda Mozzi, CERELA-CONICET, Argentina

> \*Correspondence: Paul D. Cotter<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Food Microbiology, a section of the journal Frontiers in Microbiology

Received: 19 December 2017 Accepted: 30 April 2018 Published: 23 May 2018

#### Citation:

Yeluri Jonnala BR, McSweeney PLH, Sheehan JJ and Cotter PD (2018) Sequencing of the Cheese Microbiome and Its Relevance to Industry. Front. Microbiol. 9:1020. doi: 10.3389/fmicb.2018.01020

**80**

**Abbreviations:** HTS, High Throughput Sequencing; NGS, Next Generation Sequencing; PDO, Protected designation of Origin; CIP, Cleaning-in-Place; NMC, Natural Milk Starter Cultures; NWC, Natural Whey Cultures; VOC, Volatile Components.

increasingly popular. These approaches involve DNA, and occasionally RNA, based molecular methods, with high throughput sequencing (HTS) gaining particular attention by virtue of its potential to provide an insight into the total microbial community of the cheese. The application of HTS can provide valuable information relating to the influence of geography, manufacturing processes, climatic conditions, seasonal variation, use of raw or pasteurized milk, and a variety of other factors on the cheese microbiota (**Figure 1**). The application of HTS for microbiota analysis can involve any of three primary approaches i.e., (1) Amplicon sequencing—whereby a fragment of highly conserved gene, ideally with variable regions therein, is used for sequencing, with comparison to databases allowing taxonomic assignment. The 16S rRNA gene is most frequently used and provides an insight into the bacterial composition of samples (generally to the genus level), (2) Shotgun metagenomic sequencing—which involves non-targeted sequencing of the DNA in a sample and, again on the basis of comparison with databases, can be used to classify all of the microorganisms, i.e., not just bacteria present (to species, or even strain, level) as well as information regarding the functional potential of the community, (3) Metatranscriptomics (RNASeq)—whereby total mRNA in the sample is sequenced (after first being converted to cDNA) to reveal the extent to which different genes are expressed and, in turn, the relative activity of different components of the community. These three experimental procedures differ with respect to library preparation, sequencing strategy, and the approaches taken for bioinformatic analysis. Some of the most common outputs from these analyses relate to α-diversity, β-diversity, and a determination of the relative abundance of different taxa. α-diversity is a measurement of diversity, such as total number of species, within a sample and measures such as the Simpson and Shannon diversity indices are employed. β-diversity reflects differences in diversity across different samples, with Bray-Curtis and UniFrac being among the distance metrics used. For metagenomics and metatranscriptomics, several tools are available for key steps, i.e., binning, assembly (where relevant), and mapping/assignment of sequences obtained (Di Bella et al., 2013; Walsh et al., 2017). There have been a number of research and review papers dedicated to describing methodologies used for sequencing (Di Bella et al., 2013) work flow, limitations, applications to dairy foods (Ercolini, 2013; Walsh et al., 2017), comparing sequencing platforms, and bioinformatics pipelines as well as highlighting pros and cons associated with all of the above (Kelleher et al., 2015; Clooney et al., 2016; Walsh et al., 2018) Depending on sequencing depth, HTS can facilitate the detection of taxa present at low levels, a feature of great value in the context of testing for spoilage or pathogenic microbes. HTS technology can also be employed to study the microbial diversity of protected designation of origin (PDO) cheeses to establish which components are responsible for the authentic taste, flavors, and textures of these products (Alegría et al., 2012; Fuka et al., 2013; De Filippis et al., 2014; Delcenserie et al., 2014; De Pasquale et al., 2014, 2016; Riquelme et al., 2015; Dalmasso et al., 2016; Parente et al., 2016; Giello et al., 2017; Gonçalves Dos Santos et al., 2017; Li et al., 2017). The same technology could, in the future, be employed to establish authenticity. However,

as these topics, and studies related to microbial dynamics of food ecosystems, processing environments, and production facilities have been reviewed in depth recently (Mayo et al., 2014; Galimberti et al., 2015; Bokulich et al., 2016; De Filippis et al., 2017; Doyle et al., 2017b; Macori and Cotter, 2018), we are not addressing those topics here. In this review, we specifically focus on providing an overview of studies in which HTS has been applied to provide insights into the microbiota of different cheeses that are of value from an industrial perspective as well as the factors that modulate these microbes (**Table 1**) and the identity of adventitious and other microbes (**Table 2**).

### IDENTIFICATION OF FACTORS INFLUENCING THE DEVELOPMENT OF THE CHEESE MICROBIOTA

Microbial composition and diversity differs from raw to pasteurized milk, and between curd, whey, and cheese. The raw milk microbiota is influenced by microbes present in the teat canal, the surface of teat skin, hygiene practices, animal handlers, and the indigenous microbiota of equipment and storage containers. The origin of the milk would also seem TABLE 1 | HTS detection of increase and decrease of microbiota by the action of some influencing factors.

Abundance of bacteria in different cheese samples was detected by HST, ↑ (increase) and ↓ (decrease) indicates levels of certain bacteria as influenced by factors stated

to influence the levels of diversity therein, with cow's milk appearing to be more diverse than that from goats and sheep (Quigley et al., 2012). The type of grazing system employed, i.e., extensive or semi-extensive, affects the dynamics of teat skin microbiota and thus, in turn, the microbes in raw milk cheeses. Indeed, in one study, 27% of bacteria detected in raw milk cheese were also found on the teat surface. These bacteria included species involved in the production of flavor, aroma, and color development, such as Brevibacterium linens, Staphylococcus equorum, and Lactic Acid Bacteria (LAB) such as Lactococcus lactis, Lactococcus chungangensis/raffinolactis, and Lactobacillus casei/paracasei, which can contribute to protein and fat metabolism (Frétin et al., 2018). The sensory attributes of cheese can also be influenced by the type of feed (e.g., hay, silage) given to cows. Indeed, this was recently demonstrated in a study of Caciocavallo cheese where salty, sour, bitter, and umami flavors were lower in "hay"-cheese compared to "silage"-cheese, and higher levels of tenderness and oiliness were also evident

TABLE 2 | Adventitious, previously overlooked, and spoilage bacteria identified by HTS analysis of different cheese matrix and dairy environment.

Corynebacterium casei

Hafnia alvei

#### TABLE 2 | Continued

in "silage"-cheese relative to "hay"-cheese (Giello et al., 2017). High temperature treatment (when pasteurized milk is used) and low pH also contribute to the selection of specific bacteria in some artisanal cheeses that are made in the absence of starter bacteria (De Filippis et al., 2014). These and other factors, such as salt content, degrees of ripening, addition of ingredients like herbs, and spices also influence the cheese microbiota and, in turn, organoleptic properties of cheese. The influence of such factors was demonstrated in one representative study of Irish artisanal cheeses where, for example, proportions of lactobacilli increased and lactococci decreased due to the addition of herbs, while high salt content suppressed the growth of Leuconostoc and Pseudomonas (Quigley et al., 2012).

The microbial profile of cheeses can differ between producers, according to the process of production and the inclusion, or absence, of starters. Plaisentif is a traditional Italian cheese made from fresh full-fat raw cow's milk without an added starter but, rather, milk from the previous evening (kept at <10◦C) and bovine liquid rennet are employed. HTS analyses conducted on samples from nine different producers revealed the presence of the genera Lactococcus, Lactobacillus, and Streptococcus in the core of all samples but also established that these are present in different proportions and are responsible for variations in the niche-specific characteristics of Plaisentif cheese (Dalmasso et al., 2016). Greater homogeneity was apparent in a recent study involving Herve, a PDO cheese from Belgium. In the past, Herve was made from raw milk but, for safety reasons due to the possible presence of Listeria monocytogenes, pasteurized milk has begun to be used. The 16S rDNA analysis showed that 95% of the microbial composition was same in both raw and pasteurized cheese; hence the characters of both forms of the cheese are similar. This might be because of consistency in the manufacturing process across artisanal producers (Delcenserie et al., 2014). The impact of temperature was also evident when the scalding temperature during the manufacture of Dutch-type cheeses was increased from 37 to 39◦C. This change resulted in an increase in Lactobacillus sp. and Leuconostoc sp. during ripening. This change also promoted flavor formation by lysing lactococci, which release enzymes like peptidases (Porcellato and Skeie, 2016). The influence of manufacturing processes on bacterial communities was also noted in the case of Oscypek cheese. In that instance, a decrease in the counts of enterobacteria was apparent from the curd to the smoked cheese as a result of the smoking process, thus improving the quality and safety of the traditional product (Alegría et al., 2012). Cotija is another example of a cheese in which fluctuations in composition have been observed. This handmade Mexican cheese is prepared from raw milk with no added starters. Ripening of Cotija is performed in an open environment and, thus, humidity, rain, and temperature are important parameters to be considered. Furthermore, the wooden surfaces used for kneading the curd while salting, vats, and other tables used in the process of cheese-making all act as a source of the microbes found in this cheese. About 80% of the bacterial population of Cotija was consisted of a combination of Lactobacillus plantarum, Leuconostoc mesenteroides, and Weissella paramesenteroides. These species are thought to be responsible for the development of authentic flavor compounds arising from their lipolytic and proteolytic activity during ripening (Escobar-Zepeda et al., 2016). It is thought that the wooden surfaces are most likely the source of the Leuconostoc and Weissella species present (Settanni et al., 2012). It is important to note, however, that even when cheeses are produced on a larger industrial scale, including the use of starters to minimize variability, differences in the associated microbiota can be seen. This was apparent from a study of different batches of brine-salted continental-type cheeses made during the same day of production. In this study, 16S rRNA sequencing revealed that cheese produced later in the production day had a more diverse bacterial composition than cheese produced earlier, possibly due to the accumulation of bacteria in the system before cleaning-in-place (CIP) processes are employed at the end of the manufacturing day (O'Sullivan et al., 2015a).

Cheese rinds have a very complex microbiota when compared to core samples; this microbiota varies between types of rinds, degrees of ripening, and the environmental conditions. HTS analyzes of 11 artisanal Irish cheeses (soft, hard, and semi hard) revealed the presence of 19 genera, of which Lactococcus, Leuconostoc, and Lactobacillus were in both rind and core samples. Corynebacterium, Facklamia, Flavobacterium, and Cronobacter were detected in rind samples only. The relative proportions of lactococci were found to be high in naturally ripened rinds relative to smear/washed rinds. This is an ultimate result of the washing of rind in that cheese (Quigley et al., 2012). In one seminal study, 137 different rind samples that included bloomy, natural, and washed types from 10 countries were examined using HTS to understand the mechanisms of how microbial communities form multispecies biofilms. It was noted that rind samples from cheeses made at different geographical regions showed similar patterns of microbial clusters, showing that, in this instance, distance did not have impact on microbial community. Some of the specific patterns noted were as follows; the fungus Galactomyces, which is responsible for the formation of a dense white rind on bloomy rind cheeses (e.g., Brie, Camembert), is positively correlated with moisture. Natural rind cheeses (e.g., Clothbound Cheddar, St. Nectaire, and Tomme de Savoie) are abundant in fungi and bacteria, such as Scopulariopsis, Aspergillus, Actinobacteria, and Staphylococcus, which are negatively correlated with moisture. The microbial composition of washed rind cheese (e.g., Gruyere, Epoisses) was found to be a mix between that of bloomy and natural rinds (Wolfe et al., 2014). Finally with regard to cheese rinds, it is also worth noting that HTS is beginning to be combined with other techniques, such as microscopy, to identify the fungal and bacterial interactions, and the bacterial dispersal on fungal networks (Zhang et al., 2018).

As already noted above, many artisanal cheeses are produced without the addition of starter cultures. Undefined natural milk starter cultures (NMC), produced by heat treatment of raw milk (60–63◦C, 20–30 min) followed by incubation at high temperatures (39–42◦C), are used to make some Italian PDO cheeses. Such cultures, prepared by backslopping (i.e., inoculating milk with a previous batch of the culture) using raw milk collected from pasta filata cheese plants, were found to be dominated by Streptococcus thermophilus and Lactobacillus delbrueckii (Parente et al., 2016). Natural whey cultures (NWC), produced by incubating the whey from previous batches of cheese production at high temperature (39–54◦C), is also used to make some PDO Italian cheeses. These NWCs were dominated by L. lactis, Lactobacillus fermentum, S. thermophilus, L. delbrueckii, and Lactobacillus helveticus, though the abundance of these species was found to be significantly different across the three cheese types, Mozzarella, Grana Padano, Parmigiano Reggiano (De Filippis et al., 2014). The combined pressure of temperature and pH selects for these bacteria. In the case of high-moisture Mozzarella cheese, the mode of acidification, i.e., through addition of citric acid, by starters, or by using undefined starters, greatly influences the growth of the dominant microbial components, thereby differentiating cheeses produced by different dairies (Guidone et al., 2016). This study demonstrated that HTS can be used by Mozzarella producers to identify the dominant bacteria present and, in turn, helps them to choose the mode of acidification that is best suited to obtaining specific desirable characteristics. Pico is an artisanal Azorean cheese, made from raw milk, animal rennet, and salt without any starters. HTS analysis of this cheese by Riquelme et al. (2015) established that Pico manufacture is a Lactococcus-driven process with contributions from accompanying Lactobacillus and Gammaproteobacteria populations. The availability of this information is useful in terms of selecting starter/adjunct cultures specific for Pico cheese, as well as highlighting some hygienic and safety measures that need to be taken by industries to improve the shelf life of these cheeses (Riquelme et al., 2015). From a safety perspective, HTS detection of microbes in Danish raw milk found that L. helveticus, S. thermophiles, and L. lactis are dominant (Masoud et al., 2011) and, following spiking with Listeria innocua and Staphylococcus aureus during manufacture of this cheese to detect the fate of pathogenic bacteria during ripening, it was found that these species were absent from the ripened cheese. The control of these microbes was presumed to be due to acidification, environmental conditions or antimicrobials produced by the LAB (Masoud et al., 2012).

Serpa is an artisanal PDO Portuguese cheese, made from raw ewe's milk using aqueous infusion of Cynara cardunculus L. (artichoke), without starters. The fungal communities in this cheese are thought to have a key role in developing the associated organoleptic properties. HTS analysis of this fungal community of both PDO and non-PDO forms of Serpa cheeses, highlighted diversity among the yeasts present in a manner suggested that production practices and the associated environment have a considerable influence on the final yeast population. Through these investigations it was possible to alter ripening conditions to favor the dominance of yeast species that are desirable for optimal Serpa production (Gonçalves Dos Santos et al., 2017). Another ewe's milk cheese produced by traditional techniques without starters, in this case from Croatia, has also been the subject of HTS-based analysis. This revealed the presence of pathogenic bacteria in fresh milk and cheese. However, their abundance was low at the end of ripening, highlighting the importance of ripening for 90 days to avoid problems related to the safety of the cheese (Fuka et al., 2013). This approach highlighted the merits of using HTS to amend production processes to alter quality and safety, particularly in instances where neither a starter nor pasteurized milk are used. An alternative, HTSbased approach to studying the microbiota of cheese has been to use single molecular real-time (SMRT) sequencing, a third generation technique that provides longer DNA reads. This technique has been used to sequence a Kazakhstan cheese made using NWC as starters, revealing that L. lactis, L. helveticus, S. thermophiles, and Lactobacillus bulgaricus were the dominant species in this cheese. Differences between the microbiota of this cheese and artisanal cheeses from Belgium, Italy, and Kalmykia (Li et al., 2017) likely explain the varying characteristics of these cheeses.

The house microbiota that colonize the equipment, brine tanks, vats, wooden surfaces, knives, and other surfaces in production facilities can play a key role in shaping the microbial communities of cheese. This colonization depends on the characteristics of the surfaces, nutrient availability and composition, ability of microbes to form biofilms, ecological factors, and on operators and cleaning processes (Stellato et al., 2015). Studies focusing on environmental microbiota have been conducted on cheese plants in which one company produces two different cheese types (Calasso et al., 2016) and, in another instance, two different facilities that produce the same type of cheese (Bokulich and Mills, 2013). In the former case, it was deduced that S. thermophilus colonized, to different degrees, almost all surfaces in the dairy plant and acted as an indirect source of starter inoculation (Calasso et al., 2016). In the latter study, it was revealed that the type of processing facility and selective forces in the environment, such as the temperature in the cheese plants, influenced the growth of specific house microbiota, which were in turn responsible for the chemosensory properties of the artisanal cheese (Bokulich and Mills, 2013).

One of the factors that can contribute to differences across different production facilities is brining. The microbial diversity and composition of brine depends on the type of cheese made, the specific cheese plant and on salinity concentrations. Adventitious bacteria present in brine, such as Staphylococcus equorum, which can have a strong antibacterial activity against L. monocytogenes on cheese surfaces, can play important roles in the development of the flavor and color properties of smear ripened cheeses. However, contaminated brine can in turn contaminate the cheese core and surface, with Mozzarella cheese being among the cheeses that are susceptible to such spoilage. Indeed, as a consequence, it has been suggested that sanitization of brine should be prioritized and that brine should be checked regularly for the presence of contaminants (Marino et al., 2017). Cleaning is an important factor to eliminate spoilage bacteria present on the surfaces in dairy plant, with the type of surface material having a key influence on bacterial adherence. Surface materials such as plastic, in the case of gaskets used for the molding of some cheese, is more porous and not suitable for thorough cleaning as it is sensitive to hot water, which causes corrosion and increases the possibility of bacterial adherence (Stellato et al., 2015). As a consequence, it is recommended that steel be used where possible to enhance the safety of products. The importance of proper washing during the manufacture of fermented foods was also highlighted by Lee et al. (2017), who reported that pathogens were replaced by harmless bacteria, such as LAB, as a result of washing. It has been speculated that an equilibrium may exist between the dairy food production and processing environments and the resultant products, with microbial transfer occurring in both directions, which can affect processing dynamics and the quality of the final products (Stellato et al., 2015).

### IDENTIFICATION OF PREVIOUSLY OVERLOOKED BACTERIA IN CHEESE

Depending on the depth of sequencing, HTS-based taxonomic analysis of cheese can detect microbes present at very low levels (O'Sullivan et al., 2015b; Cotter and Beresford, 2017), which would previously have been overlooked. While, in many cases, it is not yet clear what roles these microbes play, if any, in the context of cheese microbiology, an initial awareness of their presence can lead to further investigations. Here we are highlighting some such subdominant populations as revealed in a selection of representative studies. In one such study, a subdominant taxon was identified after culturebased enrichment. Enriching a food sample in broth for 24 h followed by selective agar plating is a traditional method to isolate subdominant bacteria, and is an especially relevant technique for the detection of foodborne pathogens, in the sample. However, overnight growth in non-selective broth can result in a bias toward the detection of fast-growing, relative to slow-growing, bacteria. In this instance HTS of overnight enriched Queso Fresco cheese sample did successfully reveal the presence of Exiguobacterium, which at that point had not previously been documented in cheese (Lusk et al., 2012). It should be noted, however, that the majority of HTS studies take place without an enrichment step. In one typical study, it was noted that Plaisentif cheese contains some rare genera, representing 0.01–0.0001% of total reads, that correspond with Flavobacterium, Brevibacterium, Salinicoccus, Vagococcus, Anaerobacillus, Sphingobacterium, and Klebsiella, among others (Dalmasso et al., 2016), with Klebsiella in particular being regarded as indicative of poor hygiene conditions during manufacturing. In contrast, sub-dominant genera found in Cotija cheese include the halophilic microbes Dehalobacter, Desulfohalobium, Halomonas, Thermohalobacter, and Haloquadratum (Escobar-Zepeda et al., 2016). Anaerobic bacteria, such as Prevotella and Faecalibacterium, that are more typically associated with gastrointestinal environments, were for the first time identified in Irish artisanal cheeses at low levels in 2012 (Quigley et al., 2012). Some other recent "firsts" include the detection of Idiomarina and two candidate divisions, namely GNO2 and TM7, in brine samples. Of these, Idiomarina is consider to be a contaminant of brine, with salt acting as the source for this bacteria (Marino et al., 2017). Tremellomycetes, also known as jelly fungi, were also recently found on wooden shelves used for ripening (Guzzon et al., 2017). While the importance of many of these microbes in the context of cheese microbiology is not clear, the fact that HTS has detected the presence of these and other previously-overlooked taxa means that further investigations of such taxa in cheese can now take place.

### BIOTYPE DIVERSITY AND FUNCTIONAL POTENTIAL OF CHEESE

Assessments of biotype diversity among a species enable researchers to detect variation at sub-species level. This can be achieved by HTS analysis of species-specific genes. One approach to determine biotype diversity involved analyzing a species specific gene with sequence heterogeneity across three Italian cheeses, i.e., Mozzarella, Grana Padano, and Parmigiano Reggiano. S. thermophilus was the most abundant bacterium in the three cheeses and, to assess species level diversity, the lacS gene, encoding a lactose permease, was selected for specific analysis. Ultimately, 28 different sequence types were identified among which 13 were present at relative abundance > 1%. Mutations at 60 positions identified in the promoter region upstream of lacS allowed for further differentiation among the 28 sequence types (De Filippis et al., 2014). A similar sub-species differentiation of S. thermophilus has been performed, but instead using the phosphoserine phosphatase gene (serB) as a target sequence. In this case, the approach highlighted the presence of six different sequence types in an undefined milk starter culture (Parente et al., 2016).

Shotgun metagenomic sequencing followed by bioinformatic analysis provides the opportunity to study specific pathways encoded within the cheese microbiome. Ultimately, this technology provides an insight into the functional potential of cheese microbes by studying genes that encode enzymes involved in the catabolic and conversion reactions of amino acids relating to the development of flavor, aroma, and a broad range of other features of relevance to cheese physics, chemistry, and biology. With regard to microbial function or functional potential, spatial distribution analysis of the metabolically active microbiota of three Italian PDO cheeses, Fiore Sardo, Pecorino Siciliano, and Pecorino Toscano, showed a correlation between mesophilic lactobacilli (L. plantarum) and secondary proteolysis as well as the synthesis of volatile components (VOC), such as esters, alcohols, aldehydes, and sulfur compounds. Thermophilic LAB found in Pecorino Siciliano cheese correlated with total free amino acid (FAA) concentrations and Brevibacterium present on the surface of Pecorino Toscano were found to be related to the synthesis of volatile sulfur compounds, such dimethyl trisulphide, dimethyl disulphide, methional, and S-methyl thioesters (De Pasquale et al., 2016). Shotgun sequencing also revealed a strong correlation between bacteria present in commercially available smear cultures and VOC in surface ripened cheeses. In these instances B. linens, Geotrichum candidum, and Staphylococcus xylosus correlated with levels of sulfur compounds and 2-methyl-1-butanol, Debaryomyces hansenii correlated with levels of alcohols and Glutamicibacter arilaitensis with ketones, alcohols, and acids (Bertuzzi et al., 2018). Similarly, a study on Kazak, an artisanal cheese made from fermentation of fresh cow milk through a process involving goat skin bags, identified positive and negative correlations between the core microbiota and amino acid, VOC, and fatty acid concentrations. Kazak cheese analysis revealed that the amino acids, glutamic acid (Glu), histidine (His), isoleucine (Ile), and proline (Pro) had higher positive correlation with one or more abundant bacteria, i.e., Acetobacter, Lactococcus, Bacillus, Staphylococcus, Kurthia, and Moraxella. Among fungi, Dipodascus correlated with concentrations of 9-octadecenoic acid. Pichia and Penicillium correlated with (Z)-9-hexadecenoic acid, Issatchenkia and Candida were correlated with leucine (Leu) and phenylalanine (Phe) levels. Various correlations with VOCs were also evident (Zheng et al., 2018). Through harnessing of this knowledge it may be possible to enhance flavor and fat content by choosing specific strains for cheese manufacture. Another multiomic study focussed on Canestrato Pugliese, an Italian PDO cheese made from raw ewe's milk. HTS revealed that Lactococcus dominated within 3 days of ripening and showed a high capacity to degrade carbon sources. Over the ripening period mesophilic lactobacilli, responsible for proteolysis, increased. As part of the study, the authors developed a model system using HTS together with the Biolog Eco-microplate phenotype method that was designed to select adventitious non-starter lactic acid bacteria (NSLAB) to use as adjunct cultures that guarantee high quality standards for traditional cheese (De Pasquale et al., 2014). A shotgun metagenomic sequencing approach was also employed in the pioneering study of washed rind cheese samples from Europe and USA referred to previously. The metabolism of cysteine and methionine produce volatile sulfur compounds such as methanethiol, and the degradation of valine, leucine, and isoleucine give these cheeses a sweaty putrid aroma. Bloomy and washed rind cheeses were found to contain the halotolerant γ-proteobacteria, Pseudoalteromonas, a bacterial genus originally associated with marine environments. Pseudoalteromonas were found to possess a gene predicted to encode methionine-gammalyase (MGL), an enzyme that converts L-methionine to the aforementioned methanethiol. Among cheese microbes, this enzyme had only been found in Brevibacterium linens previously. Pseudoalteromonas spp. produce cold-adapted enzymes that participate in lipolysis and proteolysis, a feature that is regarded as beneficial in cheese aged and stored at low temperatures as it leads to the development of flavor compounds (Wolfe et al., 2014). Finally, with regard to flavor, pathways involved in the production of flavor and aroma compounds were captured through metagenomic analysis of Cotija cheese (Escobar-Zepeda et al., 2016). The microbes present encode enzymes with the potential to catabolise phenylalanine and branched chain amino acids to yield flavor compounds. The study also detected genes encoding enzymes involved in free fatty acids catabolism, such as carbonyl reductase I and phenol-2 monoxygenase (an enzyme responsible for the production alcohol from xylene) which results in products that impact a characteristic aroma to traditional cheeses. The genetic determinants for antimicrobial compounds such as bacteriocins produced by strains of LAB, which can play an important role in food biopreservation, were also present (Escobar-Zepeda et al., 2016).

Finally, shotgun metagenomic analysis of cheeses has also been employed to determine the microbial basis for cheese colors. Such analyzes of cheeses with a pink discoloration defect revealed the presence of genes that codes metabolic compounds responsible for the production of carotenoids, that were absent from control cheeses. The presence of this pathway was associated with microbes from the genus Thermus, which are renowned for their ability to produce carotenoids but had not been thought to be important in the industrial context of the cheese microbiota previously (Quigley et al., 2016). The aforementioned Bertuzzi et al. study also revealed that surface ripened cheese made from commercially available smear cultures have a greater potential for the production of carotenoids responsible for red/orange colors on the surface of cheeses, highlighting the key role of these smear culture-associated microbes in color development (Bertuzzi et al., 2018).

The next logical step beyond shotgun metagenomics (to assess functional potential) is the use of metatranscriptomics to identify genes that are actually expressed. Metatranscriptomic analysis of Camembert cheese, in which G. candidum and Penicillium camemberti play key roles in ripening, and the de novo assembly of the reads was employed to determine at what stage of cheese ripening are pathways associated with the production of flavor, aroma, and texture development expressed. The study found that these fungal species, already known to contribute to the appearance of the cheese, also play a key role in the development of sensory characteristics of Camembert (Lessard et al., 2014). Metatranscriptomics has also been employed to assess the impact of environmental conditions on gene expression in Caciocavallo Silano, a traditional Italian cheese. These investigations arose from observations that increasing the temperature and decreasing the humidity of the ripening room enhanced the growth of NSLAB in the cheese. Cheese ripened at a higher temperature for 20 days had a similar metabolic profile to that of cheese ripened at standard conditions over a long period, highlighting that increase in ripening temperature results in the rapid maturation of a cheese of quality, while helping the manufacturer to reduce operating costs and to increase the turnover of ripening rooms. The explanation for this reduced ripening time became apparent when it was established that higher temperatures boosted the expression of genes involved in proteolysis, lipolysis, amino acid and fatty acid catabolism and, thus, promoted the VOC related to flavor and aroma (De Filippis et al., 2016). Another such study focussed on a Reblochon-style French cheese, rind activity in which is mostly influenced by S. thermophilus, L. delbrueckii spp.bulgaricus, G. candidum, D. hansenii, and Brevibacterium aurantiacum. Metatranscriptomic analysis of this cheese during ripening showed the up-regulation of S. thermophilus genes associated with protein folding, sorting, degradation, and signal transduction. Transcription, translation as well as nucleotide and lipid metabolism was high among L. delbrueckii. Among the yeasts (D. hansenii, G. candidum) for G. candidum a large increase in the transcription of genes associated with carbohydrate, lipid, and amino acid metabolism was evident, for D. hansenii an increase in genes associated with the metabolism of other amino acids were noted. A decrease in the expression of gene associated with nucleotide metabolism was observed for both yeasts (Monnet et al., 2016). Metatranscriptomics can also be combined with other technologies to provide an even greater insight into the microbiology of cheese. Indeed, a combination of biochemical, metagenomic, and metatranscriptomic analysis has been used to study the microbiology of surface ripened cheese. In this particular study it became clear that dominant bacteria such as L. lactis, Kluyveromyces lactis, D. hansenii, G. candidum, Corynebacterium casei, and Hafnia alvei showed high levels of expression of transcripts related to primary (lactose metabolism, lipolysis, proteolysis) and secondary (catabolism of free amino acids and fatty acids) biochemical reactions (Dugat-Bony et al., 2015). It is anticipated that metatranscriptomic-based analysis of the cheese ripening process will become increasingly common in the near future.

### HTS REVEALS THE PRESENCE OF PATHOGENIC AND SPOILAGE BACTERIA IN CHEESE

HTS can also be employed to detect pathogenic or spoilage microbes in cheeses, milk, or in the production and processing environments. Numerous sources of contamination exist. These include cow feces and teat surfaces (Doyle et al., 2017a; Frétin et al., 2018) as well as bulk tanks and milking machines can also act as sources, such as in the case of the presence of Enterococcus faecalis in raw Herve cheese (Delcenserie et al., 2014).

The consequences of the presence of different microbes are various. The appearance of pink discoloration in cheese manifests through the appearance of pink patches at various locations within the ripened cheese block. As noted above, HTS revealed the presence of spoilage bacteria Thermus thermophilus, a species associated with the production of carotenoids, within defective cheeses (Quigley et al., 2016). The genera Brevibacterium, Corynebacterium, and Microbacterium have been associated with a red-brown defect in smear-ripened cheese, such as Fontina. Deacidification by Debaryomyces facilitates the growth of these bacteria on wooden shelves used for ripening (Guzzon et al., 2017). Analysis of Grana Padano cheeses that exhibited a blowing defect resulted in the detection of seven species of clostridia, i.e., C. sporogenes, C. butyricum, C. disporicum, C. perfringens, C. difficile, C. sordelii, and C. tyrobutyricum. Among these C. tyrobutyricum and C. butyricum were more abundant in blowing defect samples (Bassi et al., 2015). In another study, the presence of genera such as Klebsiella, Morganella, Erwinia, and Acinetobacter, i.e., taxa commonly found in soil and water, in cheese and curd samples of Plaisentif cheese was taken to represent poor hygiene conditions and the contamination of boilers or tools used in the early stages of processing (Dalmasso et al., 2016). Environmental contaminants, including Escherichia sp., Enterobacter cowanii and other Enterobacteriaceae, Agrobacterium sp., Alicyclobacillus sp., and Propionibacterium acnes were also detected in NWC used to make Italian cheeses (De Filippis et al., 2014) and the detection of Pseudomonas was observed in the intermediates of Mozzarella cheese production was suggested to be related to handling and hygiene conditions (Ercolini et al., 2012).

Another mechanism via which microbes can influence cheese quality is through the production of biogenic amines. Histamine, tyramine, putrescine, and cadaverine are the most common biogenic amines that occur in food, with their accumulation causing symptoms such as hypertension, headaches, palpitations, and vomiting in sensitive individuals. The presence of species such as Lactobacillus curvatus, Enterococcus faecium, and Enterococcus faecalis can result in the production of tyramine and Lactobacillus buchneri has been associated with the production of histamine. A HTS-based approach that specifically detects the decarboxylase genes that are responsible for biogenic amine production has been developed (O'Sullivan et al., 2015b).

One obvious concern relating to the use of DNA-based approaches to microbial detection is the possibility of false

### REFERENCES

Alegría, Á., Szczesny, P., Mayo, B., Bardowski, J., and Kowalczyk, M. (2012). Biodiversity in Oscypek, a traditional Polish Cheese, determined by culture-dependent and -independent approaches. Appl. Environ. Microbiol. 78, 1890–1898. doi: 10.1128/AEM. 06081-11

positives arising from the detection of DNA from dead microorganisms. The inclusion of steps involving the use of ethidium monoazide (EMA) IUPAC name: 3-amino-8 azido-5-ethyl-6-phenylphenanthridium bromide or propidium monoazide (PMA) IUPAC name: 3-Amino-8-azido-5-{3- [diethyl(methyl)ammonio]propyl}-6- phenylphenanthridinium to inactivate any DNA which is not sourced from living microbes provides a means of addressing this issue in the context of amplicon based sequencing at least (Rudi et al., 2005; Josefsen et al., 2010; Quigley et al., 2013; Porcellato and Skeie, 2016). Ultimately, while, currently, these investigations have been carried out to identify problematic microorganisms or the source thereof in academic laboratories, it is envisaged that the technology will evolve in a manner that will allow its use in quality assurance laboratories across the food industry.

### CONCLUSIONS

HTS is a powerful technique that can be used to provide a detailed insight into the microbiology of dairy related samples, including raw and pasteurized milk, cheese curd, whey, starter cultures, and cheese. The technology also allows a determination of the impact of seasonal variations, feed given to animals, source of milk, and other environmental factors relating to milk production and processing on the milk and cheese microbiota. Once the microbial community in cheese has been formed, HTS can then be employed to determine the role of the microbial population in factors such as the development of organoleptic properties. While there are very many advantages associated with the use of HTS, some barriers, such as the detection of DNA sequences from non-viable bacteria, cost effectiveness, and experience in bioinformatic analysis need to be overcome before these benefits can be applied extensively across industry. Furthermore, there is a need for the construction of a comprehensive dairy microbial gene catalogs to improve the analysis of HTS reads in the dairy environment (Almeida et al., 2014). However, given that solutions to these issues exist or are in development, it is hoped that any delays to the industrial application of this technology will be short.

### AUTHOR CONTRIBUTIONS

BY and PC prepared the manuscript and all co-authors contributed to editing and critical reviewing thereof.

### FUNDING

This work was funded by the Teagasc Walsh fellowship programme.

gradient gel electrophoresis and pyrosequencing. Int. Dairy J. 21, 142–148. doi: 10.1016/j.idairyj.2010.10.007

Pico cheese (an artisanal Azorean food). Int. J. Food Microbiol. 192, 86–94. doi: 10.1016/j.ijfoodmicro.2014.09.031

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Yeluri Jonnala, McSweeney, Sheehan and Cotter. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# *Zygosaccharomyces bailii* Is a Potential Producer of Various Flavor Compounds in Chinese *Maotai*-Flavor Liquor Fermentation

Yan Xu, Yan Zhi, Qun Wu\*, Rubing Du and Yan Xu\*

State Key Laboratory of Food Science and Technology, Key Laboratory of Industrial Biotechnology, Ministry of Education, Synergetic Innovation Center of Food Safety and Nutrition, School of Biotechnology, Jiangnan University, Wuxi, China

### *Edited by:*

Florence Abram, National University of Ireland Galway, Ireland

#### *Reviewed by:*

Jyoti Prakash Tamang, Sikkim University, India Valentina Bernini, Università degli Studi di Parma, Italy

#### *\*Correspondence:*

Qun Wu<EMAIL_ADDRESS>Yan Xu<EMAIL_ADDRESS>
#### *Specialty section:*

This article was submitted to Food Microbiology, a section of the journal Frontiers in Microbiology

*Received:* 06 July 2017 *Accepted:* 14 December 2017 *Published:* 22 December 2017

#### *Citation:*

Xu Y, Zhi Y, Wu Q, Du R and Xu Y (2017) Zygosaccharomyces bailii Is a Potential Producer of Various Flavor Compounds in Chinese Maotai-Flavor Liquor Fermentation. Front. Microbiol. 8:2609. doi: 10.3389/fmicb.2017.02609

Zygosaccharomyces bailii is a common yeast in various food fermentations. Understanding the metabolic properties and genetic mechanisms of Z. bailii is important for its industrial applications. Fermentation characteristics of Z. bailii MT15 from Chinese Maotai-flavor liquor fermentation were studied. Z. bailii MT15 produced various flavor compounds, including 19 alcohols, six acids, three esters, three ketones, and two aldehydes. Moreover, production of acids and aldehydes were increased by 110 and 41%, respectively, at 37◦C (the maximum temperature in liquor fermentation) compared with that at 30◦C, indicating its excellent flavor productivity. Z. bailii MT15 is a diploid with genome size of 20.19 Mb. Comparative transcriptome analysis revealed that 12 genes related to amino acid transport were significantly up-regulated (2.41- to 5.11-fold) at 37◦C. Moreover, genes ARO8, ARO9, and ALDH4 involved in amino acid metabolism also showed higher expression levels (>1.71-fold) at 37◦C. Increased substrate supply and a vigorous metabolism might be beneficial for the increased production of acids and aldehydes at 37◦C. This work revealed the potential contribution of Z. bailii to various flavor compounds in food fermentation, and produced insights into the metabolic mechanisms of Z. bailii in flavor production.

Keywords: *Zygosaccharomyces bailii,* food fermentation, flavor compounds, genome, transcriptome, Chinese *Maotai*-flavor liquor

## INTRODUCTION

Zygosaccharomyces bailii is a yeast species widely present in various food fermentations, such as wine, tea, and vinegar fermentations (Teoh et al., 2004; Solieri et al., 2006; Garavaglia et al., 2015). In most food and beverage industries, Z. bailii is considered as a problematic spoilage yeast due to its high resistance to preservatives and high tolerance of various stresses (Stratford et al., 2013; Palma et al., 2015). Regardless of its association with spoilage, the potential beneficial effects of Z. bailii have also been proposed in food industries (Ciani et al., 2009; Domizio et al., 2011). In wine fermentation, owing to its high production of esters, Z. bailii in a mixed starter with Saccharomyces cerevisiae improved the production of ethyl esters (Garavaglia et al., 2015). The coculture of Z. bailii with S. cerevisiae also increased the production of polysaccharides that improved the taste and body of wine (Domizio et al., 2011). Additionally, Z. bailii formed part of the tea fungus in Kombucha and Haipao tea fermentation that were rich in crude protein, crude fiber, and lysine (Jayabalan et al., 2010). Nevertheless, the metabolic properties and genetic mechanisms of Z. bailii in food fermentations are still unclear and need to be systematically studied.

Chinese Maotai-flavor liquor is a popular alcoholic beverage even in other parts of the world (Xu and Ji, 2012). This liquor contains over 300 influential flavor compounds including alcohols, acids, esters, ketones, and aldehydes, which greatly contribute to its unique aroma and quality (Xu and Ji, 2012). Maotai-flavor liquor is produced from grains by a spontaneous and solid-state fermentation. Yeasts play essential roles in Maotai-flavor liquor fermentation (Wu, 2013; Wu et al., 2013). Among them, S. cerevisiae is one of the most important and contributes significantly to the quantity and quality of the liquor (Wu et al., 2012; Meng et al., 2015). Z. bailii has been found to be a dominant species in Maotai-flavor liquor fermentation, with proportions close to those of S. cerevisiae (Wu, 2013). Despite its large population, the metabolic activity of Z. bailii in the liquor fermentation is still unclear.

The present work is aimed to study the fermentation characteristics of Z. bailii MT15, including the genome sequencing and comparative transcriptome analysis to unravel its metabolic mechanisms at the relatively high temperature used in liquor fermentation. The results shed new light on the function of Z. bailii and provide a guide for its efficient use in food fermentation.

### MATERIALS AND METHODS

### Yeast Strains

Z. bailii MT15 and S. cerevisiae MT1 were previously isolated from the Maotai-flavor liquor fermentation process and were deposited in the China General Microbiological Culture Collection Center with accession number CGMCC 4745 (Xu et al., 2013), and the China Center for Type Culture Collection with accession number CCTCC M2014463 (Meng et al., 2015), respectively.

### Fermentation Conditions

Sorghum extract was used as fermentation medium and was prepared according to the method below. Two kilogram of ground sorghum was added to 8 L of deionized water. Then the mixture was steamed for 2 h and subsequently saccharified by glucoamylase (5 U/L) at 60◦C for 4 h. The supernatant was collected after being filtered through gauze and centrifuged at 8,000 × g for 15 min. The obtained sorghum extract was diluted with water to produce a final reducing sugar concentration of 75 ± 5 g/L before sterilization (Lu et al., 2015). Fermentation media were prepared with aliquots of 50 mL sorghum extract in 250 mL flasks, and were sterilized at 115◦C for 30 min.

Z. bailii MT15 and S. cerevisiae MT1 were pre-cultured in sorghum extract at 30◦C for 16 h to obtain the seed culture. Then they were inoculated into fermentation media with an initial concentration of 1 × 106 colony-forming units (CFU)/mL. For the comparison of fermentation characteristics between Z. bailii MT15 and S. cerevisiae MT1, the strains were, respectively, fermented at 30◦C for 48 h, with shaking at 200 rpm. To study the effects of temperature on fermentation, Z. bailii MT15 was fermented at 30 and 37◦C for 48 h, with shaking at 200 rpm. During the fermentation, 1 mL of samples were withdrawn at 8-h intervals to determine the cell numbers (**Figure 1A**). Each experiment was performed in triplicate.

### Analytical Determinations and Statistical Analysis

Yeast cells in the seed culture and during the fermentation process were obtained by centrifugation, washed three times with sterile saline solution, diluted to the applicable concentration by saline solution, and counted by hemocytometer under microscope. To analyze the contents of ethanol, flavor compounds, and amino acids, fermentation broths were centrifuged at 8,000 × g for 10 min to remove cells. Ethanol was monitored by high-performance liquid chromatography (HPLC) with a refractive index detector (Varian 355 RI) (Meng et al., 2015). Flavor compounds were analyzed by the gas chromatography-mass spectrometry method (Kong et al., 2014). Free amino acids were detected by HPLC via a pre-column derivatization method (Zhang et al., 2014).

The data are presented in terms of arithmetic averages of three replicates and the error bars indicate the standard deviations. The statistical significance of the difference between the means of samples was tested by one-way analysis of variance (ANOVA) using SPSS Statistics 22.

### Genomic Sequencing and Assembly

The Z. bailii MT15 genome was sequenced using the wholegenome shotgun sequencing approach on an Illumina Miseq platform (Quail et al., 2012). Four paired-end/mate-paired sequencing libraries were constructed with insert sizes of 450, 700 bp, 3, and 8 kb. The de novo assembler Newbler and SSPACE software packages were employed to assemble raw data into contigs and scaffolds (Boetzer et al., 2011; Nederbragt, 2014). The GapCloser program was used to close gaps (Boetzer and Pirovano, 2012). The genome sequence of Z. bailii MT15 was deposited in GenBank under the Whole Genome Shotgun project number SRR5452526.

The genes of Z. bailii MT15 were predicted with the CEGMA pipeline, combining Augustus, SNAP, and Glimmer gene prediction software (Delcher et al., 1999; Stanke and Waack, 2003; Korf, 2004; Parra et al., 2007). The functional annotation of each gene was based on the eggNOG and Swissprot databases (Boeckmann et al., 2003; Powell et al., 2012). The functional classifications were performed with the eggNOG database. Genes for tRNAs and rRNAs were predicted with tRNAscan-SE and RNAmmer 1.2 Server, respectively (Schattner et al., 2005; Lagesen et al., 2007).

### cDNA Preparation and Transcriptome Analysis Using RNA Sequencing

For RNA extraction, Z. bailii MT15 cells were cultivated with three parallel at 30 and 37, respectively, for 24 h. Then, the three parallel were mixed in the same volume and harvested after centrifugation at 8,000 × g for 5 min at 4◦C. Subsequently, the supernatant was removed and the cell pellet was washed three

times with sterile saline solution on ice. Then the washed cell pellet was immediately frozen in liquid nitrogen. Total RNA was isolated using the Trizol Reagent (Invitrogen Life Technologies, Shanghai) according to the manufacturer's instructions. Quality and integrity of total RNA were determined using a Nanodrop spectrophotometer (Thermo Scientific, USA) and Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA). Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA) was used to remove the ribosomal RNAs. The mRNA was subsequently fragmented and used as a template for oligo (dT)-primed PCR.

The cDNA libraries were prepared employing standard techniques for subsequent Illumina sequencing using the mRNA-seq Sample Prep Kit (Illumina, San Diego, CA). The cDNA libraries were sequenced on an Illumina NextSeq 500 according to the manufacturer's instructions. Sequencing raw reads were pre-processed after filtering sequencing adapters, rRNA reads, short-fragment reads, and other low-quality reads. The remaining clear reads were mapped to the reference genome of Z. bailii MT15 using Bowtie2/Tophat2 software (Langmead and Salzberg, 2012; Kim et al., 2015) based on the local alignment algorithm. Gene expression level was normalized by calculating reads per kilobase per million reads (RPKM) (Mortazavi et al., 2008). Differential expression of all of transcripts was quantified using DESeq software, and the method of FDR (False Discovery Rate) control was used to correct the results for multiple hypothesis testing (Anders and Huber, 2013). Significant DEGs were screened based on an FDR threshold of ≤0.05, and a Fold change ≥1.5. The RNA sequence data of Z. bailii MT15 at 30 and 37◦C was deposited in the DNA Data Bank of Japan (DDBJ) with accession IDs DRX082892 and DRX082893, respectively.

### RESULTS AND DISCUSSION

### Metabolic Properties of *Z. bailii* MT15

The biomass and ethanol production of Z. bailii MT15 and S. cerevisiae MT1 was determined during the fermentation process. The cell number of Z. bailii MT15 and S. cerevisiae MT1 was 9.18 × 107 and 8.18 × 107 CFU/mL at the end of fermentation (**Figure 1A**), respectively, which showed no significantly statistical difference (Supplementary Table 1). The highest ethanol production of Z. bailii MT15 was 4.03 g/L, which was 63.76% less than that of S. cerevisiae MT1 (**Figure 1B**, Supplementary Table 2). Considering the large population of Z. bailii MT15, whose proportion can reach 78% of the total yeast population (Wu, 2013), we considered that this yeast also contributes to ethanol production in Maotai-flavor liquor fermentation.

Flavor compounds play essential roles in forming the unique flavor quality of Maotai-flavor liquor (Xu and Ji, 2012). The metabolic activity of Z. bailii MT15 in flavor production was studied at the end of fermentation and compared with that of S. cerevisiae MT1 (**Figure 1C**). Production of alcohols by Z. bailii MT15 was 16,333.78 µg/L, which was approximately half that of S. cerevisiae MT1. Production of ketones by Z. bailii MT15 was 319.95 µg/L, which was also less than that of S. cerevisiae MT1 (452.98 µg/L). However, the production of acids, esters, and aldehydes by Z. bailii MT15 was 1.75, 2.28, and 3.45 times of that of S. cerevisiae MT1, respectively.

Additionally, Z. bailii MT15 produced 38 flavor compounds including 19 alcohols, six acids, three esters, three ketones, two aldehydes, and five other compounds (**Figure 1C**). Among them, 19 flavor compounds were also detected in the fermentation broth of S. cerevisiae MT1. The amounts of 10 flavor compounds, including acetic acid, benzaldehyde, phenethyl acetate, 1-phenylethyl propionate, and acetophenone, produced by Z. bailii MT15 were significantly higher than those produced by S. cerevisiae MT1 (**Figure 1C**). Furthermore, Z. bailii MT15 was able to produce 19 unique flavor compounds compared with S. cerevisiae MT1, including 13 alcohols, two acids, one ketone, one aldehyde, and two other flavor compounds. Some of these unique flavor compounds were influential flavor compounds and contributed significantly to the flavor of the Maotai-flavor liquor, such as 2-heptanol (fruity flavor), 2-nonanol (fruity flavor), 1 nonanol (fruity flavor), 1-hexanol (floral, green scent), geraniol (sweet, rose-like scent), 2-heptanone (fruity, spicy, cinnamon scent), propionic acid (vinegar-like scent), and butyric acid (cheesy-like scent) (Xu and Ji, 2012).

These results demonstrated that Z. bailii MT15 could generate ethanol and various flavor compounds including alcohols, acids, esters, aldehydes, and ketones during liquor fermentation, which would contribute to the flavor and quality of Maotai-flavor liquor. Moreover, previous studies showed that Z. bailii contributed to the flavor complexity of wine and could be used as a mixed starter with S. cerevisiae to improve the production of ethyl esters in wine fermentation (Ciani et al., 2009; Garavaglia et al., 2015). Therefore, except for the capacity of flavor production, Z. bailii would interact with S. cerevisiae during liquor fermentation and positively affect the flavor and quality of Maotai-flavor liquor.

### Effect of Temperature on Flavor Metabolism of *Z. bailii* MT15

Maotai-flavor liquor is produced by a spontaneous fermentation process with a relatively high temperature of up to about 37◦C (Wu et al., 2013). However, little is known about the metabolic properties and mechanisms of Z. bailii under such a relatively high temperature in food fermentation. We therefore studied the metabolic activity of Z. bailii MT15 at 37◦C, and compared it with that at 30◦C.

The biomass of Z. bailii MT15 decreased by 32.46% at 37◦C compared with 30◦C. Meanwhile, few differences were observed in the types of flavor compounds produced, but the amounts of flavor compounds per unit cell were substantially different between the two temperatures (Dataset 1). As shown in **Figure 2**, no significant change was observed in the production of alcohols and esters per unit cell, and the production of ketones decreased by 28.30% at the higher temperature. By contrast, the production of acids and aldehydes per unit cell at 37◦C was, respectively, 110 and 41% higher than that at 30◦C. Maotaiflavor liquor contains large amounts of acids and aldehydes, which have positive effects on its unique sensory characteristics (Fan et al., 2011; Xu et al., 2017). The results suggested that production of acids and aldehydes in Z. bailii MT15 was enhanced under higher temperature. This enhancement would probably contribute to the unique sensory characteristics of Maotai-flavor liquor. Nevertheless, the metabolic mechanisms for its vigorous activities under higher temperature remains unclear. Therefore, we used omics' technology including genomic and transcriptomic analysis to unravel the metabolic features and mechanisms of Z. bailii MT15.

### Genomic and Transcriptomic Analysis of *Z. bailii* MT15

The whole genome of Z. bailii MT15 was sequenced and compared with the genome of three other Z. bailii strains (Z. bailii CLIB 213T , Z. bailii ISA1307, and Z. bailii IST302) (Galeote et al., 2013; Mira et al., 2014; Palma et al., 2017), and the results are shown in **Table 1**. The Z. bailii MT15 genome assembly resulted in 287 contigs (>527 bp) with an N50 value of 156,420 bp, and 95 scaffolds (>2020 bp) with an N50 value of 684,448 bp. The assembled genome was 20.19 Mb with a GC content of 42.38 mol%. The genome size of Z. bailii MT15 was 1.97-fold that of Z. bailii CLIB 213T and 1.87-fold that of Z. bailii IST302, but was similar to that of Z. bailii ISA1307. A total of 9,498 genes were predicted with an average length of 1,461 bp, occupying 68.72% of the whole genome (**Table 1**).

The number of genes in the Z. bailii MT15 genome was nearly twice that of Z. bailii CLIB 213T and Z. bailii IST302,

TABLE 1 | General features of the Z. bailii MT15, Z. bailii CLIB 213T , Z. bailii ISA1307, and Z. bailii IST302 genomes.

NP, Not published.

but was close to the number in Z. bailii ISA1307. It has been proven that Z. bailii ISA1307 is an interspecies hybrid strain that was generated in a stressful environment to improve strain robustness (Sipiczki, 2008; Mira et al., 2014). We analyzed the sequences of housekeeping genes including RPB1, RPB2, TBB, and EFGM in the Z. bailii MT15 genome, which are proposed to have a high capacity to discriminate Zygosaccharomyces species (Suh et al., 2013). The RPB1, RPB2, TBB, and EFGM genes were all duplicated in the Z. bailii MT15 genome, and the sequences of alleles of the four genes were different, which indicated that Z. bailii MT15 might be an interspecies hybrid strain. Furthermore, alleles of the four genes showed high sequence identity to the orthologous genes in Z. bailii CLIB 213T and Zygosaccharomyces parabailii ATCC 60483 (Supplementary Table 3). Thus, it is likely that Z. bailii MT15 is an interspecies hybrid strain generated from Z. bailii and Z. parabailii, which was beneficial for its adaptation to the relatively high temperature environment in Maotai-flavor liquor fermentation.

RNA-Seq was employed to reveal the transcriptomic features of flavor metabolism in Z. bailii MT15 under heat stress at 30and 37◦C. The results showed that 257 genes (2.71% of the Z. bailii MT15 genome) were differentially expressed (≥1.5-fold, P < 0.05), including 126 up-regulated genes (Dataset 2) and 131 down-regulated genes (Dataset 3). These differentially expressed genes (DEGs) were clustered according to eggNOG functional categories (**Figure 3**). Among these categories, general function prediction only (8.56% of DEGs), transcription (7.39% of DEGs), amino acid transport and metabolism (7.00% of DEGs), and carbohydrate transport and metabolism (5.84% of DEGs), contained the greatest number of DEGs (**Figure 3**). Amino acid metabolism is important for the production of ethanol and most flavor compounds in liquor; therefore, we further analyzed the DEGs involved in amino acid transport and metabolism.

Comparative transcriptome analysis showed that 12 genes related to amino acid transport and metabolism were significantly up-regulated at 37◦C (**Figure 4A**). Among them, eight genes encoded amino acid permeases, including one general amino acid permease and seven permeases specific for arginine, proline, lysine, and γ-aminobutyric acid. Higher transcription levels of these genes would be important for the transport of more amino acids into yeast cells (Jauniaux and Grenson, 1990; Regenberg et al., 1999). The other four genes were found to be involved in amino acids biosynthesis, including PTR2 and OPT1 associated with the transport of small peptides, and MEP21 and MEP22 encoding ammonium transporters, would also contribute to a strengthened amino acid metabolism (Meister, 1957). Thus, we speculated that these up-regulated amino acid metabolism related genes would promote the absorption and utilization of amino acid at 37◦C. To validate this speculation, we detected the absorption of the most conventional amino acids by Z. bailii MT15 at 30 and 37◦C after 24 h of cultivation. As shown in **Figure 4B**, the absorption of 14 amino acids was obviously improved by 83% at 37◦C, including the uptake of glutamate acid, aspartate, proline, alanine, and arginine. Among these, the uptake of glutamate acid increased most prominently from 1.14 × 10−9 to 2.11 × 10−9 mg/CFU. By contrast, the absorption of valine and cysteine decreased, and no significant difference was found in the phenylalanine absorption. Therefore, we can conclude that Z. bailii MT15 could absorb and utilize amino acids more efficiently at 37 than at 30◦C, due to the strengthened amino acid associated pathway at higher temperature. The strengthened absorption and utilization of amino acids possibly contributed to the improvement of the corresponding flavor compounds production, in particular, the production of acids and aldehydes.

### Metabolic Mechanisms of Production of Flavor Compounds in *Z. bailii* MT15

Z. bailii MT15 could produce various flavor compounds, especially alcohols, acids, and esters (**Figure 2**). During food fermentation, alcohols and acids are generated by yeasts from amino acids through the Ehrlich pathway and from

sugars through the Harris pathway (**Figure 5**; Ehrlich, 1907; Chen, 1978). To unravel the molecular mechanisms of flavor metabolism in Z. bailii MT15, we compared its genome with that of strain S. cerevisiae MT1 and analyzed the transcription levels of genes involved in flavor metabolism at 37 and 30◦C.

As shown in **Figure 5**, sugars could be converted to αketo acids via the glycolysis and tricarboxylic acid (TCA) cycle pathways. In the Z. bailii MT15 genome, the copy number of genes involved in glycolysis and the TCA cycle was around twice those in the S. cerevisiae MT1 genome (Supplementary Table 4), which was expected as Z. bailii MT15 was considered to be an interspecies hybrid strain. In addition, α-keto acids could also be generated by various aminotransferases through the Ehrlich pathway (Hazelwood et al., 2008). In the Z. bailii MT15 genome, at least 16 different genes were annotated as aminotransferases including AGX1, BCA1, ARO9, ARO8, and YGD3. (Supplementary Table 5). Among them, ARO9 and ARO8

were uniquely found in the Z. bailii MT15 genome and might be associated with the differences in flavor metabolism compared with S. cerevisiae MT1.

The α-keto acids were further decarboxylated to the corresponding aldehydes by α-keto acid decarboxylases including PDC1, PDC5, PDC6, ARO10, and THI3 (ter Schure et al., 1998; Dickinson et al., 2003; Vuralhan et al., 2005). Among them, PDC1, ARO10, and THI3 existed in both of the two strains, while PDC5 and PDC6 were uniquely found in the S. cerevisiae MT1 genome. Pyruvate decarboxylases (PDCs) participate in alcoholic fermentation by converting pyruvate to acetaldehyde and in amino acid metabolism through the Ehrlich pathway (ter Schure et al., 1998). We found 13 genes encoding PDC enzymes in the Z. bailii MT15 genome, whereas only six were found in the S. cerevisiae MT1 genome. Since the ethanol production of Z. bailii MT15 was less than that of S. cerevisiae MT1 (**Figure 1B**), more PDC genes in Z. bailii MT15 did not favor its alcoholic fermentation, but might be beneficial for the conversion of amino acids and reducing sugars to aldehydes by the Ehrlich pathway.

Aldehydes could subsequently be converted to higher alcohols and acids by alcohol dehydrogenases and aldehyde dehydrogenases, respectively (Hazelwood et al., 2008). So far, at least 16 genes encoding alcohol dehydrogenases have been found to catalyze the interconversion of aldehydes and alcohols (Hazelwood et al., 2008). However, only four of these genes (ADH1, ADH4, ADH6, and AAD16) were found in the Z. bailii MT15 genome, while 11 such genes were present in the S. cerevisiae MT1 genome. Z. bailii MT15 harboring fewer alcohol dehydrogenases genes may be the reason that it produced more aldehydes and less alcohols compared with S. cerevisiae MT1. Moreover, the genome analysis showed that three genes (ALDH2, ALDH4, and ALDH5) annotated as aldehyde dehydrogenase were found in the Z. bailii MT15 genome, and ALDH2 was a specific gene compared with S. cerevisiae MT1. This would account for the higher acid productivity of Z. bailii MT15. Therefore, we can conclude that genome differences, particularly in the aspect of genes involved in amino acid metabolism between Z. bailii MT15 and S. cerevisiae MT1, would lead to the different metabolic features of the two strains (**Figure 1C**).

The transcriptome comparison of Z. bailii MT15 revealed that genes involved in the glycolysis pathway showed lower RPKM values at 37 than those at 30◦C, which indicated that the glycolysis pathway might not be related to the increased production of acids and aldehydes at higher temperature. Moreover, in the TCA cycle, ACON2 associated with the production of α-ketoglutaric acid was significantly up-regulated at 37◦C, which would promote production of α-keto acids. Meanwhile, the RPKM values of ARO8 and ARO9, which are involved in the transamination step of the Ehrlich pathway and regarded as broad-substrate-specificity aminotransferases (Iraqui et al., 1998), were higher at 37 than those at 30◦C (≥1.50-fold, P ≤ 0.25). The higher expression of these genes might be beneficial for the transamination of amino acids to corresponding α-keto acids at 37◦C. Furthermore, ALDH4 also showed higher RPKM values at 37◦C (1.78-fold, P = 0.19), which was consistent with the increased production of acids at 37◦C (**Figure 2**). Therefore, the higher expression of genes, including ACON2, ARO8, ARO9, and ALDH4 in the Harris pathway and Ehrlich pathway, might be beneficial for the increased production of acids and aldehydes at the relatively high fermentation temperature.

### CONCLUSIONS

This study revealed the fermentation characteristics and potential function of Z. bailii MT15 that produces various flavor compounds including alcohols, acids, esters, aldehydes, and ketones. Its ability to generate acids and aldehydes is improved at the relatively high temperature used in liquor fermentation,

### REFERENCES

which is beneficial for the complexity of the aroma and quality of the liquor. The genome and transcriptome analysis of Z. bailii MT15 revealed that amino acids metabolism plays important roles in flavor production. This work sheds new light on the metabolic characteristics of Z. bailii in flavor production during Maotai-flavor liquor fermentation that would be applicable to various food fermentations.

### AUTHOR CONTRIBUTIONS

YX (first author), QW, and YZ drafted the manuscript. YX (first author), QW, and RD performed the physiological studies and genome sequencing and transcriptome analysis. QW and YX (last author) participated in the design of the study. All authors read and approved the final manuscript.

### FUNDING

This work was supported by the National Natural Science Foundation of China (31371822, 31530055), the National Key R&D Program (2016YFD0400503), the Priority Academic Program Development of Jiangsu Higher Education Institutions, and the 111 Project (no. 111-2-06).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2017.02609/full#supplementary-material

during kombucha fermentation. Food Sci. Biotechnol. 19, 843–847. doi: 10.1007/s10068-010-0119-6

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Xu, Zhi, Wu, Du and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Novel RAYM\_RS09735/RAYM\_RS09740 Two-Component Signaling System Regulates Gene Expression and Virulence in *Riemerella anatipestifer*

Ying Wang1 , Ti Lu1 , Xuehuan Yin1 , Zutao Zhou1, 2, Shaowen Li 1, 2, Mei Liu1, 2, Sishun Hu1, 2 , Dingren Bi 1, 2 and Zili Li 1, 2 \*

*1 Department of Preventive Veterinary Medicine, College of Veterinary Medicine, Huazhong Agricultural University, Wuhan, China, 2 State Key Laboratory of Agricultural Microbiology, Huazhong Agricultural University, Wuhan, China*

### *Edited by:*

*Florence Abram, NUI Galway, Ireland*

#### *Reviewed by:*

*Dave Siak-Wei Ow, Bioprocessing Technology Institute (A*∗*STAR), Singapore Robin Anderson, Agricultural Research Service (USDA), USA*

> *\*Correspondence: Zili Li<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Systems Microbiology, a section of the journal Frontiers in Microbiology*

*Received: 23 January 2017 Accepted: 04 April 2017 Published: 21 April 2017*

#### *Citation:*

*Wang Y, Lu T, Yin X, Zhou Z, Li S, Liu M, Hu S, Bi D and Li Z (2017) A Novel RAYM\_RS09735/RAYM\_ RS09740 Two-Component Signaling System Regulates Gene Expression and Virulence in Riemerella anatipestifer. Front. Microbiol. 8:688. doi: 10.3389/fmicb.2017.00688*

The Gram-negative bacterium *Riemerella anatipestifer* is an important waterfowl pathogen, causing major economic losses to the duck-producing industry. However, little is known of the virulence factors that mediate pathogenesis during *R. anatipestifer* infection. In this study, RAYM\_RS09735 and RAYM\_RS09740 were predicted to form a two-component signaling system (TCS) through bioinformatics analysis. This TCS was highly conserved across the *Flavobacteriaceae*. A mutant YM1RS09735/RS09740 strain was constructed to investigate the role of the RAYM\_RS09735/RAYM\_RS09740 TCS in *R. anatipestifer* virulence and gene regulation. The median lethal dose (LD50) of YM1RS09735/RS09740 was found to be >1011 CFU, equivalent to that of avirulent bacterial strains. The bacterial abundances of the YM1RS09735/RS09740 strain in the heart, brain, liver, blood, and spleen were significantly lower than that of the wild-type *R. anatipestifer* YM strain. Pathological analysis using hematoxylin and eosin staining showed that, compared to the wild-type, the mutant YM1RS09735/RS09740 strain caused significantly less virulence in infected ducklings. RNAseq and real-time PCR analysis indicated that the RAYM\_RS09735/RAYM\_RS09740 TCS is a PhoP/PhoR system. This is a novel type of TCS for Gram-negative bacteria. The TCS was also found to be a global regulator of expression in *R. anatipestifer*, with 112 genes up-regulated and 693 genes down-regulated in the YM1RS09735/RS09740 strain (∼33% genes demonstrated differential expression). In summary, we have reported the first PhoP/PhoR TCS identified in a Gram-negative bacterium and demonstrated that it is involved in virulence and gene regulation in *R. anatipestifer*.

Keywords: *Riemerella anatipestifer*, *RAYM\_RS09735/RAYM\_RS09740*, two-component signaling system, virulence, RNAseq

## INTRODUCTION

The disease riemerella anatipestifersis, caused by the Gram-negative bacterium Riemerella anatipestifer, occurs primarily in 1–8-week-old ducks but is most common in more susceptible 2– 3-week-old ducklings. It is currently the most economically damaging bacterial infection affecting the global duck industry. Symptoms are characterized by fibrinous pericarditis, glissonitis, airbag inflammation, and meningitis (Segers et al., 1993). A total of 21 R. anatipestifer serotypes have been identified, with no significant cross-protection reported, making it difficult to control the disease through vaccination (Pathanasophon et al., 1995, 2002). R. anatipestifer serotypes 1, 2, and 10 are responsible for most major outbreaks in China. In recent years, many R. anatipestifer strains with differing virulence have been isolated from duck farms in China (Yuan et al., 2011; Wang et al., 2015; Zhang et al., 2015; Song et al., 2016).

Bacteria sense and adapt to their environment via twocomponent signaling systems (TCSs). A paired TCS typically has a sensing histidine kinase (HK) coupled to a response regulator (RR). The sensing of a signal by the HK leads to autophosphorylation on a histidine residue. Subsequent transfer of the phosphate to an aspartate residue on the cognate RR facilitates the binding of the RR to its specific DNA. Each phosphorylated RR regulates specific genes that enable individual bacteria to sense environmental factors and respond to stresses (Stock et al., 2003; Mike et al., 2014). TCSs are present in nearly all sequenced bacterial genomes, as well as in some fungal, archaeal, and plant species (Skerker et al., 2008).

TCSs are involved in the regulation of a variety of important biological functions in bacteria. In particular, PhoP/PhoQ TCSs control the transcription of key virulent genes in several bacterial pathogens, including Escherichia coli, Shigella spp. (Tobe, 2008), Yersinia pestis (O'Loughlin et al., 2010), and Salmonella typhimurium (Tran et al., 2016). TCSs also control the expression of Pneumococcal surface antigen A (PsaA), subsequently regulating virulence and resistance to oxidative stress in Streptococcus pneumonia (McCluskey et al., 2004). Related two-component system also regulate genes that are essential to virulence and complex lipid biosynthesis. For example, SenX3, and RegX3 form a TCS that is involved in Mycobacterium tuberculosis virulence (Ryndak et al., 2008). The TCS is expressed during phosphate starvation and is required for phosphate uptake and aerobic respiration (Parish et al., 2003; Haydel et al., 2012). Additionally, in TB, the TCS PrrAB is required during early intracellular infection (Haydel et al., 2012) and the TCS MprAB responds to envelope stress by regulating stress-responding and virulence-associated genes (He et al., 2006; Pang et al., 2007).

The bacterium R. anatipestifer YM was isolated in Yunmeng, Hubei province, China, and is a highly virulent strain classified as serotype 1 (Zhou et al., 2011). We have previously used in vivo-induced antigen technology (IVIAT) to identify in vivoinduced protein antigens from R. anatipestifer. This predicted the involvement of a putative TCS. To research the function of this TCS in R. anatipestifer, we constructed a mutant with the putative TCS genes, RAYM\_RS09735 and RAYM\_RS09740, deleted to investigate their biological characteristics. We found that RAYM\_RS09735/RAYM\_RS09740 form a PhoP/PhoR TCS, the first reported such TCS in Gram-negative bacteria. We confirmed that the RAYM\_RS09735/RAYM\_RS09740 TCS is an important global transcription regulator and regulates the expression of virulence-associated genes in R. anatipestifer. This may provide the theoretical basis for further study into the molecular pathogenesis of R. anatipestifer and facilitate the design of genetically engineered vaccines against R. anatipestifer.

### MATERIALS AND METHODS Analysis of RAYM\_RS09735 and RAYM\_RS09740 Homology in *R. anatipestifer*

The bacterial strains and plasmids used in this study are listed in **Table 1**. Primers used in this study are listed in **Table 2**. In vivo-induced antigen technology (IVIAT) was

#### TABLE 2 | Primers used in this study.

Wang et al. PhoP/PhoR Regulation in *R. anatipestifer*

then used to characterize potential virulence factors that are expressed in ducks during infection with R. anatipestifer. A genomic DNA library of R. anatipestifer was screened, demonstrating in vivo-induced increased expression of genes with two ORFs, RAYM\_RS09735 and RAYM\_RS09740 (data not shown). RAYM\_RS09735 and RAYM\_RS09740 were predicted to be a histidine protein kinase (HK) and a response regulator (RR), respectively. To investigate RAYM\_RS09735 and RAYM\_RS09740 homology among different R. anatipestifer strains, the open reading frames of RAYM\_RS09735 and RAYM\_RS09740 were amplified from the wild-type R. anatipestifer YM strain. The PCR product was cloned into pMD-18T vectors and sequenced. Homologous amino acid sequences were identified by searching the GenBank database using BLASTX. The resulting alignments were used for the construction of a phylogenetic tree using neighbor-joining (NJ) and were further analyzed using MEGA v6.06 software.

### Bacterial Strains and Animals

R. anatipestifer RA-YM was grown in trypticase soy broth (TSB) or on agar plates (Difco Laboratories, USA) at 37◦C with 5% CO2. Escherichia coli X7213 was cultured in Luria Bertani (LB) broth containing 50 µg/mL diaminopimelic acid (DAP), with shaking at 37◦C overnight (Roland et al., 1999). When needed, antibiotics were used at the following concentrations: 50 µg/mL spectinomycin (Spec), 25 µg/mL chloramphenicol (Cm), and 100 µg/mL ampicillin (Amp). A portion of the bacterial colonies grown in LB medium were stored at −80◦C with 15% glycerol.

### Plasmid Construction

The suicide vector used to create the R. anatipestifer YM RAYM\_RS09735/RAYM\_RS09740 mutant was based on pRE112. Briefly, sequences 800 bp upstream and 800 bp downstream of RAYM\_RS09735 and RAYM\_RS09740 were amplified by PCR using the primer pairs RAYM\_RS09735/RAYM\_RS09740L F/R and RAYM\_RS09735/RAYM\_RS09740R F/R. These contained KpnI and SacI restriction sites. A spectinomycin resistance (SpecR ) cassette (1,185 bp) was PCR amplified from plasmid pIC333 using the primers SpecR F and SpecR R. These three fragments were then purified from an agarose gel and used as a PCR template at a 1:1:2 M ratio to join overlapping PCR products with the primers RAYM\_RS09735L-F and RAYM\_RS09740R-R. The final product was digested with KpnI and SacI enzymes and ligated into the pRE112 plasmid (Edwards et al., 1998) to yield suicide plasmid pRE112-LSR. This plasmid was used to delete RAYM\_RS09735 and RAYM\_RS09740.

### Construction of the *R. anatipestifer* RAYM\_RS09735/RAYM\_RS09740 Mutant

The R. anatipestifer YM RAYM\_RS09735/RAYM\_RS09740 mutant, termed YM1RS09735/RS09740, was constructed via allelic exchange using the previously constructed suicide plasmid pRE112-LSR. The donor strain X7213 was transformed with pRE112-LSR and grown overnight in LB medium supplemented with 50 µg/mL DAP, 50 µg/mL Spec, and 25 µg/mL Cm. The recipient R. anatipestifer YM strain was cultured in TSB medium to an OD600 of 0.4–0.5. One milliliter of the donor strain and 3 mL of receptor strain were centrifuged at 3,000 rpm/min for 5 min, and the pellet re-suspended in 1 mL TSB. Mixed cultures were then incubated on a TSA plate supplemented with 50 µg/mL DAP at 37◦C for 24 h, facilitating the conjugation of the donor and receptor strains. Cultures were next streaked onto a TSA plate containing 50 µg/mL spectinomycin (Spec) to isolate the putative R. anatipestifer YM conjugants from the mixed strains. Single colonies were re-purified on TSA plates supplemented with 50 µg/mL Spec. The R. anatipestifer RAYM\_RS09735/RAYM\_RS09740 mutant strain was screened and validated by PCR.

### Characteristics of the YM1RS09735/RS09740 Strain

The wild-type strain R. anatipestifer or the mutant strain YM1RS09735/RS09740 was grown in TSB at 37◦C for 12 h with shaking, respectively. Equal amounts of YM culture were transferred into fresh TSB (without serum) at a ratio of 1:100 (vol/vol) and incubated at 37◦C with shaking at 200 rpm. Bacterial growth was measured as described previously (Hu et al., 2011) by counting the number of bacterial CFUs at 2 h intervals for 14 h.

### Determination of Bacterial Virulence and Survival *In vivo*

To determine whether deletion of RAYM\_RS09735 and RAYM\_RS09740 influenced R. anatipestifer virulence, the median bacterial lethal dose (LD50) of the mutant strain YM1RS09735/RS09740 was determined (Hu et al., 2010). Healthy Cherry Valley ducks were purchased from Chunjiang Duck Company (China) and hosted in an isolated animal room. All animal experiments and procedures were approved by the Research Ethics Committee, Huazhong Agricultural University, Hubei, China. A total of 75 ducks were divided randomly into 15 groups (five ducks per group). The respective bacterial strains were then injected into duck flippers at a dose of 105–1011 CFU. Ducks in groups 1–7 were injected with 105–1011 CFU of YM bacteria, and ducks in groups 8–14 were injected with 105–1011 CFU of YM1RS09735/RS09740 strain (Zou et al., 2015). The control group was injected with an equivalent volume of PBS.

For pathological histological examination, ducks from each group were sacrificed 48 h after injection with 5 × 106 CFU of each bacterial strain diluted in 0.5 mL PBS. Heart, liver, brain, and spleen samples were collected. Tissues were then fixed in formalin, sectioned, and stained with hematoxylin and eosin.

Twelve 12-day-old ducks were randomly divided into three groups (four ducks per group) and infected with either 5 × 106 CFU of the wild-type R. anatipestifer strain or the mutant strain YM1RS09735/RS09740. Blood, heart, liver, brain, and spleen samples were collected 48 h after injection. Tissues were homogenized in 5 mL PBS, serially diluted, and plated on TSA plates to assess the number of viable bacteria.

### Extraction of Total RNA from *R. anatipestifer In vitro*

The mutant strain YM1RS09735/RS09740 or wild-type strain YM was grown in TSB to log phase, respectively, and then harvested in no more than 3 mL culture by centrifugation at 4000–5000 × g for 5–10 min at 4◦C. Total RNA was extracted using a Bacterial RNA Kit (Omega). A NANODROP 2000c (Nanodrop) was used to measure the concentration and quality of bacterial RNA.

### RNAseq Library Construction

A total of 1 µg RNA per sample was used for RNA sample preparations. Sequencing libraries were generated using a NEBNext R UltraTM RNA Library Prep Kit for Illumina R (NEB), following the manufacturer's recommendations. Briefly, mRNA was purified from total RNA using poly-T oligoattached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). The first strand of cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H–). Second strand cDNA synthesis was then subsequently performed using DNA Polymerase I and RNase H. Any remaining overhangs were converted into blunt ends via exonuclease/polymerase. After adenylation of the 3′ ends of DNA fragments, adaptors with a hairpin loop structure were ligated in preparation for hybridization. To select cDNA fragments by length, library fragments were purified with a AMPure XP system (Beckman Coulter). Next, 3 µL USER Enzyme (NEB) was incubated with size-selected, adaptor-ligated cDNA at 37◦C for 15 min followed by 5 min at 95◦C before PCR. PCR was performed using a Phusion High-Fidelity DNA polymerase, a universal PCR primer set and an Index (X) primer. Finally, PCR products were purified and the library quality assessed using an Agilent Bioanalyzer 2100 system (Agilent).

### RNAseq Differential Expression Analysis

HTSeq v0.6.1 was used to count the number of reads that mapped to each gene and then gene expression was calculated using an RPKM method (Reads Per kb per Million reads). Differential expression analysis was performed using DESeq. A q-value (or FDR) <0.001& and a log2 fold-change >1 were set as the thresholds for determining significantly differential expression.

### Go Enrichment Analysis

Gene Ontology (GO) enrichment analysis of differentially expressed genes was performed using the Bioconductor package GOseq, with a gene length bias correction.  The Gene Ontology database used can be found at http://www.geneontology.org/.

### KEGG Pathway Enrichment Analysis

As different genes cooperate with each other to exercise biological functions, pathway-based analysis can help further understand the biological functions of genes. We used KOBAS to test for the statistical enrichment of differential expression genes in the KEGG pathway data set (http://www.genome.jp/kegg/).

### Quantitative Real-time PCR Analysis

For quantitative real-time PCR (qRT-PCR) validation experiments, 10 genes were randomly selected to assess the RNAseq data (**Table 2**). For this analysis, 1 µg of RNA was reverse-transcribed to cDNA using the PrimeScriptTM RT regent kit with gDNA Eraser (Takara), according to manufacturer's instructions. cDNA was diluted 10-fold and used for real-time PCR analysis using a Bio-Rad CFX96TM System and signal detection protocols in accordance with the manufacturer's instructions (TaKaRa). DnaB was used as an endogenous control. Primers used for the qRT-PCR are described in **Table 2**. Data analysis was performed using GraphPad Prism v 5.0 Software (GraphPad).

### Statistical Analysis and Data Records

Student's t-tests were used to compare gene expression data. Pvalues of ≤0.05 were considered significant. RNAseq original data were uploaded to the NCBI Short Read Archive (SRA) with study accession number SRP096616 (http://www.ncbi.nlm.nih. gov/sra).

### RESULTS

### Identification of the RAYM\_RS09735 and RAYM\_RS09740 Genes in *R. anatipestifer*

Homologous amino acid sequences were identified by searching the GenBank database using BLASTX. A phylogenetic tree was constructed containing 19 amino homologous acid sequences (**Figures 1A,B**). RAYM\_RS09735 and RAYM\_RS09740 share more than 98% sequence identity with seven other R. anatipestifer strains. In addition, RAYM\_RS09735 and RAYM\_RS09740 shared 70% identity with species of Flavobacteriaceae, including Cloacibacterium, Epilithonimonas, and Chryseobacterium. Our results demonstrated that RAYM\_RS09735 and RAYM\_RS09740 are not only highly conserved in R. anatipestifer, but also conserved across the Flavobacteriaceae in general. Functional assessment predicted RAYM\_RS09735 and RAYM\_RS09740 to

domains of RAYM\_RS09735 (C) and RAYM\_RS09740 (D).

be elements of a two-component signaling system (TCS). TCS are typically composed of a sensor with histidine kinase activity and a cytoplasmic transcriptional regulator. RAYM\_RS09735 was identified as a BaeS family histidine kinase, while RAYM\_RS09740 was predicted to be an OmpR family transcriptional regulator (**Figures 1C,D**). Both members of the TCS had the same promoter. Additionally, RAYM\_RS09735 was predicted to be a phosphate regulon sensor protein (PhoR).

### Characterization of a Mutant YM1RS09735/RS09740 Strain

The RS09735/RS09740 gene was deleted from the chromosome of an R. anatipestifer YM strain via allelic exchange. The gene was replaced with a SpecR cassette, allowing successful transfects to be drug selected. The mutant strain was further validated by PCR amplification of RAYM\_RS09735, RAYM\_RS09740, and 16S rRNA fragments from transconjugants (**Figure 2A**). Real-time PCR analysis confirmed that RAYM\_RS09735 and RAYM\_RS09740 transcription was completely abolished in the mutant strain. However, inactivation of the RAYM\_RS09735 and RAYM\_RS09740 genes led to significantly increased transcription of the chromosomally upstream gene RAYM\_RS09730 and the downstream gene RAYM\_RS09745 (**Figure 2B**). The RAYM\_RS09735 and RAYM\_RS09740 deletion mutant strain was designated YM1RS09735/RS09740. Growth curve measurements revealed that growth in TSB was similar between the wild-type YM and mutant strains (**Figure 2C**). Transmission electron microscopy also showed that there were no significant changes in bacterial morphology in the RA-YM RAYM\_RS09735/RAYM\_RS09740 mutant strain (data not shown).

### Pathogenicity of the YM1RS09735/RS09740 Mutant

Bacterial virulence was evaluated by median lethal doses (LD50) in 12-day-old Cherry Valley ducks. The LD50 for the mutant strain was greater than 1011 CFU, which was more than a 103 fold attenuation in virulence compared to the wild-type YM strain (4 × 107 CFU). These LD50 values demonstrate that the strain is almost avirulent in ducklings (**Table 3**). To further investigate the role of RAYM\_RS09735 and RAYM\_RS09740 in in vivo systemic infections, the bacterial loads in the liver, spleen, heart, brain, and blood from infected ducks were quantified. This showed that, in these organs, bacterial capacity was significantly decreased compared to ducks infected with the wild-type strain (**Figure 3**). Pathological histological analysis was also performed 48 h post-infection comparing the wild-type YM and YM1RS09735/RS09740 strains. Hematoxylin and eosin staining revealed that the capacity of the YM1RS09735/RS09740 mutant to damage the heart, brain, and spleen was significantly decreased relative to the wild-type strain, although mild liver damage was still observed (**Figure 4**). These results show that RAYM\_RS09735 and RAYM\_RS09740 significantly affect the virulence of R. anatipestifer.

### Function and GO Enrichment Analysis of the Differentially Expressed Genes

Differentially expressed genes (DEGs) between the mutant YM1RS09735/RS09740 and wild-type strains were identified

YM1RS09735/RS09740 strain; Lane 4, delta genes were amplified from RA-YM strain; Lane 5, SpecR cassette was amplified from YM1RS09735/RS09740 strain; Lane 6, SpecR cassette was amplified from RA-YM strain; Lane 7, 16S rRNA fragment was amplified from YM1RS09735/RS09740 strain; Lane 8, 16S rRNA fragment was amplified from RA-YM strain. (B) Real-time PCR analysis. The flanks mRNA levels of *RAYM\_RS09735/RAYM\_RS09740* genes were measured. The changes of transcription were expressed as fold expression; (C) Bacterial growth curves. The mutant strain YM1RS09735/RS09740 or wild-type strain RA-YM was grown on TSB medium, and growth of each strain was monitored by measuring the CFU/ml. \*\*\**P* < 0.001.

using RNAseq. In total, 805 genes were found to be differentially expressed, with 112 genes upregulated (13.9%), and 693 genes downregulated (86.1%) in the mutant YM1RS09735/RS09740 strain compared to the wild-type strain (**Figure 5A**). Of the ∼2,000 genes identified in the R. anatipestifer genome, more than one third were differentially expressed. Further analysis of the RNAseq data using KEGG pathways predicted that RAYM\_RS09735 and RAYM\_RS09740 are components of a PhoP/PhoR TCS (**Figure 5B**). Gene Ontology (GO) enrichment analysis showed that differentially expressed genes can be found in a wide-variety of biological processes, cellular components, and molecular functions, including 21 distinct pathways (**Figure 6**). RNAseq was validated using qRT-PCR analysis of randomly selected 10 genes. Overall, the changes in expression of these randomly genes tested by qRT-PCR agreed with the direction determined by RNAseq (**Figure 7A**). Linear regression analysis was performed examining the fold-changes of the gene expression ratios between RNAseq

and qRT-PCR, showing a significant positive correlation (**Figure 7B**).

We found that 11 genes had more than 4-fold higher expression in the mutant strain compared to the wild-type, eight of which encode hypothetical proteins. The three other genes encoded a carbohydrate-binding protein, a glycan metabolism protein (RagB), and a phosphate subunit transfer protein (PstS), respectively. We also found that several transcription factors, components of the CRISPR system, and additional putative proteins were also upregulated, with most upregulated genes involved in bacterial metabolism. In addition, 20 genes were found to be downregulated more than 4-fold. Of these, 13 genes encoded hypothetical proteins, while the other seven encoded RAYM\_RS09735, RAYM\_RS09740, a lipoprotein, a peptidoglycan hydrate (Nlp/P60), a DNA-binding protein, von Willebrand factor A, ATPase AAA, and an uncharacterized conserved protein. Other downregulated genes included several hypothetical proteins, transcription factors, and metabolic genes, in addition to multiple molecular chaperones and TonBdependent receptors.

Finally, we performed real-time PCR analysis of three predicted virulence-associated genes expression in the YM1RS09735/RS09740 mutant strain (pstS, BLP, and Nlp/P60). We found that the mRNA expression of pstS was 8-fold higher (**Figure 8A**), consistent with our RNAseq data. Bacterial lipoprotein (BLP) and hydrolase Nlp/P60 were found to be significantly downregulated (**Figure 8B**).

### DISCUSSION

Bacterial two-component regulatory systems (TCSs), consisting of a sensing histidine kinase and a response regulator, mediate gene expression in response to environmental stimuli

molecular-function, which including 21 pathways.

(Stock et al., 2003). To achieve this, histidine kinases sense environmental signals and autophosphorylate. This additional phosphate is then transferred to an aspartic acid residue on the corresponding response regulator (Hoch, 2000; Alm et al., 2006). The phosphorylated response regulator elicits a diverse range of downstream responses, including an enhancement of its DNA binding ability. This allows the TCS to modulate target gene expression (West and Stock, 2001; Lin et al., 2015). Using an in vivo antigen technology (IVIAT) to examine ducks infected with R. anatipestifer, we have identified a putative TCS involving the genes RAYM\_RS09740 and RAYM\_RS09735. Bioinformatics analysis of these two genes, and their proteins, revealed that RAYM\_RS09740 and RAYM\_RS09735 encode a histidine kinase and response regulator, respectively. We also found that, like

FIGURE 7 | Real-time PCR verification of differentially expressed genes in YM1RS09735/RS09740 mutant strain. (A) The fold expression changes of randomly selected 10 genes to validity of the RNAseq results. (B) Coefficient analysis of fold changes data between qRT-PCR and RNAseq. Ten different expression genes were selected from qRT-PCR.

most TCSs, RAYM\_RS09735 and RAYM\_RS09740 have the same promoter and co-transcribe as an operon (Aggarwal et al., 2016). RAYM\_RS09735 and RAYM\_RS09740 possess all of the characteristic TCS domains essential to their biochemical activities and responses. RAYM\_RS09735 and RAYM\_RS09740 were predicted to belong to the OmpR and BaeS families, respectively, and RAYM\_RS09735 was further predicted to be a phosphate regulon sensor protein (PhoR).

Our study aimed to investigate the RAYM\_RS09735/RAYM\_RS09740 TCS and explore its potential functions in R. anatipestifer. To better evaluate the role of RAYM\_RS09735 and RAYM\_RS09740 in R. anatipestifer, we constructed the mutant strain YM1RS09735/RS09740 in which the RAYM\_RS09735 and RAYM\_RS09740 genes were mutated by deleting a 1400-bp fragment from the wild-type RA-YM strain. RT-PCR and qRT-PCR were used to validate the mutant strain (**Figures 2A,B**). We found that upstream and downstream genes were significantly down-regulated in the YM1RS09735/RS09740 mutant, although the strain showed the same growth characteristics as the wild-type in TSB (**Figure 2C**).

LD50 analysis of the two strains demonstrated the RAYM\_RS09735/RAYM\_RS09740 mutant had significantly reduced virulence. Indeed, the strain was almost avirulent to ducklings (**Table 3**). A normal R. anatipestifer infection is characterized by septicemia and the ability to colonize and develop in the tissues of their host. We found that the bacterial loads in the liver, spleen, heart, brain, and blood of ducks infected with the mutant YM1RS09735/RS09740 strain were significantly lower than that of wild-type RA-YM infected ducklings (**Figure 3**). In addition, pathological histological analysis using hematoxylin and eosin staining showed that the mutant strain also induced significantly less damage to the heart, brain, and spleen, compared to the wild-type strain. There was still slight damage to the liver (**Figure 4**). Our results indicate that the RAYM\_RS09735 and RAYM\_RS09740 genes are important mediators of R. anatipestifer virulence during infection in ducklings. Therefore, RAYM\_RS09740 was determined to be a response regulator, with a combined promoter affecting downstream genes. The protein also altered the expression of target genes. Our study therefore indicated that virulence factors in R. anatipestifer could be regulated by RAYM\_RS09740.

RAYM\_RS09735 and RAYM\_RS09740 share up to 70% sequence identity with Flavobacteriaceae species, such as Cloacibacterium, Epilithonimonas, and Chryseobacterium. These species have been previously reported to exhibit drug resistance. Molecular docking analysis identified inhibitors of various RRs in OmpR family to affect this observed drug resistance. Therefore targeting TCSs has become an attractive option for the development of new drugs, particularly to combat M. tuberculosis (Banerjee et al., 2016). These studies also demonstrate that TCSs can be used as the theoretical basis for the study of drug resistance mechanisms in Flavobacteriaceae.

To investigate how RAYM\_RS09735 and RAYM\_RS09740 regulate gene expression, RNAseq analysis was performed to identify differentially expressed genes in the mutant YM1RS09735/RS09740 strain. This revealed 112 genes that were upregulated (13.9%) and 693 genes that were downregulated (86.1%) relative to the wild-type strain (approximately one-third of genes had differential expression). To validate the RNAseq results, we randomly assessed the differential expression of 10 genes using qRT-PCR, confirming the accuracy of the RNAseq data (**Figure 7**). Gene ontology (GO) enrichment analysis showed that differentially expressed genes were involved in a wide-variety of biological processes, cellular components, and molecular functions, including 21 distinct pathways (**Figure 6**). Therefore, we hypothesize that the RAYM\_RS09735/RAYM\_RS09740 TCS is a global expression regulator that controls the expression of virulence genes, thus affecting the pathogenicity of R. anatipestifer.

KEGG pathway analysis showed that RAYM\_RS09740 is a PhoP protein, indicating that RAYM\_RS09735 and RAYM\_RS09740 form a PhoP/PhoR TCS. PhoP/PhoR TCSs have been implicated in several biological processes in Gram-positive bacteria, such as B. subtilis (Hulett et al., 1994), Streptomyces lividans (Sola-Landa et al., 2003), and M. tuberculosis (Perez et al., 2001). Specifically, PhoP/PhoR has been shown to control respiration(Birkey et al., 1998), cell wall metabolism (Minnig et al., 2005), culture metabolism (Thomas et al., 2012), and biofilm formation (Bluskadosh et al., 2013). Most importantly, in many pathogens PhoP/PhoR TCSs have also been shown to regulate pathogenesis (Perez et al., 2001; Ryndak et al., 2008). PhoP and PhoQ TCS components have been reported in many Gram-negative bacteria, including E. coli (Eguchi et al., 2007), Edwardsvilla (Lv et al., 2012), and Salmonella (Tran et al., 2016). PhoP/PhoR is also an important TCS in M. tuberculosis and is a key player in virulence (Ryndak et al., 2008). However, PhoP/PhoR TCSs have not been previously observed in Gram-negative bacteria. Our data shows that a putative PhoP/PhoR TCS is involved in R. anatipestifer virulence, although further phosphotransfer profiling will be required to validate RAYM\_RS09735/RAYM\_RS09740 as a PhoP/PhoR TCS.

Our data also suggest that RAYM\_RS09740 can act as a transcription factor and plays a role in regulating genes expression, although there are likely other TCSs present in R. anatipestifer. This is based on our observation that, after the deletion of RAYM\_RS09735 and RAYM\_RS09740, other TCSs regulated gene expression. We hypothesize that crosstalk and -regulation may occur between a multitude of twocomponent systems in R. anatipestifer, similarly to other organisms. For example, the transcription factor NemR is a redox-regulated transcriptional repressor in E. coli (Gray et al., 2013) and a carbohydrate-responsive regulatory protein (BadR) is a transcriptional repressor of rpoS in Borrelia burgdorferi (Miller et al., 2013).

The PhoP/PhoR TCS we identified was initially predicted to be part of the phosphate (Pi) stress-response in R. anatipestifer. PhoR acts to sense environmental phosphate (Pi) levels in bacteria through the ABC-type phosphate-specific transport (Pst) system and the protein PhoU. PstS is a periplasmic protein that binds Pi with high affinity and PhoP induces the promoter activity of pstS. We found with RNAseq and qRT-PCR that the mRNA expression of pstS was 8-fold higher (**Figure 8A**) in the mutant strain and different from other Gram-negative bacteria. Therefore, RAYM\_RS09735 and RAYM\_RS09740 may act more like the PhoP/PhoR system of M. tuberculosis. Additionally, M. tuberculosis PhoP/PhoR does not respond to phosphate starvation and is similar to the PhoP/PhoQ TCS in Salmonella. In Salmonella PhoP/PhoQ TCS, PhoQ senses low Mg2+ levels, and PhoP activates expression of genes encoding high-affinity Mg2+ transport systems (Walters, 2006). We aim to investigate the environmental factors sensed by R. anatipestifer PhoP/PhoR in future work.

Further, DEG analysis showed that the expression of several other signal transduction system genes was altered in the YM1RS09735/RS09740 strain. It is likely that PhoP can influence the expression of transcription factors in other signal transduction system and there is cross-talk or cross regulation between PhoP/PhoR and other TCSs in R. anatipestifer. As more than one third of genes demonstrated differential expression, it is unlikely that the PhoP/PhoR TCS directly regulates all of these DEGs. However, the TCS we identified may affect wider gene expression by influencing the expression of other transcription factors in a variety of signal transduction systems.

Additionally, the CRISPR-Cas bacterial immune system has been shown to be present in R. anatipestifer (Chamnongpol et al., 2003; Zhu et al., 2016) and we found that the expression of cas9, cas1, and cas2 were all different in the YM1RS09735/RS09740 mutant, compared to the wild-type. Recent studies have reported that quorum sensing has a role in controlling the Pseudomonas aeruginosa CRISPR-Cas adaptive immune system (Høyland-Kroghsbo et al., 2016) and the TCS BfiSR affects the production of proteins involved in virulence, post-translational modification, and quorum sensing (Petrova and Sauer, 2010). We therefore reasoned that TCS might also affect the CRISPR-Cas immune system of R. anatipestifer.

We found 13 genes encoding hypothetical proteins that were downregulated more than 4-fold. At present, M949\_RS01915, M949\_1360, AS87\_01735, AS87\_03730, M949\_1556, and AS87\_04050 are all hypothetical proteins but have been reported to associate with virulence in R. anatipestifer. Bacterial lipoprotein (BLP) and the hydrolase Nlp/P60 were both found to be significant down-regulated in the mutant YM1RS09735/RS09740 strain (**Figure 8B**). Nlp/P60 is involved in cell growth and division, autolysis, and invasion (Xu et al., 2015). Bacterial lipoproteins on the cell surface include an important class of virulence factor that is expressed by many pathogens. BLP can promote TLR2 expression, TLR2-induced NF-κB activation, and IL-6 production, leading to immune inflammation (Buddelmeijer, 2015). The clinical symptoms of R. anatipestifersis are characterized by fibrinous pericarditis, glissonitis, airbag inflammation, and meningitis. Our pathological results demonstrated that the damage induced by the mutant YM1RS09735/RS09740 strain to the heart, brain, and spleen was significantly less severe than the wild-type strain (**Figure 4**). We therefore hypothesize that R. anatipestifer reduced expression of BLP and the hydrolase Nlp/P60, subsequently affecting inflammation in ducklings.

Other downregulated genes we identified included transcription factors, metabolic genes, multiple molecular chaperones, and TonB-dependent receptors. The primary role of DNA binding proteins and transcription factors is to regulate the expression of the genes and RAYM\_RS09740 may regulate transcription by modifying the expression of these transcription factors. Molecular chaperone proteins assist in the folding of nascent polypeptide chains and prevent the aggregation of denatured proteins (Hartl and Hayerhartl, 2002). We found that transcript for the molecular chaperones dnaJ, dnaK, grpE, clpP, tir, and skp were all downregulated in the mutant YM1RS09735/RS09740 strain, suggesting that the RAYM\_RS09735/RAYM\_RS09740 TCS may regulate the expression of the chaperone genes to influence the severity of R. anatipestifersis. The DnaK molecular chaperone system includes DnaK, DnaJ, and GrpE (Mayer et al., 2000; Ben-Zvi and Goloubinoff, 2001; Genevaux et al., 2007) and is related to stress tolerance and virulence in Salmonella. DnaK can also regulate the expression of virulence islands and a DnaK mutation strain was found to be largely avirulent (Takaya et al., 2004). ClpP was found to be a virulence-related factor in Actinobacillus pleuropneumoniae, with a clpP gene deletion strain demonstrating decreased biofilm production (Xie et al., 2013). We also found that TonB-dependent receptors were

### REFERENCES

downregulated. These receptors have been reported to associate with hemin iron acquisition and R. anatipestifer virulence.

In this study, RAYM\_RS09735 and RAYM\_RS09740 were predicted to be components of a PhoP/PhoR TCS in R. anatipestifer using bioinformatics analysis. Furthermore, deletion of the RAYM\_RS09735 and RAYM\_RS09740 genes significantly decreased R. anatipestifer virulence. RNAseq analysis showed that RAYM\_RS09740 is a global expression regulator in R. anatipestifer. We also demonstrated that the RAYM\_RS09735/RAYM\_RS09740 TCS contributes to virulence, signal translation, and the CRISPR-Cas adaptive immune system of R. anatipestifer. BLP and the hydrolase Nlp/P60 were predicted to be two virulence factors involved. Although common in Gram-positive bacteria, the predicted PhoP/PhoR TCS we have putatively characterized is the first such TCS identified in Gram-negative bacteria. Our study is also the first to report this TCS in R. anatipestifer, providing new insight into the role of TCSs in regulating R. anatipestifer virulence.

### ETHICS STATEMENT

All the animal experiments were carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals from Research Ethics Committee, Huazhong Agricultural University, Hubei, China. All procedures performed in studies involving animals were in accordance with the ethical standards of the institution or practice at which the studies were conducted.

### AUTHOR CONTRIBUTIONS

Funding acquisition: ZL. Investigation: YW. Methodology: YW, TL, and XY. Project administration: SL, ML, SH, and DB. Resources: ZZ. Supervision: ZL. Validation: ZL. Visualization: YW. Writing—original draft: YW. Writing—review's editing: ZL.

### FUNDING

This work was supported by grants from the Wuhan Science and Technology Bureau (No.<PHONE_NUMBER>010070 to ZL) and Natural Science Foundation of Hubei Province (No. 2015CFB268 to ZL).

RstA/RstB and PhoP/PhoQ systems. Biochim. Biophys. Acta 1864, 1686–1695. doi: 10.1016/j.bbapap.2016.09.003

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.