Physiological and Transcriptomic Analysis Provide Insight into Low Temperature Enhancing Hypericin Biosynthesis in Hypericum perforatum

Hypericin (Hyp), well-known as an antidepressant, is mainly extracted from Hypericum perforatum. Although Hyp accumulation and biomass are greater at lower compared to higher temperature, the regulation mechanism has not been reported. Here, the physiological characteristics and transcriptome of H. perforatum grown at 15 and 22 °C were determined and analyzed by HPLC and de novo sequencing. The results showed that the stomatal density and opening percentages were 1.1- and 1.4-fold more, and the Hyp content was 4.5-fold greater at 15 °C compared to 22 °C. A total of 1584 differentially expressed genes (DEGs) were observed at 15 versus 22 °C, with 749 characterized genes, 421 upregulated (UR) and 328 downregulated (DR). Based on biological functions, 150 genes were associated with Hyp biosynthesis, plant growth and the stress response, including photosynthesis, carbohydrate metabolism, fatty acids metabolism, cytochrome P450 (CYPs), morpho-physiological traits, heat shock proteins (HSPs), cold-responsive proteins (CRPs) and transcription factors (TFs). The differential expression levels of the master genes were confirmed by qRT-PCR and almost consistent with their Reads Per kb per Million (RPKM) values. This physiological and transcriptomic analyses provided insight into the regulation mechanisms of low temperature enhancing Hyp biosynthesis in H. perforatum.

Hyp biosynthesis is started from glucose to acetyl-and malonyl-coenzyme A (CoA) in the cells of green tissues and from acetyl-and malonyl-CoA to Hyp in the cells of dark glands ( Figure 1). Briefly, in green tissue, photosynthesis firstly produces glucose, and pyruvic acid is produced through glycolysis. Acetyl-CoA can be formed through pyruvate dehydrogenase. Malonyl-CoA can be derived from acetyl-CoA by acetyl-CoA carboxylase. In addition, acetyl-CoA can also be recycled as the result of the fatty acid metabolism. In the dark gland, one molecule of acetyl-CoA and seven molecules of malonyl-CoA form an octa-β-ketoacyl chain by polyketide synthase (PKS) or octaketide synthase (OKS); emodin anthrone was formed after several steps, including aldolic condensation, thioesterase (TER), decarboxylic and dehydration reactions. Emodin dianthrone can be generated through two pathways, including the oxidation of emodin bianthrone, which is formed by the free radical coupling of emodin anthrone, and the synthesis of emodin anthrone and emodin by the phenoloxidative coupling protein (POCP). Then, Hyp is subsequently yielded via the catalysis of the berberine bridge enzyme (BBE) [9][10][11]22,[24][25][26]. Experimental evidence shows that a low temperature improves Hyp production. Specifically, Brunáková et al. [27] reported that a −4 • C treatment improved the Hyp content on a dry weight (DW) basis. Tavakoli et al. [20] found that the Hyp content on a DW basis in a root suspension culture reached the highest amount at 4 • C than at 8, 16 and 25 • C. Velada et al. [28] determined that the genes (phenolic oxidative coupling protein (Hyp-1), PKS1 and PKS2) involved in Hyp synthesis were differentially expressed at different temperatures. Our previous studies found that 15 • C could enhance aerial parts of a biomass and Hyp contents of dry weight (DW) on a per plant basis, as well as mRNA expressions of the hyp-1 and PKS genes, compared to 22 • C [22]. For the molecular function of hyp-1 gene, Bais et al. [29] found that a conversion of emodin to Hyp was solely catalyzed by the hyp-1 gene, while Košuth et al. [30] reported that the hyp-1 gene is not a limiting factor for Hyp biosynthesis in the genus Hypericum or that additional factors are necessary for completion.
To date, the mechanism of low temperature regulating plant growth and Hyp biosynthesis in H. perforatum is still limited. With the aim of unraveling the genes involved in growth enhancement and Hyp accumulation upon low temperature induction, the seedlings grown at 15 and 22 • C were examined by RNA-seq, the mRNA expression levels were confirmed by qRT-PCR and the Hyp content was quantified by HPLC.

