Transcriptomic Analysis of the Influence of Methanol Assimilation on the Gene Expression in the Recombinant Pichia pastoris Producing Hirudin Variant 3

Hirudin and its variants, as strong inhibitors against thrombin, are present in the saliva of leeches and are recognized as potent anticoagulants. However, their yield is far from the clinical requirement up to now. In this study, the production of hirudin variant 3 (HV3) was successfully realized by cultivating the recombinant Pichia pastoris GS115/pPIC9K-hv3 under the regulation of the promoter of AOX1 encoding alcohol oxidase (AOX). The antithrombin activity in the fermentation broth reached the maximum value of 5000 ATU/mL. To explore an effective strategy for improving HV3 production in the future, we investigated the influence of methanol assimilation on the general gene expression in this recombinant by transcriptomic study. The results showed that methanol was partially oxidized into CO2, and the rest was converted into glycerone-P which subsequently entered into central carbon metabolism, energy metabolism, and amino acid biosynthesis. However, the later metabolic processes were almost all down-regulated. Therefore, we propose that the up-regulated central carbon metabolism, energy, and amino acid metabolism should be beneficial for methanol assimilation, which would accordingly improve the production of HV3.


Introduction
Leeches (or Hirudinea) which belong to the phylum Annelida are closely related to oligochaetes (e.g., the best-known earthworm), and all have muscular, soft, and segmented bodies [1]. More than 700 species of leech living in marine or freshwater environments have been described [2,3]. Leeches are well-known to us for its blood-feeding animal features due to the secretion of hemotoxic proteins from the salivary glands into the blood streams after biting the host [4][5][6]. Leeches can usually consume over 10 times their weight in blood during feeding periods because of the successful counteraction against host hemostatic processes by those hemotoxic proteins as anticoagulants [7,8]. Given the powerful function in preventing the clotting of blood, the leech anticoagulants have been widely used in clinics as modern medicine, such as for venous congestion treatment and as post-surgical instruments following reconstructive and plastic surgery [9][10][11]. Among those anticoagulants from leech saliva, hirudin and its variants are the most famous and well-studied for its high efficiency in inhibiting thrombin with K I and K D values being in the pico-to femtomolar range [12].
Hirudin, as one of the most potent naturally-occurring anticoagulants, was firstly isolated from the saliva of leech Hirudo medicinalis [13]. It can efficiently and specifically bind to thrombin, which converts fibrinogen to fibrin in blood clot formation. Shortly after the release of the amino acid sequence of hirudin, three main subtypes of hirudins were discovered in succession and designated as hirudin variant 1 (HV1), HV2, and HV3, sharing a similar core domain of N-terminal sequence and inhibitory activity against thrombin [14][15][16][17].
Given the huge therapeutic value in clinical application and the limited availability of natural hirudin from leech, many researchers have put effort into improving the productivity and production of hirudin variants by genetic recombination with host-microorganisms, including Escherichia coli [18], Saccharomyces cerevisiae [19], and Pichia pastoris [20]. Among them, HV3 and its mutein have been successfully expressed in E. coli [21,22]. However, lack of a mature secretory system is still a critical defect of E. coli for heterologous protein secretion although additional signal peptide can direct HV3 to secrete into the medium [23]. Therefore, an eukaryotic host should be more appropriate for the production of HV3 due to the existence of post-translational modifications and a secretory system.
The methylotrophic yeast P. pastoris has been widely applied in the heterologous protein expression owing to its many advantages, such as achieving high cell density, post-translational modifications and a secretory system [24]. Three intramolecular disulfide bridges located in the hirudin N-terminal core domain are vital for its antithrombin activity through binding to the active center of thrombin to block its hydrolysis activity [25]. Therefore, the presence of post-translational processing system in P. pastoris is much more favored for the expression of hirudin and its variants.
Few reports covering the production of recombinant HV3 using P. pastoris have been published up to now. Although P. pastoris is a well-developed heterologous protein expression host, some disadvantages about this system are still noteworthy. The clonal variation is likely to be a major bottleneck for the industrial application of this host system as hundreds of clones per phenotype need to be screened to get the highest secretor [26,27]. Previous studies have usually employed two carbon sources, and adopted a two-stage cultivation strategy when culturing P. pastoris: Glycerol for biomass production and methanol for foreign protein expression. However, it was difficult for carbon source switching from one to the other at a precise or proper time point [24,28]. Furthermore, a deep and comprehensive insight on the influence of methanol assimilation on the overall metabolism in the recombinant P. pastoris is still not yet clear and is necessary to be explored further.
In this study, the recombinant P. pastoris GS115/pPIC9K-hv3 was successfully constructed, and the expression of hv3 regulated by promoter of AOX1 encoding alcohol oxidase (AOX) was realized. However, the growth of this recombinant strain slowed down after the initiation of methanol induction. Moreover, our work indicated that the cell growth should be closely related to the production of HV3 in the recombinant P. pastoris during the methanol induction phase. Thus, it is quite necessary to investigate the influence of methanol assimilation on the gene expression in general for further developing an effective strategy regarding the enhanced HV3 production by this strain. Hence, our research would provide a new approach for finding out or solving the bottlenecks for enhanced heterologous protein production in P. pastoris.

