De Novo RNA Sequencing and Transcriptome Analysis of Monascus purpureus and Analysis of Key Genes Involved in Monacolin K Biosynthesis

Monascus purpureus is an important medicinal and edible microbial resource. To facilitate biological, biochemical, and molecular research on medicinal components of M. purpureus, we investigated the M. purpureus transcriptome by RNA sequencing (RNA-seq). An RNA-seq library was created using RNA extracted from a mixed sample of M. purpureus expressing different levels of monacolin K output. In total 29,713 unigenes were assembled from more than 60 million high-quality short reads. A BLAST search revealed hits for 21,331 unigenes in at least one of the protein or nucleotide databases used in this study. The 22,365 unigenes were categorized into 48 functional groups based on Gene Ontology classification. Owing to the economic and medicinal importance of M. purpureus, most studies on this organism have focused on the pharmacological activity of chemical components and the molecular function of genes involved in their biogenesis. In this study, we performed quantitative real-time PCR to detect the expression of genes related to monacolin K (mokA-mokI) at different phases (2, 5, 8, and 12 days) of M. purpureus M1 and M1-36. Our study found that mokF modulates monacolin K biogenesis in M. purpureus. Nine genes were suggested to be associated with the monacolin K biosynthesis. Studies on these genes could provide useful information on secondary metabolic processes in M. purpureus. These results indicate a detailed resource through genetic engineering of monacolin K biosynthesis in M. purpureus and related species.

Red yeast rice (RYR) is produced by M. purpureus fermentation and is also known as red koji [17,18], Hung-Chu, Hong Qu, Ang-kak, Ankak rice, Red Mold rice, or Beni-Koji. RYR has long been recognized for its effects in lowering cholesterol [19,20] and improving resistance of bacteria [21]. In East Asian countries, RYR is traditionally used as a food flavoring agent, coloring agent, preservative, and as a health product to maintain lipid profiles and blood pressure, for instance by normalizing c-reactive protein levels in hypertension [22]. M. purpureus is a commercially important organism owing to its production of the medicinal ingredient monacolin K [23]; the biosynthesis pathway for this compound along with those of citrinin [24] and pigments (fungal polyketone compounds) has been linked to polyketide synthase.
Currently, only about 889 nucleotide sequences and 37 expressed sequence tags (ESTs) can be found for M. purpureus in GenBank (National Center for Biotechnology Information). Moreover, there is little information on functional genes in M. purpureus. The transformation efficiency of M. purpureus has recently improved, which is expected to aid in the generation of genetically modified M. purpureus cell lines. However, the lack of genomic sequence information especially information regarding those of functional genes remains a major obstacle.
High-throughput sequencing has revolutionized genomics research with its low cost and ultra-high data output, including the parallel production of millions of short cDNA reads. RNA sequencing (RNA-seq) is a high-throughput approach that can be used to determine transcript abundance [25] and identify novel transcriptionally active regions [26] even in single cells. This technique is especially important for analyzing the transcriptomes of non-model species [27,28] because it does not require prior knowledge of transcript sequences.
In this study, we investigated the transcriptome of M. purpureus by RNA-seq [29][30][31], which was then annotated using publicly available databases and tools to identify genes involved in the biosynthesis of fungal polyketone compounds that may be related to the tricarboxylic acid cycle [32]. To this end, mRNA was isolated from liquid medium during vegetative growth of M.purpureus, and was used to generate an RNA-seq library. Unigenes were then assembled de novo from sequenced short reads and annotated using the Basic Local Alignment Search Tool for proteins (BLASTX). Gene Ontology (GO) classification and pathway analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was also performed [33][34][35][36][37][38]. The results in our study provide a foundation for the functional characterization of M. purpureus genes.
We performed quantitative real-time PCR (qRT-PCR) to detect the expression of genes related to monacolin K (mokA-mokI) at different phases (2,5,8, and 12 days) of M. purpureus M1 and M1-36 to find the key genes related to monacolin K biogenesis in M. purpureus. The qRT-PCR technology provides a powerful and efficient method for gene function analysis in M. Purpureus.