Effect of Low Temperature on Biomass, Dark Glands, Leaf Stomata and Hyp Accumulation
A significant increase of the DW of whole plants was observed at 15 • C than 22 • C (Figure 2A,C). Although 1.1-fold more of the number of foliar dark glands was observed at 15 • C than 22 • C, there was no significant difference between them ( Figure 2B,D). Four point six and 4.5-fold increases of the Hyp content, on a DW and per plant basis, were observed at 15 • C compared to 22 • C, respectively ( Figure 2E,F). Additionally, there were significant differences in the leaf stomata, with the stomatal density and opening percentage 1.1-and 1.4-fold more ( Figure 3A,C and Table 1 and a more rounded shape with the ratio of length/width (2.7) at 15 • C compared to 22 • C ( Figure 3B,D and Table 1).

Global Gene Analysis
To dissect the mechanisms of low temperatures enhancing plant growth and Hyp accumulation, RNA-seq for plants grown at 15 and 22 • C was performed. Although robust data were generated ( Figure S1), 39.4 and 45.5 million high-quality reads were collected after data filtering, and 2.5 and 2.9 million multiple reads were mapped from the 15 and 22 • C plants, respectively, from the 44,776 compiled genes and annotated against the databases (Tables S2 and S3 and Figures S2-S5). The correlation between 15 and 22 • C was assessed by a principal component analysis (PCA) with R packages (Figures S6 and S7). A total of 1584 differentially expressed genes (DEGs) were identified at 15 versus 22 • C ( Figure S8). Of these 1584 DEGs, 768 genes were identified to match with the databases ( Figure 4B). Among the 768 genes, 749 genes with known functions were partitioned into 421 upregulated (UR) and 328 downregulated (DR) ( Figure 4C,D).

Discussion
While low temperatures improving the plant growth and Hyp accumulation by stimulating related gene expressions has been observed in H. perforatum [20,22,27,28], the mechanism responsible for temperature-dependent growth and Hyp biosynthesis has not been dissected.
Here, it is shown that there is a greater plant biomass and Hyp accumulation, more leaf dark glands, stomatal density and rounded stomatal shape at 15 compared to 22 • C. By treating seedlings at 15 and 22 • C, ca. 1600 genes were differentially expressed, with over 56% of the 749 characterized genes being UR at 15 • C compared to 22 • C ( Figure 4). By grouping genes based on biological functions, 150 genes were associated with Hyp biosynthesis, plant growth and stress responses, including photosynthesis, carbohydrate metabolism, fatty acid metabolism, CYPs, morphophysiological traits, HSPs, CRPs and TFs ( Figure 5).
Experimental evidence shows that H. perforatum plants can adapt to low temperatures by increasing the biomass and morphogenesis (e.g., leaf area, number of stem and its internode and root), and low temperatures can improve Hyp biosynthesis by increasing the number of dark glands and related gene expressions [10,12,20,22,27,28]. The current study exhibited greater biomass by increasing the stomatal density and opening percentages, as well as a more rounded stomatal shape with the ratio of length/width and greater Hyp content via the overexpression of 61 genes involved in Hyp biosynthesis (Figures 2 and 3 and Table 1).
As is shown in Figure 1, the Hyp biosynthetic pathway depends on photosynthesis and is involved in carbohydrate metabolism and fatty acid metabolism in the upstream steps, as well as PKS, TER, POCR and CYPs. For the photosynthesis, the CABs are the most abundant protein complex on the thylakoid membrane and are light receptors that capture and deliver excitation energy to photosystems [31]. The UR of CAB genes that are involved in photosynthesis provides energy and fixation and a reduction of CO 2 into carbohydrate [32].
For the fatty acid metabolism, for example, the phospholipase D zeta 1 (PLDZETA1) is involved in the lipid catabolic process by the hydrolysis of glycerol-phospholipids to produce phosphatidic acids [39]. The acetyl-CoA-benzyl alcohol acetyltransferase (BEAT) is involved in the biosynthesis of benzyl acetate [40]. The 2,4,6-trihydroxybenzophenone synthase (BPS) is involved in the malonyl-CoA metabolic process [41]. 3-ketoacyl-CoA thiolase 2 (PED1) is involved in long chain fatty acid beta-oxidation [42]. The fatty acid metabolic process will provide abundant acetyl-CoA, which can combine with its product malonyl-CoA to produce an octa-β-ketoacyl chain.
Extensive studies have demonstrated that CYPs can catalyze several reactions, such as hydroxylation, dealkylation and oxidation reactions. Studies have demonstrated that the CYPs are a large and diverse family involved in secondary metabolism [43][44][45][46][47]. During the downstream of the Hyp biosynthetic pathway in dark glands, most of the steps are catalyzed by the oxidation, reduction, dehydration and decarboxylic reactions, as well as free radical coupling and electron-paired donors (Figure 1). For example, flavonoid 3 -monooxygenase (CYP75B2) catalyzes the 3 -hydroxylation of flavonoid B-ring to a 3 , 4hydroxylated state [48], isoflavone 2 -hydroxylase (CYP81E1) catalyzes the hydroxylation of isoflavones to 2 -hydroxyisoflavones, daidzein to 2 -hydroxydaidzein and formononetin to 2 -hydroxyformononetin [49]. The UR of genes that encode CYPs may play critical roles in Hyp biosynthesis in the downstream pathway.
However, the genes involved in downstream steps (e.g., PKS, OKS, TER, POCP and BBE) were not examined to be differentially expressed, except for the CYPs class ( Figure 1). While several studies were devoted to dissecting the Hyp biosynthetic pathway [9][10][11]22,[24][25][26], limited genes involved in the downstream steps were submitted to the public databases (e.g., NR, Swiss-Prot, KEGG, KOG and GO), which led to differentially expressed contigs that could not be mapped and were listed as unidentified or uncharacterized genes in this study.
Adaptation to environmental stresses results from integrated events occurring at all levels of organization, not only from the anatomical, cellular, biochemical and molecular levels but, also, from the morphological level [32]. For example, the alkane hydroxylase MAH1 (CYP96A15) is involved in wax biosynthesis during plant growth [50]. The Leucine-rich repeat receptor protein kinase (EMS1) can mediate signals that control the fate of reproductive cells and their contiguous somatic cells [51]. The GLABRA2 expression modulator (GEM) is involved in the spatial control of cell division and differentiation of root epidermal cells [52]. MIZU-KUSSEI 1 (MIZ1) plays an important role in lateral root development [53]. Subtilisin-like serine proteases (SBTI1.1) can promote plant cell differentiation, organogenesis and somatic embryogenesis, as well as cell proliferation [54]. Zeatin O-glucosyltransferase (ZOG1) catalyzes the formation of O-xylosylzeatin in plant development [55]. The differential expression of the 15 genes involved in the morphophysiological traits may be an adaptation to low temperatures.
There is considerable evidence that HSPs are essential for plant survival under abiotic stresses [56,57]. The expression of HSP genes is observed to be induced by high temperatures and inhibited by low temperatures [58,59]. For example, CRPs are induced to enable plants to survive from chilling and freezing stress [60]. Three CRPs that were identified in this study-for example, BON1-associated protein 1 (BAP1), cold-regulated 413 inner membrane protein 1 (COR413IM1) and cold-responsive protein kinase 1 (CRPK1)-were observed to be differentially expressed in response to cold stress [61][62][63].
TFs can regulate plant growth and metabolism in response to abiotic stresses [64]. For example, ethylene-responsive transcription factors (ERFs) are involved in regulating gene expression and participating in signal transduction pathways [65,66]. MYB transcription factors play regulatory roles in developmental processes and defense responses, the MYB4 is involved in scavenging reactive oxygen species during cold stress [67] and MYB14 induces the expression of CBF genes to increase their freezing tolerance [68]. The NAC domain-containing protein family may act downstream of N-rich protein (NRP)-A or NRP-B in the integration of endoplasmic reticulum stress and osmotic stress cell death signals [69]. WRKY domain transcription factors play an important role in plant responses to both biotic and abiotic stress by regulating ABA signaling [70].
Based on the above functional analysis, the schematic representation of the proposed pathways of genes for regulating Hyp biosynthesis, morphophysiological traits and the stress response are shown in Figure 7. Briefly, low temperatures induce the differential expression of genes involved in photosynthesis, carbohydrate metabolism and fatty acid metabolism, which improves the accumulation of primary metabolites, as well as the metabolic processes of carbohydrate and fatty acids. Hyp biosynthesis is improved by the overexpression of genes encoding CYPs, as well as other unidentified enzymes ( Figure 7A). Genes involved in the morphophysiological traits (e.g., root development, leaf and stem development and cell differentiation) are also upregulated, which improves the plant growth to some extent under low temperatures ( Figure 7B). Meanwhile, low temperatures inhibit the expression of HSPs and induce the overexpression of CRPs, which subsequently regulates the differential expression of TFs to adapt H. perforatum plants to low temperature conditions ( Figure 7C).

Plant Materials
Mature seeds of Hypericum perforatum L. were obtained from plants cultivated in Kangxian County (33 • 16 20" N, 105 • 31 50" E) of Gansu Province, China in July 2016. Seedlings were induced based on a previous reported protocol [22] with slight modifications. Specifically, seeds were rinsed with tap water for 10 min and successively immersed in 70% ethanol for 2 min and 0.1% HgCl 2 for 1 min. After each treatment, seeds were rinsed with sterile water 4 times. Sterilized seeds were inoculated on MS + 20.0 g/L sucrose + 4.0 g/L agar (pH 5.8) and germinated at 22 • C with a 24-h/d photoperiod (500 µmol·m −2 ·s −1 flu). After 25 days, germinated seeds were transferred to the abovementioned basal medium (pH 5.8) supplemented with 0.5-mg/L NAA and 1.0-mg/L 6-BA for rooting and sprouting. After 20 days of growth, seedlings were treated at a constant 15 or 22 • C in a growth chamber, with each treatment containing 40 flasks (3 seedlings per flask). After 20 days of treatment, seedlings were harvested with 20 flasks used for physiological measurements and Hyp quantification and 20 flasks used for RNA-seq and gene expression.

Observation of Dark Gland and Leaf Stomata
The number of dark glands was calculated on a per leaf basis, with 25 leaf replicates. Stomata in the middle adaxial part of the leaf were observed by a scanning electron microscope (S-3400N, Hitachi, Tokyo, Japan), and this step was applied following a previous protocol [71].

Hyp HPLC Quantification
Hyp content was quantified according to previous protocols [22,72]. Air-dried seedlings were ground into powder, and a 0.1-g sample was used to extract Hyp. Ethanol extract (10 µL) was analyzed by HPLC (symmetry C18, 250 mm × 4.6 mm, 5 µm and column temperature 60 • C; Waters ACQUITY Arc, Waters Corporation, Milford, MA, USA). Hyp content was evaluated with a peak area comparison with a reference standard (hypericin, 56690; Sigma Chemical Co., St. Louis, MO, USA).

Total RNA Extraction and Illumina Sequencing
Total RNA samples of 15 or 22 • C with three biological replicates were extracted using a RNA kit (R6827, Omega Bio-Tek, Inc., Norcross, GA, USA) according to the manufacturer's protocols. To ensure the quality for sequencing, the integrity of total RNA was evaluated by 1.0% agarose gel electrophoresis ( Figure S9), the purity of total RNA was determined by a NanoDrop spectrophotometer 2000 (Thermo Fisher Scientific, Waltham, MA, USA), the concentration of total RNA was quantified by a Qubit2.0 Fluorometer (Thermo Scientific) and the accurate integrity of the total RNA was detected on an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA). The processes of enrichment, fragmentation, reverse transcription, synthesis of the second-strand cDNA and purification of cDNA fragments was applied following the previous protocols [47]. RNAseq was conducted using an Illumina HiSeqTM 4000 platform (Gene Denovo Biotechnology Co., Ltd., Guangzhou, China).

Sequence Filtration, Assembly, Unigene Expression Analysis and Basic Annotation
Raw reads were filtered according to previous descriptions [47]. Specifically, read pairs were used to filter with the principles: (1) remove reads containing adapters, (2) remove reads containing more than 10% of unknown nucleotides (N) and (3) remove low-quality reads containing more than 50% low-quality (Q-value ≤ 20) bases. De novo assembly of clean reads was performed using Trinity [73]. The expression level of each transcript was calculated and normalized to Reads Per kb per Million (RPKM) [74], and DEGs were identified according to a criteria of |log 2 (fold-change)| ≥ 1 and p ≤ 0.05 by DESeq2 software and the edgeR package [75,76]. Unigenes were annotated against the databases, including NR, Swiss-Prot, KEGG, KOG and GO [77].

qRT-PCR Validation
The primer sequence (Table S1) was designed on the primer-blast of the NCBI website. First-strand cDNA was synthesized using a RT kit (KR116, Tiangen, China). PCR amplification was performed using a SuperReal PreMix (FP205, Tiangen, China). Melting curve was analyzed at 72 • C for 34 s. Actin gene (Accession no. CP002685.1) was used as a reference control. Relative expression level (REL) was calculated using the 2 − Ct method [78].

Statistical Analysis
All the measurements were performed using three biological replicates. A t-test was applied for independent samples, with p < 0.05 considered significant.

Conclusions
From the above observations, lower temperatures improve the plant growth and Hyp accumulation in H. perforatum. The DEGs observed in H. perforatum grown at different temperatures strongly indicate there is a transcription-based regulation of Hyp biosynthesis at lower temperatures. The master genes involved in Hyp biosynthesis, the morphophysiological traits and the stress response were identified, and their causative roles should be investigated by gene editing, like CRISPR/Cas 9 systems.
Supplementary Materials: The following are available online, Figure supplemental legends: Figure S1: Length distribution of unigenes; Figure S2: Basic annotation for all unigenes; Figure S3: Top 10 plant species distribution of the total homologous sequences; Figure S4: Distribution of unigenes in the transcriptome with KOG functional classification; Figure S5: Functional classifications of GO terms of all assembled unigenes; Figure S6: The PCA analysis between 15 and 22 • C samples; Figure S7: The correlation analysis of heat map between 15 and 22 • C samples; Figure S8: Number of differentially expressed genes and volcano plot when plants grown at 15 versus 22 • C; Figure S9: The quality of total RNA at 15 and 22 • C was evaluated by 1.0 % agarose gel electrophoresis. Table supplemental legends: Table S1: Sequences of primer employed in qRT-PCR analysis; Table S2: Summary of sequencing data of H. perforatum transcriptome;   Institutional Review Board Statement: Not applicable. This study did not involve any humans or animals.
Informed Consent Statement: Not applicable. This study did not involve any humans or animals.