Materials and Reagents
The complete CDS of HV3 (Accession No. KR066905) was obtained from the NCBI website [29]. The gene hv3 was synthesized by referring to this reported sequence and the codon usage bias in P. pastoris, and this optimized sequence was deposited in GenBank (Accession No. MK341489). E. coli DH5α and P. pastoris GS115 were applied in the plasmid propagation, gene cloning, and the construction of recombinant for hv3 expression. Empty plasmid pPIC9K (Invitrogen, San Diego, CA, USA), which contains four restriction enzyme sites (EcoRI, SnaBI, AvrII, and NotI), was used for the plasmid engineering. E. coli DH5α/pUC57-hv3 obtaining the gene of hv3 was purchased from the Beijing Genomics Institute (BGI, Beijing, China). The analytical reagents used in the culture and fermentation experiments were purchased from local suppliers.

Fermentation of P. pastoris GS115/pPIC9K-hv3
We performed the cultivation of P. pastoris following a reported method and made some minor modifications [31]. Briefly, P. pastoris GS115/pPIC9K-hv3 and P. pastoris GS115/pPIC9K were precultured in YPD medium at 30 • C for 24 h, respectively. Afterwards, 2.5 mL of the above culture was pipetted into 50 mL of BMGY medium in a 250 mL Erlenmeyer flask and cultivated at 30 • C for 24 h with shaking at 250 rpm. The cells harvested by centrifugation at 4000× g for 5 min were resuspended in BMMY medium to achieve an OD 600 of 1.0, and were continuously cultured at 30 • C for 5 days to express HV3, with a supplementation of methanol to a final concentration of 1% (v/v) every 24 h. Then the broth of P. pastoris GS115/pPIC9K-hv3 and P. pastoris GS115/pPIC9K were collected for testing the antithrombin activity.
The fermentation of P. pastoris GS115/pPIC9K-hv3 was conducted in a 5 L fermenter (Baoxing, Shanghai, China). The initial fermentation procedure was as follows: 300 mL of P. pastoris GS115/pPIC9K-hv3 culture grown in BMGY medium for 24 h was inoculated into the fermentation tank containing 2.7 L of minimal salts medium: H 3 PO 4 (85%), 13 mL/L; CaSO 4 2H 2 O, 0.93 g/L; K 2 SO 4 , 18.2 g/L; MgSO 4 , 7.27 g/L; KOH, 10.6 g/L; sodium citrate 2H 2 O, 1.47 g/L; and 2 mL/L PTM1 trace metals solution [32]. Glycerol was used as the sole carbon source and its final concentration was 4% (w/v). In the initial phase, the fermentation parameters were set as follows: Temperature 30 • C, 30% of dissolved oxygen maintained by automatically controlling the agitation speed, airflow at 4 slpm, pH at 5 kept by automatically adjusting with ammonium hydroxide. The yeast cells were grown in this fermenter for about 18 h until the exhaustion of glycerol, observing a dramatic increase of dissolved oxygen (DO) and slowdown of agitation speed. Afterwards, 250 mL of 50% glycerol was gradually fed into the fermenter with a speed of 50 mL/h, and the pH was maintained at 3. The methanol induction phase was initiated when a sharp increase of DO was observed, a sign of glycerol depletion. And the methanol induction conditions were as follows: The methanol-fed rate was increased from 1 to 6 mL/h in the first 6 h. Afterwards it was set as: 6-16 h, 0.09 mL/min; 16-40 h, 0.18 mL/min; 40-90 h, 0.26 mL/min; 90-136 h, 0.18 mL/min. 270 mL of minimal salts medium was supplemented at 60 and 96 h, respectively.