M. purpureus culture preparation
M. purpureus was obtained from The Chinese General Microbiological Culture Collection Center (Strain Number, CGMCC 3.0568), China. M. purpureus M1 was a wild-type strain. M. purpureus M1-36 was a high producer of monacolin K strain, based on the M1 by mutagenesis. M1 and M1-36 were reactivated by streaking on potato dextrose agar plates from a frozen glycerol stock and were grown for 2-3 days.
Liquid cultures were grown for 2-3 days in a 250 ml flask at 30˚C with shaking at 220 rpm. For liquor fermentation, cells were cultured in a 250 ml flask at 30˚C with shaking at 150 rpm for 2 days for followed by 25˚C with shaking at 150 rpm for 10 days.

M. purpureus collection and RNA isolation
High monacolin K-producing M. purpureus (M1-36) and wild-type M. purpureus (M1) were collected at two growth phases namely rowing and maturation during which secondary metabolites are produced. Cells were immediately frozen in liquid nitrogen and stored at -80˚C until use. Total RNA was isolated using the RNAprep Pure Plant kit (Polysaccharide & Polyphenolic-rich) and DNase I was used to remove contaminating DNA. RNA quality was analyzed using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Oligo (d) T beads (Qiagen, Hilden, Germany) were used to isolate poly(A) mRNA from total RNA.

RNA-seq library preparation and sequencing
To construct the RNA-seq library, 2.5 μg of RNA from four biological replicates was pooled. mRNA enrichment, fragment interruption, addition of adapters, size selection, PCR amplification, and RNA-Seq were performed by LifeTech Co. (Beijing, China). Poly(A) mRNA was isolated using oligo d(T) beads and broken into short fragments. A paired-end RNA-seq library was prepared and sequenced on the Illumina HiSe 2000 platform according to the manufacturer's instructions. In total, 5 G of raw data was obtained.

Analysis of illumina sequencing results
"Dirty" raw reads (i.e., those with adaptors of low quality, with > 5% unknown nucleotides) were discarded. De novo assembly of short reads (from the paired-end RNA-seq library) into a transcriptome was performed with previously described software using Trinity as the default parameter [39]. We first removed low-quality reads from the raw data set and used the Trinity assembler to obtain transcript references; we then aligned the clean reads to the de novo assembly reference. The EdgeR module was used to evaluate differences in expression. Transcripts were aligned with BLASTX, and annotations were performed using the Swiss-Prot database. Only transcripts with E < 1e -5 were selected.

Annotation and classification of unigenes
All unigenes were annotated by a BLASTX search against NCBI nonredundant protein(Nr), Swiss-Prot, Translated European Molecular Biology Laboratory Nucleotide Sequence, COG, NCBI nucleotide database (Nt) and KEGG databases by using a cut-off E-value of 1025. Based on the results of the database annotation, we selected the protein database with the highest sequence similarity to the function annotations. Unigenes were assigned to different pathways by searching against the KEGG databases. Based on the results of the Nr database annotation, the Blast2 GO program was used to obtain the GO annotation of assembled unigenes. WEGO software was then used to plot the GO functional classification for unigenes with a GO term hit to visualize the distribution of gene functions at the macro level.

