Transcriptomic analysis of thermotolerant yeast Kluyveromyces marxianus in multiple inhibitors tolerance

During pretreatment of lignocellulosic biomass, toxic compounds were released and inhibited the growth and fermentation of microorganisms. Here the global transcriptional response of K. marxianus to multiple inhibitors including acetic acid, phenols, furfural and HMF, at 42 °C, was studied, via RNA-seq technology. Genes involved in the glycolysis pathway, fatty acid metabolism, ergosterol metabolism and vitamin B6 and B1 metabolic process were enriched in the down-regulated gene set, while genes involved in TCA cycle, respiratory chain, ROS detoxification and transporter coding genes were enriched in the up-regulated gene set in response to the multiple inhibitors stress. Further real time-PCR results with three single inhibitor stress conditions showed that different transporters responded quite differently to different inhibitor stress. Coenzyme assay results showed that the level of NAD+ was increased and both NADH/NAD+ and NADPH/NADP+ ratio decreased. Furthermore, genes involved with transcription factors related to carbohydrate metabolism, sulfur amino acids metabolism, lipid metabolism or those directly involved in the transcriptional process were significantly regulated. Though belonging to different GO categories or KEGG pathway, many differentially expressed genes were enriched in maintaining the redox balance, NAD(P)+/NAD(P)H homeostasis or NAD+ synthesis, energy production, and iron transportation or metabolism. These results suggest that engineering these aspects represents a possible strategy to develop more robust strains for industrial fermentation from cellulosic biomass.


Introduction
Lignocellulosic biomass has the potential to contribute substantially to future global energy demands, because of its low cost, large-scale availability, not competing with food production, and high potential to reduce greenhouse gas emission. [1][2][3][4] Pretreatment is essential for releasing fermentable sugars from lignocellulose biomass. However, during harsh pretreatment processes, various small molecules such as weak acids, furan aldehydes, and phenolic compounds (referred to as "fermentation inhibitors"), are generated from partial overdegradation of lignocellulose and inhibit consequent cell growth and microbial fermentations. 5 The inhibitor tolerance is one of the important parameter for effective fermentation process. While prior studies are mostly focused on characterization of genetic mechanisms for yeast stress response to individual inhibitory compounds, cellulosic hydrolysates contain multiple fermentation inhibitors with distinct toxicity mechanisms rather than a single inhibitor. It was found that different mechanisms could be adopted by the yeast to resist hydrolysates inhibitors, e.g. acetic acid, furfural, and 5-hydroxymethylfurfural (HMF). 6 However, there is still limited information on what genetic perturbation targets could be elicited to improve yeast resistance to mixed fermentation inhibitors. Therefore, a better understanding of genetic regulatory networks underlying the resistance to multiple lignocellulosicderived fermentation inhibitors is needed to develop strains with enhanced tolerance to cellulosic hydrolysates. Considering the fact that various inhibitors oen coexist in the hydrolysate and could cooperate with each other to become even more toxic to yeast than existing alone (i.e., synergistic stress), the knowledge on how yeast cells reprogram their metabolism in response to lignocellulosic-derived fermentation inhibitors is of particular interests to biofuel and biochemical production.
Kluyveromyces marxianus is considered as a 'generally regarded as safe' (GRAS) microorganism. Though the genome of K. marxianus was much smaller (less than 5000 open reading frames) 7 than that of S. cerevisiae (over 6000 genes), 8 it has advantages such as short generation time and high growth rate at elevated temperatures (0.86-0.99 h À1 at 40 C), with an upper growth limit of 52 C of some strains. 9 K. marxianus also has the intrinsic fermentation capacity to utilize various substrates including xylose. [10][11][12] Therefore, there are increasing applications of K. marxianus in high temperature fermentation with lignocellulosic hydrolysates. However, the knowledge of its stress physiology is scarce. Moreover, K. marxianus natively exhibited higher assimilation rates for aldehydes such as furfural, HMF, vanillin etc., compared to glucose-fermenting microorganisms such as Klebsiella pneumoniae, Saccharomyces cerevisiae, and Zymomonas mobilis with no genetic modication. 13 Our study also showed that K. marxianus could ferment with non-detoxied diluted acid pretreated corncob to produce ethanol and xylitol and possess considerate inhibitors tolerance especially to furfural and HMF. 14 However, compared with vast information of various inhibitors tolerance in S. cerevisiae, there is very limited information on K. marxianus with the resistance mechanism to the fermentation inhibitors. Therefore, transcriptomic analysis of the tolerance response of lignocellulosic hydrolysates inhibitors or fermentation inhibitors will be much helpful in K. marxianus fermentation study.
Although genome sequences of several K. marxianus strains have been published, [15][16][17] detailed reports on the transcriptional analysis of K. marxianus with various fermentation perturbations are still very limited. Lertwattanasakul et al. conducted transcriptome analyses of K. marxianus DMKU 3-1042 to identify genes related to growth with glucose at 45 C and with xylose at 30 C. Gao et al. reported the transcriptional analysis of K. marxianus for ethanol production from inulin. 18 Up to now, no detailed transcriptional analysis of K. marxianus is available with lignocellulosic-derived fermentation inhibitors at elevated temperature (>30 C). Comparing with the vast transcriptional analysis reports on S. cerevisiae, the study of K. marxianus is very limited which hindered the future development of K. marxianus application in industry.
Here we conducted transcriptomic analysis of K. marxianus at elevated temperature (42 C) with or without three main lignocellulosic-derived fermentation inhibitors including acetic acid, furfural, HMF and phenols by next-generation sequencing technology for RNA (RNA-seq). The transcriptional comparison provides useful information on the molecular basis of genomewide microbial responses to the mixed fermentation inhibitors, including the molecular basis of the central carbon metabolism, mitochondrial respiratory chain, redox homeostasis, MSN2/4 mediated stress response element (STRE)-controlled genes, fatty acid and ergosterol metabolism, alanine, aspartate and glutamate metabolism, vitamin B6 and B1 metabolism, together with various transporters genes which would facilitate the development of K. marxianus in the industrial application. Results of this study will aid dissection of lignocellulosic hydrolysate inhibitors tolerance mechanisms in yeast and metabolic engineering efforts for more tolerant strain development.