Antithrombin Activity Analysis
The antithrombin activity of HV3 in the broth of P. pastoris GS115/pPIC9K-hv3 was analyzed following a previous method [33,34]. Briefly, the activity of HV3 was quantitatively measured by titrating a solution of thrombin and expressed as antithrombin units (ATU). One ATU of test sample is defined as the activity to neutralize one NIH unit of thrombin (Sigma, Darmstadt, Germany). The activity analysis was performed as follows: 0.2 mL of 0.05% bovine fibrinogen solution prepared with 50 mM Tris-HCl buffer (pH 7.4) was thoroughly mixed with 0.01-0.1 mL of sample, and then 5 µL of 100 NIH/mL thrombin solution (0.5 NIH unit) was added progressively and mixed gently. The above mixture was incubated at 37 • C for 1 min, and the end point of the titration was considered to be reached when a fibrin clot formed within 1 min. Otherwise, another 5 µL of thrombin solution was added continuously until a fibrin clot was observed.

cDNA Library Construction and Illumina Sequencing
Three biological replicates were set for each sample. Total RNA was isolated and purified using the RNeasy Mini Kit (QIAGEN, Alameda, CA, USA) and the RNA integrity of each sample was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). And a total of 3 µg RNA per sample was used for the cDNA library preparation. Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The cDNA fragments of preferentially 150-200 bp in length were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Lastly, the PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent).
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) following the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform in Novogene (Beijing, China) and 125 bp/150 bp paired-end reads were generated.

De Novo Assembly and Functional Annotation Analysis of Illumina Sequencing
Raw reads were firstly processed through in-house perl scripts. The transcriptome raw data were deposited in SRA database (Accession No. PRJNA509983). Clean reads obtained by removing reads containing adapter, ploy-N, and low-quality reads from raw reads were used for de novo assembly. Reference genome and gene annotation were downloaded directly from the database (NCBI, Komagataella phaffii GS115, assembly ASM2700v1). Index of the reference genome was built using Bowtie v2.2.3 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12. HTSeq v0.6.1 was used to count the reads numbers mapped to each gene, and the expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) of each gene was calculated based on the length of the gene and the count of reads mapped to this gene [35]. Prior to differential gene expression analysis, the read counts for the sequenced library were adjusted by the edgeR program package through one scaling normalized factor [36]. Differential expression analyses of comparison groups (three biological replicates per condition) were performed using the DEGSeq R package v. 1.18.0 [37,38]. The p-value was adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate [39]. An adjusted p-value (p adj ) of 0.01 and log 2 (fold change) of 1 were set as the threshold for screening the significantly expressed genes by DESeq. Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) was implemented by the GOseq R package, in which gene length bias was corrected [40].

HV3 Fermentation by P. pastoris GS115/pPIC9K-hv3 in 5L Fermenter
The sequencing result indicated the successful recombination of gene hv3 in P. pastoris GS115/pPIC9K-hv3 (data not shown). To further validate the successful translation of hv3 and its antithrombin activity, we cultivated P. pastoris GS115/pPIC9K-hv3 and P. pastoris GS115/pPIC9K in BMMY medium in which methanol served as the sole carbon source and an inducer for hv3 expression ( Figure 1A). The antithrombin activity in the broth of P. pastoris GS115/pPIC9K-hv3 reached the maximum value of 25 ATU/mL, while the fermentation liquor of P. pastoris GS115/pPIC9K did not show any antithrombin activity. Fermentation of P. pastoris GS115/pPIC9K-hv3 was conducted in a 5 L fermenter to further improve the production of HV3, and the highest antithrombin activity reached 5000 ATU/mL after 143.4 h of cultivation. Additionally, the synthesis of HV3 was almost synchronous with cell growth during the methanol induction phase ( Figure 1B). In order to obtain a deep insight on the transcriptional change after carbon source switching from glycerol to methanol, a comparative transcriptomic study of three timepoint samples (collected at 24 h, before the methanol induction initiation, designated as BI; at 96 h and 132 h, named as PI1 and PI2, respectively) was performed ( Figure 1B).