Identification of transcription factors and target genes in M. purpureus
To identify the transcription factors, all unigenes were searched using Transdecoder analysis tool (http://www.transcriptionfactor.org). The high-throughput sequencing results were predicted by Open Reading Frame (ORF) to obtain the protein, using the Transdecoder software [40]. The HMM motif of the transcription factor was downloaded and predicted by hmmscan in HMMER software to obtain the transcription factor information. Only transcripts with E < 1e -5 were selected.
To extract the sequence from the upstream of transcription initiation site (100-5000 bp) as the promoter region for further search. The site weight matrix motif of the transcription factor were downloaded from http://cisbp.ccbr.utoronto.ca, and the upstream sequences were predicted by the mast program in the MEME software to obtain the target genes [41,42].
qRT-PCR analysis of key genes in monacolin K biosynthesis M. purpureus M1-36 and M1 gene expression profiles were determined by qRT-PCR on a CFX96 qRT-PCR detection system (Bio-Rad, Hercules, CA, USA) with associated software.
cDNA reverse-transcribed from total RNA was used as the template along with SuperReal PreMix Plus with SYBR Green (Tiangen Biotech, Beijing, China) and unigene-specific primers (Table 1) designed with Beacon Designer 8 software (Premier Biosoft). The amplification program was as follows: 95˚C for 15 min, followed by 40 cycles of 95˚C for 10 s, 60˚C for 30 s, and 72˚C for 30 s. Relative abundance of the transcripts was calculated by the comparative cycle threshold method with GAPDH used as an endogenous control. Reactions were prepared in triplicate for each sample.

Results and Discussion
De novo transcriptome assembly and generation of unigenes Samples were collected from rowing and maturation phases. Five cDNA libraries were constructed (M0, M1-1, M1-2, M3, and MW) from approximately 10.6, 12.8, 20.6, 12.7, and 17.9 million raw reads of 200 bp, respectively. In total, 74,557,640 original reads were generated using the HiSeq2000 sequencing platform (Illumina, San Diego, CA, USA). The sequencing error rate (Q20 = 1%) was > 94.00% and the average GC content was 50.33%. Adaptor and duplicate sequences and ambiguous and low-quality reads were removed, and high-quality reads for each cultivar were separately assembled using the NGS QC Toolkit v.2.3 [43] and a series of stringent criteria. After this quality control step, 60,348,338 high-quality reads with a length of 101 bp (average length: 86.62bp) were retained (Fig 1): these were used to assemble contigs and obtain transcripts in default parameter settings by the Trinity program [44]. In total, 41,002 contigs were generated with a k-mer of 25, which was pre-defined to avoid misassembly caused by a too-short k-mer while retaining a reasonable number of reads [45,46]. In the step of assembling, 41,002 transcripts had a median size (N50) of 900 bp and 29,713 unigenes had an N50 of 789 bp from all of the contigs; 39% of the unigenes were not shorter than 500 bp and approximately 200 (0.80%) were longer than 2 kb. In our study, the total length of the assembled transcriptome was 17.8 Mbp (Table 2).

Annotation and classification of M. purpureus unigenes
To identify the putative functions of the assembled unigenes, we determined all of the unigenes by a BLAST search of various databases [47,48]. The unigenes were first searched against the National Center for Biotechnology Information (NCBI; Bethesda, MD, USA) non-redundant protein database (Nr) using BLASTX with a cut-off of P<0.05 and | log2 | !1. The possibility of misassembly was negligible because of the relatively long k-mer (k = 25). In the end, a high percentage of unigenes, i.e., 21,331 of 29,713(71.8%), had a match in this database. A portion of the 8,382 (28.2%) unigenes that did not have a match were presumed to be unique to M. purpureus. Of the 21,331 unigenes with an orthologous match, 15,394 (72.2%) had a match with an annotated function, whereas the remainder had a match that was classified only as a predicted protein.
Additional BLAST searches revealed 5,351, 8,806, 21,331, 18,972, and 12,550 unigenes with matches in the Clusters of Orthologous Groups (COG), Expressed Sequence Tags (EST), Nr, Protein Family (Pfam), and Uniprot databases, respectively (Table 3). In total, 5,351 unigenes had a least match term with the COG databases, whereas 21,331 had a hit against the Nr database (Fig 2).

Differential gene expression in the M. purpureus libraries
We measured the transcript abundance in each sample and identified genes that were differentially expressed between the two M. purpureus cultivars. Clean reads from each sample were mapped to the constructed reference genes, and mapped reads were counted to obtain reads  (Fig 6A-6E).
Heap map of top 50 differentially expressed genes in the five samples prepared using the MeV (MultiExperiment Viewer) program (Fig 6F). The left panel of Fig 6F shows the calculated log2 for clustering of the top 50 differentially transcribed genes at different stages of M1 and M1-36. It shows significant transcriptional differences. Compared with M1-2, the function of up-regulation genes in M3 were mainly concentrated in aspects: metabolic pathways, microbial metabolism in diverse environments, aminoacyl-tRNA biosynthesis, purine metabolism, alanine, aspartate and glutamate metabolism, cysteine and methionine metabolism. Compared with M1-2, the function of down-regulation genes in M3 were mainly concentrated in aspects: metabolic pathways, biosynthesis of secondary metabolites, starch and sucrose metabolism, pyruvate metabolism, fatty acid biosynthesis, methane metabolism.

Identification of transcription factors and target genes in M. purpureus
Transcription factors (TFs) can affect the expression of many genes, also can guide a lot of different biological processes, including cell differentiation, secondary metabolism and so on. Therefore, the identification of transcription factors and target genes will help us obtain a better understanding of genes regulatory networks on monacolin K biosynthesis in M. purpureus. 15,037 protein sequences were obtained after ORF prediction and the default filter screening. 125 transcription factors were finally selected with E < 1e -5 filtering. Finally, 52 transcription factors and 27 target genes were obtained, which were screened according to the p <1e -8 (Fig  7). In Fig 7, the genes in M3 sample were used as the experimental group, and the genes in M1-2 sample as the control group. It could be seen that 4 transcription factors were down-regulated and 1 transcription factors was up-regulated in M3 strain (Table 4). There were also 4 down-regulated target genes and 1 up-regulated target gene in M3 strain (Table 4).

Monacolin K contents in five M. purpureus cultivars determined by highperformance liquid chromatography (HPLC)
We determined from the growth curve of M. purpureus that the adjustment, logarithmic, and stabilization phases were 2, 5, and 12 days, respectively. Monacolin K biosynthesis [50,51] was M1-1 and M1-2 were produced in M1 adjustment and stabilization phases, respectively, whereas M0 and M3 were from the adjustment and stabilization phases of M1-36, and MW was the maximum monacolin K biosynthesis quality phase of M1. No monacolin K was detected in M0 and M1-1, and the content was higher in M3 (295.6101 mg/l) than in M1-2 (117.7263 mg/l) and MW (83.9650 mg/l) among the five M. purpureus cultivars. Monacolin K contents were higher in M1-2 and M3 than in M1-1, M0, and MW (Fig 8).
The expression of nine genes encoding different enzymes involved in monacolin K biosynthesis were verified by qRT-PCR (Fig 9). The expression patterns of the nine genes differed between M1 and M1-36. Moreover, mokC, mokD, mokE, mokF, mokH and mokI levels showed strong positive correlations with monacolin K concentration as compared to the observed expression profiles and previously measured monacolin K contents. These results suggest that differences in the expression patterns of M.purpureus genes account for monacolin K diversity.
To better understand the biological functions of genes on monacolin K biosynthesis in M. purpureus, monacolin K biosynthetic pathway were drawn (Fig 10). This result shows that M. purpureus M1-36 contained a higher monacolin K proportion than M1, suggesting that monacolin K synthesis was more active in M1-36 than in M1. The transcriptome analysis indicated  that under high monacolin K production conditions, most genes in this polyketide synthase (PKS) cluster, especially the PKS gene (mokD), are expressed at high levels, which further supports the key role of this gene cluster in monacolin K synthesis. MokD and its homolog gene ctnB are belonged to the PKS gene clusters, and play as oxidoreductase in the synthesis pathway of their respective monacolin K and citrinin. Li et al found that the ctnB-deficient mutant in M. aurantiacus barely produced citrinin, indicating that ctnB gene is directly involved in citrinin biosynthesis [53]. Therefore, it is speculated that the mokD gene plays an important role in monacolin K synthesis. Currently, the function of mokD gene is being validated at our lab.

Conclusions
In this study, we obtained the transcriptome information of M. purpureus, including functional annotation and classification, which could provide insight into the role of intracellular metabolic pathways [54]. In this study, we examined the functions of M. purpureus unigenes by BLAST searches against various databases [55]. Illumina HiSeq 2000 sequencing generated 60 million high-quality short reads, which were assembled into 29,713 assembled unigenes in M. purpureus. In summary, 122,599 unigenes showed a match with the GO annotation. Meanwhile, biosynthetic pathways involving some of these unigenes have been identified. These results contribute to existing sequence resources of M. purpureus and will be useful for studies on the biosynthesis of secondary metabolites such as aminoacyl-tRNA, nitrogen metabolism, and terpenoids in M. purpureus.
Owing to the economic and medicinal importance of M. purpureus [56], most studies on this organism have focused on the optimization of fermentation conditions, the physiological activity of secondary metabolites, and the PKS clusters genes of monacolin K biosynthesis. However, few M. purpureus genes have been characterized. One study found that mokF modulates monacolin K biogenesis in M. purpureus. In this study, nine genes (mokA-mokI) were found to be associated with the biosynthesis of monacolin K. The research of these genes could give us a deeper understanding on monacolin K biosynthesis and other secondary metabolic processes in M. purpureus.