Reagents and strains
All chemicals used were of analytical grade or higher. D-Glucose was obtained from Sangon Biotech Co. (Shanghai, China). The yeast extract and peptone were purchased from Oxoid (Oxoid Ltd., Basingstoke, Hampshire, England). K. marxianus YHJ010 was the TRP1, LEU2 and URA3 auxotrophic strain of NBRC1777 (ref. 19) and was cultivated in YPD medium (1% w/v yeast extract, 2% w/v bacto peptone, 2% w/v glucose) at 42 C.

Samples preparation and transcriptome analysis
Cell growth conditions. K. marxianus YHJ010 was precultivated in 5 ml of YPD medium at 42 C overnight. Then the cells were transferred into 500 ml Erlenmeyer asks containing 100 ml of the YPD medium with or without mixed inhibitors containing 0.7 g l À1 furfural, 0.7 g l À1 HMF, 3.0 g l À1 acetate acid and 0.28 g l À1 phenols (4-hydroxybenzaldehyde, syringaldehyde, catechol and vanillin with 0.07 g l À1 each compound) with initial OD 600 of 0.5 and cultivated at 42 C with shaking at 250 rpm in an orbital shaker until OD 600 of 6 (the mid-exponential phase of growth). In the case of individual inhibitor condition, acetic acid stress condition with 2.5 g l À1 acetic acid, furfural stress condition with 1.5 g l À1 furfural, and phenols stress condition with 0.8 g l À1 phenols (4-hydroxybenzaldehyde, syringaldehyde, catechol and vanillin with 0.2 g l À1 each compound), respectively. Yeast cells were then recovered when OD 600 reached 6 and rapidly frozen in liquid nitrogen and kept at À70 C until the subsequent RNA extraction step. Cell growth was monitored by determining the optical density (OD 600 ) with a Hitachi UV-2550 Spectrophotometer.
RNA isolation, preparation of cDNA library and sequencing. Total RNA from yeast was extracted following the standard protocol of TRIZOL Reagent (Invitrogen, Carlsbad, CA, USA). mRNA was isolated from total RNA using Magnetic Oligo-dT beads, fragmented into short fragments and then used to synthesize rst-strand cDNA with random primers. RNase H, buffer, dNTPs and DNA polymerase I (TaKaRa) was used to synthesize the second-strand cDNA. Sequencing adapters were ligated to short fragments and resolved by agarose gel electrophoresis. Suitable fragments were puried and subsequently amplied by PCR to create the cDNA library.
Genome annotation and bioinformatics analysis. Adaptor sequences, empty reads, and low-quality sequences were removed from the raw reads, and the resulting clean reads were mapped to the reference genome of Kluyveromyces marxianus NBRC1777 from GenBank with accession no. AP014599-AP014607 (ref. 7) using TopHat (http://tophat.cbcb.umd.edu). The whole-genome sequences above were annotated according to the Gene Ontology (GO) database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database.
For gene function annotation, obtained unigene sequences were annotated by searching in various protein databases, including the National Center for Biotechnology information (NCBI) nonredundant protein (Nr) database, the NCBI nonredundant nucleotide sequence (Nt) database, Cluster of Orthologous Groups of proteins (COG), Search Tool for the Retrieval of Interacting Genes (STRING), Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). In addition, information for the differentially expressed genes (DEGs) was subjected to GO and KEGG signicant enrichment analyses to identify biological functions and metabolic pathways in which these genes participate.
For differential gene expression analysis, reads per kilobase of exon model per million mapped reads (RPKM) was used as a value of normalized gene expression. Statistical comparison of RPKM values between the samples was conducted using a web tool Cuffdiff (http://cole-trapnell-lab.github.io/cufflinks/ cuffdiff/index.html). Genes were considered differentially expressed in a given library when p-value < 0.05 and a greater than two-fold change in expression across libraries observed.

Real-time PCR analysis
Total RNA was isolated using a yeast total RNA extraction kit (Sangon Biotech Co. Shanghai, China). Isolated RNA was treated with RNase-free DNase I (Toyobo, Osaka, Japan) at 37 C for 15 min to remove the potentially contaminated genomic DNA. cDNA was synthesized by the ReverTra Ace qPCR RT Master Mix kit (Toyobo, Osaka, Japan). The reverse transcription reaction was performed in an Arktik thermal cycler (Thermo Fisher Scientic, Waltham, MA, USA) at 37 C for 15 min, 50 C for 5 min, and denaturing at 98 C for 5 min. The synthesized cDNA was quantitatively determined by Nanodrop 2000 (Thermo Fisher Scientic, West Palm Beach, Florida, USA). Real-time PCR was conducted on Bio-Rad iCycler iQ (Bio-Rad, Hercules, CA, USA) with THUNDERBIRD SYBR qPCR mix kit (Toyobo, Osaka, Japan). Gene ACT1 of K. marxianus NBRC1777 was used as an internal control. The primers for RT-PCR are shown in Table S4. † The cycle threshold values (C T ) were determined and the relative fold differences were calculated by the 2 ÀDDCT method 20 using ACT1 as the endogenous reference gene. Samples were run in triplicate in a 96-well plate, and each experiment was repeated three times. The value of qPCR presented is the mean of the triplicate results.

Measurement of the intracellular coenzyme contents
Intracellular NAD(P)H and NAD(P) + were extracted using EnzyChrom™ NAD(P) + /NAD(P)H assay kit (BioAssay Systems, Hayward, California, USA) following the manufacturer's instruction. A 10 ml sample of yeast culture was withdrawn and sprayed into quenching solution (60% methanol and 70 mM HEPES). Then, the quenched pellets were resuspended in acid extraction buffer or base extraction buffer (BioAssay Systems), aer which they were incubated at 65 C for 5 min in a water bath to extract oxidized pyridine nucleotides or reduced pyridine nucleotides, respectively. The relative amounts of NAD + , NADH, NADP + , and NADPH in each extract were then quantied by enzymatic methods using a NADP + -specic glucose dehydrogenase cycling reaction and a NAD + -specic lactate dehydrogenase cycling reaction, in which the formed NADPH or NADH reduces a formazan (MTT) reagent (BioAssay Systems). 21,22 At OD 600 ¼ 1, the concentration of the cells was equivalent to 0.411 g l À1 dry cell weight (DCW). 23

Overview of transcriptional data with mixed fermentation inhibitors by RNA-seq
Yeasts reacted differently with various inhibitors in pretreated hydrolysate. To mimic the real fermentation procedure, here we used the mixed three main group inhibitors in the lignocellulosic pretreated hydrolysate (furfural, acetic acid, phenols) and the compound concentrations were used as the previous report [24][25][26] with minor modication according to the growth of K. marxianus YHJ010 which derived from K. marxianus NBRC1777. 19 With the treatment of multiple inhibitors, the cells growth was slower than those without inhibitors, as shown in Fig. 1.
The alteration of genome-wide gene expression was analyzed by RNA-seq analysis of K. marxianus YHJ010 with or without multiple inhibitors treatment. A total of 13 622 794 and 17 697 296 clean reads were obtained from the RNA of K. marxianus with or without inhibitors treatment, respectively, and over 91% were uniquely mapped to the reference genome (Table S1 †).

Identication of differentially expressed genes and validation by qPCR
The levels of gene expression, normalized as reads per kilobase of exon model per million mapped reads (RPKM), were applied to the fold changes of the differentially expressed genes (DEGs, with absolute fold changes (FC) $ 2; p # 0.05). Under the stress with multiple inhibitors, 451 transcripts were identied to have different expression levels compared to those without stress. Among them, 279 genes were up-regulated and 172 genes were down-regulated (statistical data from the differentially expressed genes, data not shown). We then performed GO enrichment analysis on these DEGs. As shown in Table S2 and Fig. S1, † most DEGs in glycolytic process (GO:0006096), monocarboxylic acid metabolic process (GO:0032787), pyruvate metabolic process (GO:0006090) and NADP metabolic process (GO:0006739) etc. were down-regulated under the multiple inhibitors stress; on the other hand, most DEGs in tricarboxylic  (GO:0006099), transmembrane transport (GO:0055085), single-organism transport (GO:0044765), oxidative phosphorylation (GO:0006119), transporter activity (GO:0005215), respiratory chain (GO:0070469), substrate-specic transporter activity (GO:0022892), cellular response to chemical stimulus (GO:0070887) and response to stress (GO:0006950) etc. were up-regulated.
We also performed KEGG enrichment analysis on these DEGs. The ratio of DEGs/background genes indicates the effects of the DEGs in the specic KEGG pathway (background genes). As shown in Table S3 and Fig. S2, † the expression of most genes in degradation of aromatic compounds (Ko0120), retinol metabolism (Ko00830), drug metabolism (Ko00983), glycolysis/ gluconeogenesis (Ko00010), methane metabolism (Ko00680), carbon xation in photosynthetic organisms (Ko00710), citrate cycle (Ko00020) etc. were signicantly regulated under the inhibitors stress, suggesting that the yeasts boot up the detox-ication response to deal with the fermentation inhibitors in environment.
To conrm the reliability of data from RNA-seq, sixteen genes involved in various pathways were selected for a quantitative real-time PCR (qPCR) comparison of their expression levels. As illustrated in Table 1, although the relative expression levels (the fold changes, shown as the sign of the log 2 transformed fold change (FC) values, log 2 FC(I/C)) of each selected gene were different between RNA-seq and qPCR, the trends of up-or down-regulation of all genes chosen were the same, which consequently demonstrated the accuracy of the trends of gene expression change obtained by RNA-seq ( Table 1).

Characterization of DEGs involved in central carbon metabolism
In this study, most DEGs related to glycolysis/gluconeogenesis and by-products formation were down-regulated except FBP1, GUT2, TDH2, ADH3/4/6 and PCK1 (Table 2 and Fig. 2), which is consistent with relatively slow growth speed with inhibitors stress condition (Fig. 1). In the glycerol formation pathway, GUT2 encoding glycerol-3-phosphate dehydrogenase was upregulated with log 2 FC value of 3.61 under multiple inhibitors condition, while DAK1 coding for dihydroxyacetone kinase was down-regulated with log 2 FC value of À2.27( Table 2). Among the most dramatically regulated transcripts are those for glyceraldehyde-3-phosphate dehydrogenase (TDH1, 2 and 3), whereas TDH1 and TDH3 were among the most abundant transcripts of all three isoforms and heavily down-regulated with the inhibitors stress, TDH2 is strongly induced with the same stress. On the other hand, all DEGs related to TCA were up-regulated, and notably, quite a few genes are related to the generation of NADH, such as IDH1/2, KGD1, 2 and MDH2 (Fig. 2, Table 2).
It was noticeable that the isoforms of ADH were dramatically regulated with the multiple inhibitors stress condition. As shown in Table 2, the transcript of ADH2 (KMAR_40226), ADH4 (KMAR_20152) and ADH6 (KMAR_80326) was the three highest abundant isoforms. Under the multiple inhibitors stress, ADH4 and ADH6 were up-regulated with log 2 FC value of 3.18 and 3.62, respectively. ADH3, another ADH isoform which encodes an ethanol-acetaldehyde redox shuttle involved in the transfer of redox equivalents from the mitochondria to the cytosol, was upregulated with log 2 FC value of 2.23 ( Fig. 2 and Table 2), consistent with previous studies that ADH3-disrupted K. marxianus was more sensitive to the reactive oxygen species and the null mutant of ADH6 was hypersensitive to vanillin, a major phenolic aldehyde compounds derived from lignocellulosic biomass, in S. cerevisiae. 18,27 On the other hand, ADH2 was down-regulated with log 2 FC value of À6.73. Unlike other ADHs, ADH2 catalyzes the reaction of ethanol to acetaldehyde and is repressed in the presence of glucose, and the repressed expression in our study might be explained that the addition of    furfural inhibited the glucose consumption and led to higher glucose concentrations and this in turn repressed the expression of ADH2. 28

DEGs related to mitochondrial respiratory chain and ATPases
The mitochondrial respiratory chain, which forms membrane potentials to produce ATP by ATPase, consists of vital components to transfer electrons, and those electron carriers function in the form of multienzyme complexes. As shown in Table 2, under the stress of the multiple inhibitors, NDI1 coding rotenone-insensitive NADH-ubiquinone oxidoreductase, and all of the SDHs coding rotenone-insensitive NADH-ubiquinone oxidoreductases were dramatically up-regulated more than 8fold (log 2 FC > 3). And almost half genes of those coding cytochrome b-c1 complex units were up-regulated (Table 2). However, to those coding core subunits of cytochrome c oxidase, there was no quite difference with the stress of the mixed inhibitors (data not shown). Previous report showed that mutation in SDH increased ROS production to nuclear and mitochondrial genomic instability. 29 Our ndings indicated that to some extent SDH-bc1 complex-cytochrome c peroxidase played a role in scavenging ROS released from the tolerance to multiple inhibitors in K. marxianus. F 1 F 0 ATP synthase is a large, evolutionarily conserved enzyme complex required for ATP synthesis. Among vast genes encoding subunits of F 1 F 0 ATP synthase complex (F-type ATPase), only ATP1, ATP14 and ATP16 which encoding alpha subunit of the F 1 sector, subunit h of the F 0 sector, and delta subunit of the central stalk of mitochondrial F 1 F 0 ATP synthase, respectively, were up-regulated more than 4 fold (log 2 FC $ 2) under multiple inhibitors ( Table 2), suggesting that the energy production was important to the tolerance to inhibitors stress.
V-ATPase maintains the acidity of the vacuole and generates the electrogenic potential that is used to drive the accumulation of ions and small molecules, amino acids and metabolites. V-ATPase-depleting mutants exhibited sensitivity to the acids. 30 Interestingly, novel roles of V-ATPase in the regulation of cellular receptors and their trafficking via endocytotic and exocytotic pathways were recently uncovered. 31 Also, defects in acidication, through defects in the vacuolar H + -ATPase, will lead to defective assembly of the high affinity iron transport  Further details are given in Table 2.
system. 32 In this study, ATP6c coding v-type proton ATPase subunit c was down-regulated under multiple inhibitors stress, meanwhile, iron transporters coding genes such as ARN2, SIT1 were up-regulated (Table 2).

DEGs involved in genes related to ROS detoxication
Reactive oxygen species (ROS) are a group of molecules derived from molecular oxygen and have toxic effects that can damage a wide variety of cellular components resulting in lipid peroxidation, protein oxidation, and genetic damage through the modication of DNA. Inhibitors like acetic acid, furfural, and phenols have been reported to be related to the redox state inside cells, inducing reactive oxygen species (ROS) generation. [33][34][35] Genes related to ROS detoxication including those coding for superoxide dismutases (SODs) and their chaperones, catalases and peroxidases, glutathione and thioredoxin systems. These proteins remove excess ROS such as cOH, H 2 O 2 , and O 2 c À etc. by participating in oxidation-reduction reactions to maintain normal cell metabolism and to ensure a high rate of cell viabilities by their activated dimers. Interestingly, only several oxidative stress-response genes were found to be up-or down-regulated under the mixed inhibitors treatment ( Table 2). As shown in Table 2, SOD1 and SOD2, coding cytosol Cu/Zn superoxide dismutase and mitochondrial superoxide dismutase [Mn] respectively, were up-regulated under the stress of multiple inhibitors, suggesting that SODs played an important role in multiple inhibitors tolerance in K. marxianus. CTA1 and CTT1, corresponding to two isoforms of catalase with different sub-cellular locations, peroxisomal-mitochondrial matrices and cytosol, respectively, 36 while CTT1 was down-regulated signicantly to À5.52 of the log 2 FC under the stress of inhibitors (Table 2), CTA1 had no signicant change (data not shown). As to those genes encoding thioredoxin peroxidases, PRX1 was up-regulated with log 2 FC value of 3.21 while AHP1 was downregulated with log 2 FC value of À3 under the stress of mixed inhibitors (Table 2). It was reported that expression of TSA1 of K. marxianus which encoding peroxiredoxin TSA1 enhanced the tolerance to a mixture of formic acid, acetic acid and furfural in S. cerevisiae. 37 In this study, however, there was no signicant change with TSA1 (data not shown), indicating the different role of this TSA1 in K. marxianus from that in S. cerevisiae. The essential coenzymes nicotinamide adenine dinucleotides, NAD(P) + and NAD(P)H, participate in key redox reactions and contribute to maintaining cell tness and genome stability. 38 Those genes such as ADH3, ALD6, IDH1/2, GDH1 and NDI1 etc. coding for NAD(P)H/NAD(P) + shuttle systems which play a key role in the maintenance of the mitochondrial redox balance by redox transformation from NAD(P) + to NAD(P)H were up-regulated in our RNA-seq result ( Table 2).
The ratio between reduced and oxidized co-factors is thought to play a major role in metabolism since several enzymes are regulated by this ratio. 39 In the present study the NADH/NAD + and NADPH/NADP + ratio were used to determine the change of redox balance. As shown in Table 3 and Fig. 3, with 2 h multiple inhibitors treatment, the concentration of NAD + was dramatically increased from 394.64 nmol g À1 DCW to 887.63 nmol g À1 DCW, while the concentration of NADH was only a little less than that of no stress, leading to the ratio of NADH/NAD + decreased from 0.74 to 0.28 (Table 3 and Fig. 3). The concentration of NADH and NAD + pool was increased from 686.96 nmol g À1 DCW to 1132.14 nmol g À1 DCW. On the other hand, with the multiple inhibitors stress, the concentration of NADP + was increased from 37.88 nmol g À1 DCW to 58.28 nmol g À1 DCW, while the concentration of NADPH was decreased from 26.84 nmol g À1 DCW to 13.40 nmol g À1 DCW, leading to the ratio of NADH/NAD + decreased from 0.71 to 0.23 (Table 3 and Fig. 3.). The concentration of NADPH and NADP + pool was increased from 64.72 nmol g À1 DCW to 71.68 nmol g À1 DCW, not changed so much like NADH + NAD + pool ( Table 3). As a result, with the multiple inhibitors stress, the concentration of total coenzymes was dramatically increased from 751.69 nmol g À1 DCW to 1203.81 nmol g À1 DCW (Table 3). Consistently, SDT1 encoding suppressor of disruption of TFIIS which was reported to be responsible for production of precursors in NAD + synthesis in cells, 40 was up-regulated in our study ( Table 2).

MSN2/4 mediated STRE related DEGs
MSN2/4 regulated the expression of a wide variety of genes in response to multiple types of stress in S. cerevisiae. 41 Table 2 lists DEGs which are reported to be controlled via the stress responsive element (STRE) in S. cerevisiae. 41,42 As shown in Table  2, quite a few stress-responsive genes were down-regulated under the stress of multiple inhibitors, such as HXK1, GPH, TDH1, TDH3, PGM, CTT1 etc., while there were also some genes up-regulated such as SSA3, SOD2, TDH2, HSP26, HSP31, ALD, MDH2 etc. Our RNA-seq results showed that there was only MSN2 in K. marxianus and with no signicant change under the multiple inhibitors condition (data not shown), indicating that this MSN2 might be different to that in S. cerevisiae. Due to the few study of this protein in K. marxianus, the specic function of MSN2 needs to be further investigated.

DEGs related to fatty acid and ergosterol metabolism
Ergosterol and fatty acids are considered as two critical membrane components associated with tolerance to multiple stresses. 43,44 As shown in Table 2, there were 7 DEGs with related to fatty acid metabolic process, and 3 DEGs with related to the ergosterol biosynthetic process, and almost all of them were down-regulated except LipA (KMAR_50026) coding lipoyl synthase. Among those down-regulated genes, OLE1 encodes acyl-CoA desaturase 1 involved in desaturation of fatty acid and with greatest transcriptional abundance ( Table 2).
DEGs related to alanine, aspartate and glutamate metabolism, vitamin B1 and B6 metabolism In our study, AGX1 encoding alanine-glyoxylate aminotransferase 1 related to alanine synthesis, UGA1 encoding 4-aminobutyrate aminotransferase and GDH1 encoding NADP-specic glutamate dehydrogenase 2 which are related to glutamate synthesis, respectively, were up-regulated with the multiple inhibitors stress condition, while ADSS encoding adenylosuccinate synthetase and ASN1 encoding asparagine synthetase 1 related to aspartate synthesis were down-regulated with the same stress condition ( Table 2).
As to vitamin B6 and B1 metabolism, all the DEGs related to this category were down-regulated except one gene encoding a putative pyridoxal reductase with the multiple inhibitors stress (Table 2). Interestingly, though genes encoding probable pyridoxine biosynthesis protein SNZ3 and probable pyridoxal 5 0phosphate synthase SNO3 were dramatically down-regulated (Table 2), in our result, however, there were no SNZ1 and SNO1 corresponding to the counterparts of S. cerevisiae in K. marxianus, which suggests that SNZ3 and SNO3 of K. marxianus might be quite different from those in S. cerevisiae.

DEGs related to transcription factors
Transcription factors (TF) play an important role in regulatory mechanisms underlying various stresses resistance mechanisms. As shown in Table 2, DEGs related to TF are widely distributed in the regulatory network for diverse biological processes, including carbohydrate metabolism (GCR1, GCR2), sulfur amino acids metabolism (MET32), lipid metabolism (OAF1), cell proliferation (HCM1), and transcriptional process (YNG1, TFC7, MTF1, MED19, TFIIF2) etc. of those TFs with increased transcriptional expression under inhibitors stress condition, HCM1, as a transcription factor involved in cell cycle regulation, is also involved in promoting mitochondrial biogenesis and stress resistance in S. cerevisiae; 45 OAF1 was reported to regulate genes involved in beta-oxidation of fatty acids, peroxisome organization and biogenesis, activating transcription in the presence of oleate, 46 though the transcript abundance in this study is very low. As to the carbohydrate metabolism related TF genes, GCR1 and GCR2 encoding glycolytic genes transcriptional activator GCR1 and GCR2, respectively, were down-regulated with multiple inhibitors stress, and this was consistent with a previous report that GCR1 and GCR2 mutants showed lower glycolytic activities and enhanced the expression of TCA and respiratory genes to produce more energy. 47 We also noticed that several DEGs related to the transcriptional factors that are directly involved in the transcriptional process were regulated with the multiple inhibitors stress. TFIIF2 encoding transcription initiation factor IIF subunit beta, MTF1 encoding a mitochondrial transcriptional factor that confers selective promoter recognition on the core subunit of the yeast mitochondrial RNA polymerase, and YNG1 encoding a component of the NuA3 histone acetyltransferase complex that post-translationally modies histones, 48 were up-regulated, while TFC7 encoding a component of the initiation complex which functioned in RNA polymerase III recruitment and MED19 encoding a subunit of mediator were down-regulated with multiple inhibitors stress ( Table 2). Mediator binds transcription activation domains and Pol II, allowing activatordependent Pol II recruitment. 49,50 These results indicated that the processes of transcription initiation, transcription activation, the promoter recognition were selectively regulated by the multiple inhibitors stress.

DEGs related to transporters and the transcriptional response with individual inhibitor stress
Transporters play key roles in the response of the fermentation inhibitors. As an important component of transporters, major facilitator superfamily (MFS) is a large and diverse group of secondary transporters that includes uniporters, symporters, and antiporters. 51 As shown in Table 2 (Fig. 4). Previous report showed that myo-inositol improved tolerance of S. cerevisiae to the mixture of furfural, acetic acid and phenol, and deletion of gene in myo-inositol synthesis weakened strain tolerance against this stress. 52 Our result further indicated that ITR2 may play an important role against furfural in lignocellulose-derived inhibitors stress.
As to the amino acid transporters, TAT2 (KMAR_10514) encoding tryptophan permease and GAP1 encoding general amino-acid permease were down-regulated with log 2 FC value of À2.89 and À2.52, respectively, under the stress of mixed inhibitors. Tryptophan can be converted to quinolinic acid (QA), an important precursor in NAD + synthesis. 53 On the other hand, YCT1 encoding high affinity cysteine transporter and KMAR_40340 encoding cystine transporter were up-regulated about 8-fold than that with no stress condition ( Table 2).
Under the stress of multiple inhibitors, TNA1, encoding high-affinity nicotinic acid transporter which was essential for the NAD + homeostasis, 54 was down-regulated, while DAL5 encoding an allantoate and ureidosuccinate permease subjected to nitrogen catabolite repression 55,56 and SUL2 encoding sulfate permease 2 were up-regulated, but the transcript abundance was too low ( Table 2). In S. cerevisiae, it was observed that the genes involved in sulfur metabolism are mainly regulated by the cellular cysteine pool. 57 The ARN family encodes proteins involved in the uptake of siderophore-iron chelates. Genome-wide analysis showed that the acidic condition affects metal metabolism. 30 From our RNAseq results, ARN2 and SIT1 were signicantly up-regulated with log 2 FC value of 5.19 and 3.31, respectively, under the stress of mixed inhibitors (Table 2). FTR1 encoding plasma membrane iron permease was also up-regulated, though FET4 encoding low-affinity Fe 2+ transport protein was down-regulated. Meanwhile, a gene FSF1 encoding a probable mitochondrial transporter which was reported to be necessary to maintain the homeostasis of iron, 58 was down-regulated with multiple inhibitors condition (Table 2). ARN2 was found to be induced under the acid adaptation and acid affects metal metabolism. 30 In our individual inhibitor stress experiment, however, though ARN2 was induced under the acetic acid condition, the mostenhanced expression was with furfural stress, and so was that of SIT1 (Fig. 4), indicating that furfural may affect iron transportation more than acidic condition in K. marxianus. High affinity copper transporter coding gene, CTR1, was repressed with the multiple inhibitors stress (Table 2). Interestingly, low affinity copper uptake can be mediated by FET4, which was also low affinity iron transporter, 59,60 and the coding gene FET4 was also repressed in this study ( Table 2), indicating that the multiple inhibitors stress inhibited the copper uptake.
In our study, MEP3 encoding ammonium transporter MEP3 was up-regulated, while OPT1 encoding oligopeptide transporter 1 and CTP1 encoding a tricarboxylate transport protein were down-regulated with the multiple stress (Table 2).
Under the stress of multiple inhibitors, PET9 encoding a mitochondrial ADP, ATP carrier protein was dramatically upregulated, which was consistent with the up-regulation of those genes coding for ATP synthase (Table 2).
Efflux system of living cells is an efficient mechanism for detoxication of external toxic compounds and internal damaging intermediates. Two multidrug permease gene, ATR1 encoding aminotriazole resistance protein and FNX1 encoding multidrug resistance protein were up-regulated under the mixed fermentation inhibitors (Table 2). FNX1 also showed increased transcriptional expression under furfural or phenols  stress (Fig. 4). These results were consistent to previous reports that ATR1 deletion mutant S. cerevisiae showed increased sensitivity to lignocellulosic inhibitors and FNX1 mutant S. pombe presented impaired uptake of vacuolar amino acid. 61,62 In addition, KMAR_40425 encoding an uncharacterized polyamine transporter 4 was up-regulated under multiple inhibitors stress (Table 2), consistent to a recent report that higher spermidine was able to enhance tolerance of S. cerevisiae against lignocellulose-derived inhibitors. 33 Carboxylic acid transporter protein JEN1 was found to be involved in the acids efflux and the transport of the substrate is bidirectional. [63][64][65] Our RT-PCR results showed that JEN1 was signicantly up-regulated with log 2 FC value of 8.47, 6.85 and 3.84, under the furfural, acetic acid and phenols stress respectively, compared with no stress condition (Fig. 4). Meanwhile, another gene JEN2 encoding putative sialic acid transporter with 34.7% identity with JEN1 of S. cerevisiae and 74.3% identity with JEN2 of Kluyveromyces lactis, was up-regulated with log 2 FC value of 4.91, 6.28 and 3.07, under the furfural, acetic acid and phenols stress respectively, compared with no stress condition (Fig. 4). Both RNA-seq and RT-PCR results showed that these two genes were up-regulated under the multiple inhibitors stress ( Table 2 and Fig. 4). This suggests that JEN1 and JEN2 respond with different stress and play an important role against the mixed inhibitors stress.
PDR12, an ATP-binding cassette (ABC) transporter and a member of the Pleiotropic Drug Resistance (PDR) family, was demonstrated to be essential to the acquisition of tolerance to weak acid stress, being involved in the extrusion of the carboxylate anions and participating in cellular detoxication. 66 Our RT-PCR results also showed that PDR12 was induced with log 2 FC value of 5.03 in respond to acetic acid stress compared with no stress condition, the most up-regulated among three stress conditions (Fig. 4). This suggests that PDR12 be an interesting protein especially against acid stress. Another ABC transporter gene YCF1, encoding metal resistance protein YCF1 which was reported to function in the detoxication of furfural and/or HMF 67 and mediated transport of GSH-conjugated metals for metal tolerance, 68 was also up-regulated under mixed inhibitors stress (Table 2).

Discussion
Exploring the toxicity of lignocellulose-derived inhibitors to yeasts and developing the excellent strains with enhanced tolerance are becoming a more critical component of producing chemical products from lignocellulosic materials. Transcriptomic data obtained from inhibitors tolerance experiments is an extremely important step in the construction of large scale, system-based models that can be used to predict the cellular response to this stress. The integration of this data with pathway information is crucial to improve the accuracy in the prediction capabilities of the models. In the case of unconventional thermotolerant yeast K. marxianus, however, there were only several reports of transcriptome analysis on yeast responses from different carbon source including glucose, inulin and xylose, or the tolerance to high temperature or ethanol stress. 18,69,70,71 To the best of our knowledge, genome-wide transcription analysis under the multiple inhibitors stress condition with elevated temperature (42 C) has not been reported for this yeast.
Carbon central metabolism plays an important role in carbon source and energy production to yeast cells. From our results, differentially expressed genes related to the carbon central metabolism were selectively regulated by multiple inhibitors stress. Though previous report showed that the genes and proteins associated with glycolysis were over-expressed under acetic acid stress, 72,73 we noticed that DEGs related to the glycolysis were depressed in response to the multiple inhibitors ( Table 2, Fig. 2), while those related to TCA and a gluconeogenesis specic gene FBP1 were up-regulated, and consistently, most DEGs encoding the respiratory chain component functioned in the oxidative phosphorylation in mitochondria were up-regulated ( Table 2 and Fig. 2), together with the up-regulation of mitochondrial ADP/ATP carrier gene PET9, suggesting that inhibitors stimulate cells to produce more ATP. Cells need to choose the most efficient route to generate energy or reduce ATP consumption to maintain energy reserves under environmental stress condition. We speculate that cells choose to slow down the metabolic ux in glycolysis pathway while turn to enhance TCA cycle to obtain more ATP production and more NADH, since detoxication of furfural or phenolic compounds is an energy-consuming process. Coincidently, the carbohydrate metabolism related TF genes GCR1 and GCR2 were also down-regulated with multiple inhibitors stress ( Table 2). GCR1 and GCR2 mutants were reported to show lower glycolytic activities and enhanced the expression of TCA and respiratory genes to produce more energy, 47 in addition, overexpression of GCR1 increased transcription levels of HXT1 and ribosomal protein genes in S. cerevisiae. 74 Combined with our results, these studies indicated that in K. marxianus GCR1 and GCR2 may play a role with tolerance to the hydrolysates inhibitors by regulating carbon central metabolism process to produce more energy.
On the other hand, as a protective mechanism responding to environmental stress, glycerol played a key role in keeping high cell viabilities during ethanol fermentation. In accordance with this, up-regulation of GUT2 and down-regulation of DAK1 in favor of the glycerol formation pathway was observed ( Table 2).
As three main lignocellulose-derived inhibitors, acetic acid affects cell metabolism and stabilities of proteins by a drop in intracellular pH and membrane potential, furfural inhibits glycolytic and fermentative enzymes essential to central metabolic pathways, and phenolic compounds alter the permeability of biological membranes and caused irreversible damages to the cells. 75 All these inhibitors have been reported to be related to the redox state inside cells inducing ROS generation. [33][34][35] Furfural and HMF were reported to inhibit alcohol dehydrogenase (ADH), pyruvate dehydrogenase (PDH), aldehyde dehydrogenase (ALDH), hexokinase (HXK) and glyceraldehyde-3phosphate dehydrogenase (GPDH) in S. cerevisiae, 75 in our study, however, at least at the transcriptional level, only HXK, ADH2 and TDH1/3 were down-regulated, ADH3/4/6 and ALD6 were up-regulated, and there was no obvious change on genes coding for pyruvate dehydrogenase (Table 2 and Fig. 2).
Previous study reported that NADPH-dependent oxidoreductases comprise the main resistance mechanism for high concentrations of furfural. 76 Expression of some oxidoreductases could enhance the tolerance of cells to furfural, acetic acid and phenolic compounds in lignocellulosic hydrolysates, and intracellular ROS in cells with an increased tolerance has been reported to be decreased. 37,77,78 In our study, in response to multiple inhibitors stress, the transcripts for the genes encoding known NAD(P)H/NAD(P) + shuttle systems, including ADH3/ 4/6, ALD6, TDH2, GUT2, IDH1/2, GDH1 and NDI1 showed high levels of enhanced expressions, and transcripts for enzymes involved in the malate-oxaloacetate shuttle or malate-pyruvate shuttle, encoded by MDH2 and MAE1, were also induced (Fig. 2, Table 2). Previous report showed that MDH could be regarded as a transhydrogenase-like shunt, which regulated the redox state in S. cerevisiae. 79 ROS overproduction in response to the inhibitors is another reason for redox imbalance in yeast. ROS scavenging proteins remove excess ROS such as cOH, H 2 O 2 , and O 2 c À etc. generated from the multiple inhibitors by participating in oxidationreduction reactions and this requires the reducing power. Detoxication of lignocellulosic inhibitors like furfural, HMF or phenolic compounds is a process of converting them into less toxic corresponding alcohols in NAD(P)H-dependent reduction, 5,80 which requires the supply of sufficient amounts of the involved co-enzymes. This was consistent to the decrease of NAD(P)H/NAD(P) + ratio in our study and others report. 28 In this study, the total amount of NAD + and NADH increased nearly one fold when the strain exposed to inhibitors (Table 3), whereas the total amount of NAD + and NADH was decreased in the case of S. cerevisiae. 28 One possible reason of the increased intracellular concentration of NAD + might be that the NAD + synthesis was increased under the multiple stresses, based on the up-regulation of SDT1 in our study (Table 2), which was reported to be responsible for production of precursors in NAD + synthesis. 40 Combined with the up-regulation of those genes involved in NAD(P)H generation, such as ADH3/4/6, ALD6, TDH2, IDH1/2, GDH1, KGD1/2 etc. (Fig. 2, Table 2), indicating that more NAD(P)H production could be provided. This distinctive character of enhancing NAD level in response to the multiple inhibitors may endue K. marxianus intrinsic considerate inhibitors tolerance especially to furfural and HMF, which was reported in our and other previous study. 13,14 These results give us a hint that improving the amount of NAD + and NADH may enhance the yeast tolerance to lignocellulosic inhibitors.
In addition, YCT1 encoding high affinity cysteine transporter and KMAR_40340 encoding cystine transporter were upregulated with the multiple inhibitors stress. YCT1 was reported to be the principal cysteine transporter in S. cerevisiae. 81 It is well known that cysteine with reductive SH is required for the synthesis of glutathione, an essential antioxidant molecule involved in oxidative stress response and detoxication. 82 Combined with the enhanced expression of enzyme genes at NADH/NAD + shuttle sites in our study and the increased amount of NAD + and NADH pool (Table 3), it once again suggest that the regeneration or conserve actual cofactors was important to remain the cellular redox balance to K. marxianus under the lignocellulosic inhibitors stress.
We also noticed that DEGs related to alanine, aspartate and glutamate metabolism were signicantly regulated in response to the multiple inhibitors stress (Table 2). This may be explained by a previous report that the alanine, aspartate and glutamate metabolism was important for yeast cells to resist furfural, acetic acid and phenol (FAP) stress. 83 Though SNZ1 and SNO1 were required for conditions in which vitamin B6 (pyridoxal) is essential for growth, SNZ2/SNO2 and SNZ3/SNO3 pairs seemed more related with vitamin B1 (thiamine) biosynthesis during the exponential phase in S. cerevisiae. 84 In our RNA-seq results, however, there were only 2 genes encoded putative proteins showing close amino acid sequence similarity to SNZ3 and SNO3 of S. cerevisiae, and the transcript abundance of SNZ3 was very high, while both SNZ3 and SNO3 was dramatically down-regulated in response to multiple inhibitors stress ( Table 2), suggesting the encoded protein pairs may play an important role to the inhibitor tolerance in K. marxianus, though their precise functions in inhibitors tolerance remain to be elucidated. Furthermore, thiamine can affect metabolic functions through thiamine pyrophosphate (TPP)-dependent enzymes, such as pyruvate decarboxylase and alpha-ketoglutarate dehydrogenase which are important in the carbon central metabolism pathway and TCA cycle, respectively. In agreement with this, PDC coding for pyruvate decarboxylase and KGDs coding for alphaketoglutarate dehydrogenase were down or up-regulated with multiple inhibitors stress (Table 2). Meanwhile, it was reported that addition of thiamine decreased production of reactive oxygen species in yeast cells and decreases transcription of stress response genes as well. 85 Taken together with the only pair of SNZ3 and SNO3 and high transcript abundance in our study, different from those in S. cerevisiae, SNZ3/SNO3 may have multiple functions in K. marxianus.
MSN2 and its close homolog MSN4 (referred to as MSN2/4) were identied as transcriptional activator required for expression of a wide variety of genes in response to multiple types of stress via interaction with the consensus sequence known as the stress responsive element (STRE) in their promoter regions. 41 Overexpression of MSN2 of S. cerevisiae confers furfural resistance in S. cerevisiae and expression of MSN2 of K. marxianus promoted cell growth and ethanol production in S. cerevisiae. 86,87 We speculate that the function of MSN2 in K. marxianus might not relate to the inhibitors tolerance or phosphorylation of MSN2 was more important in regulating the genes in response to the multiple inhibitors stress.
Our study reveal that the fatty acid metabolic process and ergosterol biosynthetic process were depressed by multiple fermentation inhibitors, based on the 10 DEGs involved in these two biological processes (Table 2). Previous study also showed that overexpression of OLE1 improved the acetic acid tolerance in S. cerevisiae. 88 These results pointed a hint that regulating the expression of those DEGs involved in these two processes may increase the tolerance of K. marxianus to the lignocellulosic inhibitors.
Interestingly, a recent report showed that iron and copper are transition metals involved in redox reactions that are essential for all eukaryotes, but whose intracellular concentrations must be carefully monitored, as they are potentially toxic. 89 Our results showed that genes involved in iron homeostasis such as ARN2, SIT1 and FTR1were induced while those involved in copper uptake such as CTR1 and FET4 were repressed under multiple inhibitors conditions. In S. cerevisiae, most of these genes were regulated by AFT1, a transcription factor that responds to intracellular iron. 90,91 In our study, however, there was no change of AFT1 expression detected in transcriptional level (data not shown). Another gene KMAR_40422, encoding a probable mitochondrial transport protein FSF1 was repressed either. Interestingly, FSF1 was reported to belong to an ancient mitochondrial protein and necessary to maintain the homeostasis of iron within mitochondria. 58 Meanwhile, from our results, up-regulation of glutamate synthesis related genes UGA1 and GDH1 was consistent to previous report that regulation of glutamate synthesis was dependent on the iron availability. 92 Furthermore, the integrative analysis of the transcriptome with metabolome data revealed that the glucose metabolism, amino acid synthesis, ergosterol, and lipid biosynthesis biological processes were all affected due to the loss in the activities of specic iron-dependent enzymes under iron deprivation, 92 and the change of expression in transcriptional level were also identied in our study (Table 2). There are several mechanisms reported on iron uptake and the regulation on overall iron homeostasis is complicated. The results in present study give us a hint that there is relationship between iron transportation and the inhibitors tolerance though the mechanism remains unclear.
Besides the multiple inhibitors stress condition to mimic the lignocellulosic biomass fermentation to study the global transcriptional response of K. marxianus, we also investigated some transporters transcriptional response to the three individual inhibitors stress by RT-PCR. As predicted, these genes responded quite differently to different inhibitor. For example, ITR2 and JEN1, encoding myo-inositol transporter 2 and carboxylic acid transporter protein respectively, were dramatically upregulated especially with furfural stress, even more than that with the multiple stress condition (Fig. 4). Like MSN2 and SNZ3/ SNO3, ITR2 is another example of the only one isoform in K. marxianus in comparison to the 2 counterparts in S. cerevisiae. A previous study showed the essential role of ITR2 for Shizosaccharomyces pombe growth. 93 However, the regulation role of ITR2 in response to various stresses was not clear. We speculate that overexpression of ITR2 might enhance yeast to the lignocellulosic derived inhibitors tolerance. For PDR12, with acetic acid and mixed inhibitors, it seemed to have similar upregulation on the transcriptional expression level. In the case of ARN2, however, the most up-regulated stress condition was with mixed inhibitors treatment (Fig. 4). It should be noted that it is intrinsically complex and challenging to engineering yeast resistance to mixed fermentation inhibitors because each type of inhibitor may have distinct toxic effects and cellular stress response mechanisms. 94,95 Numerous metabolic pathways and regulatory genes have been reported affecting yeast tolerance to environmental stress. 96 It should be noted that some differentially expressed genes from RNA-seq dataset could be just passively up-or downregulated and may not contribute to eliciting stress responses. The future work will systemically evaluate the highly ranked differentially expressed genes and identify their effects on yeast stress responses to individual inhibitor in addition to mixed fermentation inhibitors. It is known that engineering microbial resistance to fermentation inhibitors becomes even more challenging and complex as the types of inhibitors expanded in the mixture. With the transcriptomic-guided metabolic engineering approach, our future work will concentrate on characterizing the highly ranked targets functions to elicit improved resistance to multiple inhibitors.

Conclusion
This study provides the rst transcriptomic analysis of K. marxianus at elevated temperature (42 C) to the multiple fermentation inhibitors stress. The results revealed that lignocellulosic hydrolysates inhibitors had effects on multiple aspects of cellular metabolism at the transcriptional level. Differentially expressed genes were enriched in maintaining the redox balance, NAD(P) + /NAD(P)H homeostasis or NAD + synthesis, energy production, and iron transportation or metabolism. Our results suggest that redox balance and NAD(P) + /NAD(P)H homeostasis play an important role in tolerance to lignocellulosic derived inhibitors. Besides engineering of the redox system to relieve the stress caused by HMF and furfural, improving the level of NADH and NAD + pool may represent another putative target to reduce the stress caused by the lignocellulosic derived inhibitors. Based on quite a few iron transporters were up-regulated in response to the multiple inhibitors, it hints that there is relationship between iron transportation or metabolism and the tolerance to the inhibitors though the mechanism remains unclear. The energyconsuming detoxication of furfural and HMF drove yeast to up-regulate those genes related to TCA cycle and respiratory chain to obtain more ATP. Some highly ranked differentially expressed genes were predicted as potential regulatory genes, which need further investigation. The results obtained in this study provide insights about the mechanisms of K. marxianus that are involved with lignocellulosic derived inhibitors response, enabling future metabolic engineering approach to obtain more robust strains for industrial fermentation with cellulosic biomass.

Conflicts of interest
The authors declare that they have no conict interests.