Overall Evaluation of RNA-Seq Data
The total RNA of sample BI, PI1, and PI2 were extracted and evaluated in triplicates, named as BI_1, BI_2, and BI_3; PI1_1, PI1_2, and PI1_3; PI2_1, PI2_2, and PI2_3, respectively. Then equal amounts of total RNA of each sample was employed for Illumina sequencing, and the overall quality of the high-throughput sequencing data was evaluated and shown in Table S1. The proportion of clean reads for each sample was approximately above 94%. Phred quality score (Q score) was employed to measure the base calling accuracy during sequencing [41]. Q 30 is equivalent to the probability of an error base call 1 in 1000 times (Table S1), meaning the base calling accuracy is 99.9%. Similarly, Q 20 stands for the 99% accuracy of base calling. The proportion of Q 20 was larger than 95%, and the proportion of Q 30 was no less than 88.92%. Therefore, the quality of sequencing data was good enough for the following sequence analysis. The percentage of reads uniquely mapped to the reference genome was distributed in the range from 92.28% to 94.42% (Table S2), indicating the validity of the employed reference genome. By comparing the gene expression level of sample BI with that of PI1 and PI2, we could conclude that the proportion of high abundant mRNA was increased in the recombinant P. pastoris during the methanol induction phase (Table S3). In order to check the sample quality of biological replicates for each sample, the Pearson correlation coefficients R 2 were calculated. The results showed that the R 2 of replicates for each group was above 0.997, indicating no significant difference between the biological replicates of each sample ( Figure S1).

Overview of Transcriptomic Analysis
1239 differently expressed genes (DEGs), including 690 up-regulated and 549 down-regulated genes, were determined in PI1 in comparison to BI ( Figure S2A). It indicated the significant change of gene expression pattern in P. pastoris after the carbon source switching from glycerol to methanol. When being continuously cultivated for 47 h, the number of DEGs in PI2 became 1365-696 up-regulated genes, 669 down-regulated ones ( Figure S2B). Compared with PI1, 14 genes were up-regulated in PI2, while the down-regulated genes were 107 ( Figure S2C). DEGs was clustered, and main gene expression pattern of sample BI was almost completely contrary with that in PI1 and PI2 ( Figure S3). The gene clusters with the low expression level (blue area) in BI almost entirely became up-regulated (red part) in PI1 and PI2.

GO (Gene Ontology) Enrichment Analysis of DEGs
The functional distributions of those DEGs of three groups (group A: PI1 versus BI; B: PI2 versus BI; C: PI2 versus PI1) were annotated and analyzed. The abundant genes were categorized into 30 major functional terms according to the GO category, and metabolic process, catalytic activity, single-organism metabolic process, oxidoreductase activity, and small molecule metabolic process were the top five in group A ( Figure S4A). DEGs belonging to the category of catalytic activity, single-organism metabolic process, and oxidoreductase activity were the most abundant parts in group B ( Figure S4B). However, the most abundant categories in group C were membrane and related parts belonging to cellular components, indicating cellular structure alteration as the proceeding of methanol induction ( Figure S4C).

Carbon Source Switching from Glycerol to Methanol Resulted in a Significant Slowdown of Energy Metabolism
Nearly all the detected genes encoding proteins distributed in almost all units of oxidative phosphorylation process were down-regulated in PI1 and PI2 in comparison to BI, except for few genes whose expression level didn't change significantly in PI2, such as PAS_chr3_1188 (encoding Ndufa9), PAS_chr3_0808 (Ndufs2), PAS_chr4_0535 (Ndufa6), PAS_chr2-1_0850 (ISP), PAS_chr1-4_0313 (QCR8), PAS_chr4_0520 (QCR9), and PAS_chr3_0576 (α) (Figure 3 and Table 2). PAS_chr1-3_0070 encoding mitochondrial inorganic pyrophosphatase, which converts diphosphate to phosphate for ATP production at the catalysis of ATPase, was also down-regulated in PI1 and PI2. Our results indicated that the energy metabolism in the recombinant P. pastoris declined comprehensively. Moreover, the expression level of genes involved in energy metabolism in PI2 did not show significant difference with that in PI1, according to our screening criteria (Table S5).

Expression Analysis of Genes Involved in the Biogenesis of Peroxisome and Peroxisomal Proteins
Nearly all the detected genes involved in peroxisome biogenesis were up-regulated in PI1 and PI2, except for PAS_chr4_0759 (encoding PEX12) whose expression level did not change significantly in PI1 in comparison to BI (Figure 4, Table 3). Some genes encoding peroxisomal proteins were also up-regulated in PI1 and PI2 in comparison to BI. PAS_chr2-2_0131 and PAS_chr4_0786 encoding the important and representative antioxidant enzymes, such as catalase (CAT) and superoxide dismutase (SOD), were remarkably up-regulated in both PI1 and PI2. On the contrary, PAS_chr1-4_0074 encoding outer mitochondrial carnitine acetyltransferase (CRAT), an enzyme involved in the fatty acid oxidation process in peroxisome, was down-regulated. But the other three genes (PAS_chr1-4_0538, PAS_chr4_0352, and PAS_chr2-1_0785) encoding fatty-acyl coenzyme A oxidase (ACOX) and long chain fatty acyl-CoA synthetase (ACSL) were up-regulated in this process. Besides, PAS_chr2-2_0272 encoding peroxisomal ATP-binding cassette transporter complex was also up-regulated in PI2, not in PI1. Genes encoding mevalonate kinase (MVK), isocitrate dehydrogenase (IDH), xanthine dehydrogenase (XDH) involved in sterol precursor biosynthesis, amino acid metabolism, and the purine metabolism process, respectively, were all down-regulated in both PI1and PI2. We also found that the expression of gene PAS_chr2-2_0207 encoding PEX13, an integral peroxisomal membrane protein, in PI2 increased significantly in comparison to PI1 according to our screening criteria (Table S5).

Transcriptomic Analysis of Genes Participating in the Protein Production and Degradation Processes
The expression level analysis of genes regarding to the protein production and degradation processes is vital and necessary for getting a deep understanding of the recombinant protein/peptide production in P. pastoris. We focused on four closely related aspects of the protein production and degradation, including basal transcription factors, protein processing in endoplasmic reticulum and ER-associated degradation (ERAD), SNARE interaction in vesicular transport, and ubiquitin mediated proteolysis ( Figure 5). Only a few DEGs were detected in both PI1 and PI2, in comparison to BI, during methanol induction stage (Table 4, Figure 5). The expression level of PAS_chr4_0204 and PAS_chr4_0745, which encode transcription factors MAT1 and TFIIH1, varied in contrast in PI1 and PI2 in comparison to BI-the former one increased, but the latter one decreased. As regards to SNARE interaction in vesicular transport in Golgi body, only PAS_chr1-4_0294 encoding Stx1-4 was down-regulated in both PI1 and PI2. PAS_chr1-4_0462 (Table 4) was solely down-regulated in PI1, not in PI2. PAS_chr2-2_0066 encoding Hsp40 was the only one with an increased expression level, but the others were all down-regulated. The expression level of genes involved in ubiquitin mediated proteolysis also varied in PI1 and PI2 in comparison to BI. PAS_chr3_0044 and PAS_chr1-3_0148, which encode E3 ubiquitin ligase (ARF-BP1) and ubiquitin conjugating enzyme (Apc3) respectively, were remarkably up-regulated in both PI1 and PI2. However, PAS_chr3_0924 and PAS_chr4_0429 encoding ubiquitin-conjugating enzyme (UBE2G2 and HIP2) was down-regulated. Besides, PAS_chr1-4_0609 encoding cullin, a structural protein of SCF complexes, a multi-protein E3 ubiquitin ligase complex, was also up-regulated in PI2, not in PI1. The expression level of genes involved in protein production and degradation processes in PI2 did not show a significant difference with that of PI1 (Table S5). A dash line indicates no significant change in comparison to BI.

Discussion
The methylotrophic yeast P. pastoris (syn. Komagataella spp.) has become one of the most important and widely-used expression systems for producing heterologous proteins, especially after the release of its genome sequence in 2009 [42]. High secretory capacity, a strong methanol inducible promoter and the existence of a glycosylation pathway are three attractive advantages over other eukaryotic expression systems, which allows a broad spectrum of recombinant proteins to be produced in this host, especially for those from higher eukaryotic cells [24,[43][44][45]. However, some problems of this system should not be ignored. The production of heterologously expressed protein is improved by increasing the gene copy number. But an optimal gene copy number seems to exist for each recombinant protein [46,47]. The clonal variation still seems to be the major bottleneck for the industrial application of this host system because hundreds of clones per phenotype need to be screened to get the highest secretor [26,27]. The productivity would be varied even if the recombinants own the same gene copy number. However, this problem cannot be attributed to only one level (transcriptional or translational level) [26]. We have often encountered similar problems, and the same operation of this host for different kinds of genes has usually brought in different results. The chance of successful expression and a high production of heterologous proteins have usually been very low. Moreover, the employment of two carbon sources (glucose or glycerol for cell growth and methanol for protein expression) and switching from one to the other at a precise or proper time point was difficult to realized [24,28]. Fortunately, we have the recombinant P. pastoris GS115/pPIC9K-hv3 which can produce HV3 efficiently, and have also found that the time for switching methanol induction has a profound influence on the final yield of HV3-an early or late time was not appropriate (data not shown). Moreover, we found that the function of methanol was not only in inducing exogenous genes expression, but in supporting cell growth as the carbon source. In a word, we found that the biosynthesis of HV3 was actually closely related to the growth of P. pastoris. If the cells hardly grew during the late methanol induction stage ( Figure 1B), the biosynthesis of HV3 also stopped. Shan Zhou et al. found that the specific growth rate and methanol concentration also influenced the degradation rate of hirudin produced by P. pastoris. The maximum specific hirudin production rate was achieved when a specific growth rate was maintained at 0.02 h −1 [48]. In order to get a full insight on the effects of methanol on the cell growth and expression of heterologous proteins in recombinant P. pastoris, and develop an appropriate strategy to further improve the production of HV3, we investigated the whole cellular gene expression of samples collected before and after methanol induction by transcriptomic method.
AOX encoded by gene AOX1 is the first and most important enzyme in initializing the methanol utilization pathway in methylotrophic yeast P. pastoris. Given that the promoter of AOX1 is strong-methanol inducible, the expression vector for heterologous protein is usually constructed by using the promoter sequence, and the induction of transcription is easily realized by using methanol as an inducer [24]. We adopted the same strategy for producing HV3 in the recombinant P. pastoris.
The expression of PAS_chr4_0821, encoding alcohol oxidase (EC 1.1.3.13), was greatly up-regulated in P. pastoris GS115/pPIC9K-hv3 during the methanol induction phase (Table 1, Figure 2, PI1 vs. BI, PI2 vs. BI), the log 2 (fold change) reached 6.5268 and 6.6696 in sample PI1 and PI2 in comparison to BI, respectively. Kim et al. investigated the regulation of AOX1 promoter activity and peroxisome biogenesis by visualizing the expression and localization of green fluorescent protein (GFP)-fused proteins under different fermentation processes, including a methanol-limited condition at normoxy (ML), switched feeding of carbon sources (glucose and methanol) under carbon-limited condition at normoxy (SML), and an oxygen-limited (OL) condition. The experiment results indicated that the yield of a foreign protein was highly dependent on the methanol consumption rate influenced by the availability of methanol and oxygen molecules [49]. But we could not agree with this viewpoint completely. Excessive methanol in the broth is harmful to the yeast cell [50]. Afterall, the production of heterologous protein in P. pastoris is a quite complex process, which is related to transcription, translation, protein transport and export, protein modification, excretion processes, and so on. Once the normal physiological state is affected by excessive toxic methanol, the productivity of protein would probably decrease in cells. Therefore, the key point for improving methanol consumption rate should be placed on the enhancement of methanol assimilation ability by P. pastoris, not on the strategy for the fed-batch methanol supplementation.
How can the methanol assimilation rate be fundamentally improved? To address this question, it is necessary to view the overall expression variation of genes involved in the central carbon metabolism in P. pastoris during the methanol induction phase ( Figure 2). AOX (EC1.1.3.13, Figure 2), dihydroxyacetone synthase (EC2.2.1.3, Figure 2, Table 1), and catalase (CAT, Table 3, Figure 4) were the three essential enzymes responsible for the conversion of methanol to non-toxic intermediates which subsequently entered into the central carbon metabolism pathway and supported the biomass production [51,52]. Kim et al. investigated the effect of methanol in inducing the biogenesis of peroxisome in the recombinant P. pastoris, which showed a high dependency on methanol consumption [49]. The transcriptomic results of our study also supported this viewpoint, and the detected DEGs involved in peroxisome biogenesis process were all up-regulated during the methanol induction phase (Figure 4). But the consumption of methanol should not merely depend on the above three enzymes present in peroxisome ( Figure 2). Actually, three branched pathways for methanol metabolism exist in yeast P. pastoris. In one branched pathway, methanol was finally converted into CO 2 with the catalysis of a series of enzymes, including EC1.1.1.284, EC3.1.2.12, and EC1.2.1.2, whose encoding genes were significantly up-regulated in the recombinant P. pastoris during the methanol induction phase. Glycerone-P and glyceraldehyde-3P which were converted from formaldehyde subsequently entered into the glycolysis/gluconeogenesis pathway and TCA cycle. However, the genes involved in these pathways, especially those encoding rate-limiting enzymes, were nearly all remarkably down-regulated. Obviously, this was not favorable for biomass production as the intermediates in these pathways were the main materials for protein, lipid, and polysaccharide biosyntheses. These results also coincided with our fermentation study which showed that the cells obviously grew slowly during the methanol induction phase ( Figure 1B). Besides, the expression level of those genes did not increase significantly with the induction time (PI2 vs. PI1, Table S5). Therefore, if we want to improve the methanol consumption rate and the productivity of foreign proteins, the overall enhancement of the central carbon metabolism, including glycolysis/gluconeogenesis, TCA cycle, should be feasible and effective. Joel Jordà et al. tried using a mixed carbon source (methanol and glucose) in the induction phase, and found a significant redistribution of carbon fluxes in the central carbon metabolism, reflecting not only in increased glycolytic, TCA cycle, and NADH regeneration fluxes but also as higher methanol dissimilation rates [53]. Hence this study also backed our proposal that increased carbon fluxes from methanol to glycolysis, TCA cycle and so on would not only benefit the growth of the recombinant P. pastoris, but also facilitated the methanol consumption and foreign protein production. Thus, a balanced central carbon metabolism should be focused on for the production of heterologous protein in P. pastoris.
Research focusing on the energy metabolism comparison in P. pastoris cultivated using glycerol or methanol as a carbon source are still very few. Joel Jordà et al. asserted that the oxidation of methanol requires large amounts of ATP, while Veeresh et al. held that the oxidation process from GSCH 2 OH to formate, and the final product CO 2 , actually released large amounts of energy in the form of NADH which eventually entered into the oxidative phosphorylation process [24,53]. Our results showed that the genes involved in the oxidative phosphorylation process were actually down-regulated significantly ( Figure 3, Table 2). As a result, ATP production should decline in the recombinant P. pastoris during the methanol induction stage in comparison to that cultivated in glycerol medium (PI1 vs. BI, PI2 vs. BI). This phenomenon could be attributed to the down-regulation of genes participating in TCA cycle and glycolysis ( Figure 2). As a result, the total amount of NADH entering into the electron transfer chain decreased accordingly due to the slow-down of the central carbon metabolism. Our conclusion is that the amounts of mitochondrial respiratory chain and the ATP production rate declined in the recombinant P. pastoris after carbon source switching from glycerol to methanol. Moreover, this situation did not change with the methanol induction process ( Figure 3, Table 2, PI2 vs. BI). This could also explain the slowdown of cell growth during the methanol induction phase due to the decrease of ATP production.
The expression of genes related to protein production and degradation were also analyzed ( Figure 4, Table 4). Only the expression of two transcription factors varied significantly in P. pastoris during the methanol induction phase. Moreover, two genes encoding Sec20 and Stx 1-4, which were responsible for secreting foreign proteins involved in SNARE interaction in the vesicular transport pathway, were down-regulated in P. pastoris. Besides we also investigated the expression level of genes involved in protein processing in the endoplasmic reticulum and ER-associated degradation pathway (ERAD). However, most of them were down-regulated. Yeast and other eukaryotic cells would activate the Unfolded Protein Response (UPR) pathway if the unfolded or misfolded proteins accumulated in ER lumen by inducing the expression of genes involved in protein folding and ERAD. The misfolded secretory proteins will be retro-translocated to the cytoplasmic side of the ER, poly-ubiquitinated, and delivered to the proteasome for degradation [54]. Our results showed that the ERAD pathway was not up-regulated and we inferred that the possible reasons were the low molecular weight of HV3 and its relatively simple conformation. As a result, the rate for misfolding would be very low. Besides, gene encoding HSP40, a protein chaperon, was also up-regulated in P. pastoris during the methanol induction stage (Table 4), and it possibly played an important role in HV3 folding. Most of the genes involved in amino acids biosynthesis were also down-regulated ( Figure S5, Table S4). According to the above results, we proposed that the general cellular protein abundance should be reduced even though HV3 was successfully produced by P. pastoris GS115/pPIC9K-hv3 during the methanol induction phase. Hence, little stress of HV3 synthesis and protein secretion occurred in this recombinant yeast but there is still a large space for improving the yield of this recombinant peptide by this strain.
Previous studies also illustrated the profiles of gene expression in P. pastoris cultivated in medium with different carbon sources (glucose, glycerol, or methanol) [55,56]. Two to three folds of genes were highly and selectively expressed in P. pastoris cultivated in methanol medium in comparison to those in glucose-or glycerol-based cultivation. However, genes involved in central metabolism, amino acid metabolism, lipid metabolism, and transport process were all similarly highly-expressed in P. pastoris no matter under which type of cultivation model (glucose-, glycerol, or methanol cultivation). But the genes related to peroxisomal, protein folding, and stress response were more enriched during cultivation in methanol [55]. Liang et al. (2012) also investigated the variation of gene expression in P. pastoris cultivated in methanol medium in comparision to that in glycerol medium. The genes involved in methanol metabolic pathway and proteasome were significantly up-regulated with the methanol induction [56]. However, the above studies all collected the batch culture of P. pastoris in methanol or glycerol to conduct the comparative transcriptomic analysis. We assert that the growth characteristics of P. pastoris in those two mediums are totally different. Thus, a comparison of gene expression in those cultures collected at the same time was not very scientific and rational. Methanol has been used as a carbon source to induce the expression of heterologous protein in P. pastoris for decades [57]. Nowadays, the two carbon source cultivation strategy has been widely applied for the production of heterologous proteins in P. pastoris, and the shift from glycerol to methanol was an important process [52,56]. Hence, it would be more valuable for guiding the host cell engineering for improved heterologous protein production by investigating and analyzing the influence of a carbon source switch on the gene expression in P. pastoris. Actually, our experiment results were partially consistent with the previous studies, such as an enhanced methanol utilization pathway during the methanol induction phase. Moreover, we almost analyzed the expression of every gene involved in the central carbon metabolism, amnio acid metabolic, mitochondrial oxidative phosphorylation process, and so on. Our research outcome presents a promising way to improve the heterologous protein production by enhancing the central carbon metabolism which will eventually improve the utilization of methanol. This would be beneficial for both the cell growth of P. pastoris and heterologous protein production due to the accelerated methanol consumption rate.
In conclusion, we successfully constructed a recombinant P. pastoris GS115/pPIC9K-hv3 for producing HV3, and the antithrombin activity in the fermentation broth reached 5000 ATU/mL. We observed that the biosynthesis of HV3 also became slowed down when the methanol consumption and cell growth rate declined. Our preliminary transcriptomic results showed that glycerate-P converted from methanol enters into the central carbon metabolism, which together with energy metabolism and amino acid biosynthesis, were almost all down-regulated completely in P. pastoris GS115/pPIC9K-hv3 during the methanol induction phase. Thus, we propose that the up-regulation of glycolysis, pentose phosphate pathway, TCA cycle, energy metabolism, amino acid, and lipid metabolism should be beneficial for both improving methanol consumption and HV3 production in P. pastoris GS115/pPIC9K-hv3. Besides, genes involved in EARD, the ubiquitin mediated proteolysis pathway and proteasome were not up-regulated significantly, indicating that the expression of HV3 did not give rise to foreign protein stress and a huge potential still exists for improving the production of HV3 in P. pastoris GS115/pPIC9K-hv3.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/8/606/s1, Figure S1: Reproducibility and reliability analyses of each sample (BI, PI1, PI2). The log 10 (FPKM + 1) were applied to calculate the Pearson correlation coefficient R 2 . The sequencing data with a coefficient R 2 > 0.92 was regarded as high quality, Figure S2: Overview of DEGs for each comparison groups. (A), comparison group: PI1 versus BI; (B), comparison group: PI2 versus BI; (C), comparison group: PI2 versus PI1. The horizontal axis displays the fold change of expression levels of DEGs in different groups, and the vertical axis shows the statistical significance of this variation. Red dots represent genes which are up-regulated, and green dots means the down-regulated parts, Figure S3: Cluster analysis of DEGs in different comparison groups. Sample names are marked at the bottom. Clustering with log 10 (FPKM + 1); red, high-expression level genes; blue, low-expression level genes, Figure S4: Functional annotations of DEGs between different sample groups. (A) Comparable group: PI1 and BI; (B) Comparable group: PI2 and BI; (C) Comparable group: PI2 and PI1. The green bars represent biological process; orange bars represent cellular component; purple bars represent molecular function, Figure  S5: Description of DEGs involved in the amino acid biosynthesis process in PI1 and PI2 in comparison to BI. Up-regulated or down-regulated genes in comparison to BI (p adj < 0.01) are marked with up or down arrow (black arrow stands for sample PI1, red arrow indicates PI2), respectively. A dash line indicates no significant change in comparison to BI, Table S1: Overall quality of the high-throughput sequencing data for each sample, Table S2: Proportion of clean reads mapping to the reference genome, Table S3: Proportion of detected genes at different expression levels, Table S4: List of DEGs involved in amino acid biosynthesis in PI1 and PI2 in comparison to BI, Table S5: Information of DEGs related to the central carbon metabolism, amino acid biosynthesis and function of peroxisome in PI2 in comparison to PI1.