*Article* **Physiological, Metabolic and Transcriptional Responses of Basil (***Ocimum basilicum* **Linn. var.** *pilosum* **(Willd.) Benth.) to Heat Stress**

**Lei Qin 1,2,†, Chengyuan Li 1,3,†, Dongbin Li 4, Jiayan Wang 1, Li Yang 1, Aili Qu 1,\* and Qingfei Wu 1,\***


**Abstract:** As a medicinal and edible plant, basil (*Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth.) has rich nutrition and significant economic value. The increase in heat stress caused by global warming adversely affects the growth and yield of plants. However, the response mechanism of basil to heat stress is poorly understood. This work investigated the changes in phenotype, metabolome, and transcriptome in basil under heat stress. The results showed that heat stress triggered severe oxidative damage and photosynthesis inhibition in basil. Metabonomic analysis showed that, compared to the control group, 29 significantly differentially accumulated metabolites (DAMs) were identified after 1 d of heat treatment, and 37 DAMs after the treatment of 3 d. The DAMs were significantly enriched by several pathways such as glycolysis or gluconeogenesis; aminoacyl-tRNA biosynthesis; and alanine, aspartate, and glutamate metabolism. In addition, transcriptomic analysis revealed that 15,066 and 15,445 genes were differentially expressed after 1 d and 3 d of heat treatment, respectively. Among them, 11,183 differentially expressed genes (DEGs) were common response genes under 1 d and 3 d heat treatment, including 5437 down-regulated DEGs and 6746 up-regulated DEGs. All DEGs were significantly enriched in various KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, most dominated by glyoxylate and dicarboxylate metabolism, followed by starch and sucrose metabolism, and by the biosynthesis and metabolism of other secondary metabolites. Overall, all the above results provided some valuable insights into the molecular mechanism of basil in response to heat stress.

**Keywords:** basil; heat stress; metabolome; transcriptome; molecular mechanism

#### **1. Introduction**

*Ocimum basilicum* (basil) is an annual aromatic plant belonging to the Lamiaceae family, widely cultivated in many countries, and has significant economic value [1]. Many important compounds in basil have been isolated and identified, such as volatile oils, and various polyphenols and flavonoids with antioxidant, anti-inflammatory, and other pharmacological activities [2–4]. Therefore, basil plants are often used as seasoned vegetables in cooking and as ingredients in some medical products industry [5]. The herb *Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth., as a variety of *Ocimum basilicum*, is widely distributed throughout China [6,7], contains many kinds of essential oils, and has been used as medicine to treat some diseases in folk prescriptions [7]. To the best of our knowledge, however, very few studies have focused on the heat tolerance of this plant.

As one of the common and serious abiotic interferences, heat stress imposes negative effects on the normal growth of plants, and damages agricultural production, which is

**Citation:** Qin, L.; Li, C.; Li, D.; Wang, J.; Yang, L.; Qu, A.; Wu, Q. Physiological, Metabolic and Transcriptional Responses of Basil (*Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth.) to Heat Stress. *Agronomy* **2022**, *12*, 1434. https:// doi.org/10.3390/agronomy12061434

Academic Editors: Sara Álvarez and José Ramón Acosta-Motos

Received: 26 May 2022 Accepted: 14 June 2022 Published: 15 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

becoming a worldwide problem [8,9]. It is estimated that the global air temperature will rise by about 0.2 ◦C per decade and thus increase another 1.8–4 ◦C by the year 2100 [10–12], so the world has to face more severe challenges of heat stress. By limiting various physiological and biochemical responses, heat stress could adversely impact plant development. In general, heat stress reduces photosynthetic efficiency and injures photochemical reactions and carbon metabolism, thereby diminishing the plant productivity and shortening the life cycle [13–15]. The functions and integrity of cell membranes are also highly sensitive to heat stress, and cellular membrane permeability and cellular electrolytes would be increased under heat stress, leading to membrane instability [16–19]. In addition, heat stress could induce oxidative stress to produce excessive reactive oxygen species (ROS) in plant cells, while ROS formation would give rise to various negative effects, such as lipid peroxidation, nucleic acid damage, enzyme inactivation, and protein oxidation [20]. All these factors would degenerate physiological and biochemical processes of a plant under heat stress, limit plant growth, reduce its yield, and even cause its death.

On the other hand, plants have evolved a series of strategies in adaption to heat stress. There are two ways to mitigate deleterious effects due to heat stress. One is to prevent heat stress by reducing the absorption of solar radiation; the other is to enhance plant tolerance to heat stress [11]. Generally, plants can avoid short-term heat stress through transpirational cooling, altering membrane lipid compositions, or changing leaf orientation. For the plants grown in hot climates, small hairs (tomentose) have always been found to form a thick coat on the surface of the leaf and cuticles, contributing to reduce the absorption of solar radiation. Moreover, plants can deal with heat stress by regulating a lot of physiological processes. For instance, plants under heat stress would accumulate many more osmotic protection substances including proline, polyols, and soluble carbohydrates, which could lower the osmotic potential of cells and attenuate the stress effects [21]. In addition, plants have formed a defense system of antioxidation to scavenge excessive ROS, which contain nonenzymatic antioxidants, such as ascorbic acid (vitamin C), phenolic compounds, flavonoids, carotenoids, and other secondary metabolites, as well as enzymatic antioxidants including superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) [22]. Furthermore, the responses of a plant to heat stress are also mediated by its own hormones, including abscisic acid (ABA), jasmonic acid (JA), ethylene (ET), cytokinins (CKs), brassinosteroids (BRs), and salicylic acid (SA) [23], and these hormones participate in responses to heat stress mainly through triggering stress-adaptive signaling cascades [24,25].

At the molecular level, it has been reported that a large number of genes such as heat shock protein (HSP) genes, osmosensors, mitogen-activated protein kinase (MAPK) cascades, transcription factors, and second messengers would participate in the plant response to heat stress [26,27]. Currently, RNA-sequencing (RNA-Seq) is becoming a fast and effective way to uncover gene expression profiles, so it has been widely used to explore stress response mechanisms. A series of studies have revealed basil's response to high temperatures. For example, when subjected to a short-term heat stress, the contents of phenolic compounds and activities of antioxidant enzymes were remarkably increased in *Ocimum basilicum* L. cv. 'Genovese' [28]. However, when 'Genovese' was grown under heat stress and elevated CO2 concentrations for a long time, the content of phenolic acid in the leaves decreased, the wax decreased, and the photosynthesis increased [9]. Due to basil displaying a significant number of varieties with different genetic backgrounds, their resistance to heat stress may also be different, and the molecular mechanism of basil's response to heat stress also needs further investigation. In this study, the changes in physiological features and metabolites in basil (*Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth.) under heat stress were investigated and analyzed, and RNA-Seq was used to identify differentially expressed unigenes. The results will be helpful for cloning and analyzing heat-resistant genes in basil, theoretically providing some insights for molecular breeding in the future.

#### **2. Materials and Methods**

#### *2.1. Plant Materials and Treatment*

The seeds of *Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth. were provided by Zhejiang Academy of Agricultural Science (Hangzhou, China). Seeds were sown in plastic pots filled with turfy soil, which was rich in incompletely decomposed plant residues, humus, and some minerals, and the pots were placed in an artificial illumination incubator. The cultivation conditions were set as follows: a photoperiod of 12 h light (30 ◦C)/12 h dark (24 ◦C), air relative humidity of 60%, and an illumination of 300 μmol m−<sup>2</sup> s−1. To ensure the healthy and normal growth of plants, half-strength Hoagland nutrient solution was irrigated daily. For heat treatment, the seedlings of eight weeks were exposed to a day/night temperature of 42/36 ◦C, and the leaves were collected as heat-stressed samples at 0, 1, 3, and 5 d after treatment. For each treatment, three biological replicates were performed, and each biological replicate contained 15 individual plants. Then, the collected samples were immediately frozen in liquid nitrogen and stored at −80 ◦C before use. To ensure the consistency of the experiment, the leaves of each sample were taken from the top to bottom of the third and fourth leaves, respectively.

#### *2.2. Physiological Analysis*

After heat treatment with different times (0, 1, 3, and 5 d), the contents of oxidizing substances including malonaldehyde (MDA), H2 O2, and O2 − were determined according to the instructions of reagent kits (Beijing Solarbio Science and Technology, Beijing, China). In brief, for MDA content measurement, 0.1 g of leaf tissue was weighed and 1 mL of extract buffer was added for ice bath homogenization. Centrifuging at 8000× *g* for 10 min at 4 ◦C, the supernatant was taken and placed on ice for testing. Finally, the absorbance of the supernatant at 532 nm (A532) and 600 nm (A600) was measured, and MDA content = 32.258 × (ΔA532–ΔA600) ÷ weight. The H2O2 content was calculated based on the absorbance of the supernatant at 415 nm, and the content = Δ Asample ÷ Δ Astandard ÷ weight. As the purple-red azo compound has a characteristic absorption peak at 530 nm, the O2 − content in the sample can be calculated according to the absorbance at 530 nm (A530) and standard curve. Photosynthetic parameters were measured using a portable photosynthesis system (Licor-6400, LI-COR, Lincoln, NE, USA), and all measurements were performed at a photon flux density of 1000 μmol m−<sup>2</sup> s−1, CO2 concentration of 400 μmol mol−1, and leaf temperature of 30 ◦C. Chlorophyll fluorescence parameters were measured using a Dual-PAM 100 chlorophyll fluorescence analyzer (Heinz Walz, Effeltrich, Germany) by the manufacturer's instructions. Before measurement, plants need to be dark-adapted for at least half an hour. The measurement was repeated six times and averaged.

#### *2.3. Extraction and Detection of Metabolites*

The extraction of metabolites was conducted following a previous method [29]. In detail, the sampled leaves were quickly ground into powder in a mortar cooling with liquid nitrogen, and the accurately weighed 50 mg powder was put into a 2 mL centrifuge tube. Next, 0.5 mL of the mixed solution at acetonitrile: isopropanol: water (3:3:2, *v*/*v*/*v*) and several zirconium beads of 2 mm were added in turn into the tube, which was placed into the high-throughput tissue grinding instrument, shaken at 30 Hz for 20 s, and allowed to stand for 10 s. This was repeated 8 times and then sonicated in an ice-water bath for 5 min. Once again, 0.5 mL of the mixed solution of acetonitrile: isopropanol: water (3:3:2, *v*/*v*/*v*) was added and then sonicated in an ice-water bath for 5 min. After centrifugation at 14,000 rpm for 2 min, the supernatant solution of 500 μL was pipetted into a new 2 mL tube and concentrated to dryness with a vacuum concentrator for 8–10 h, while the rest of the supernatant was placed into a refrigerator at −80 ◦C for backup. The as-concentrated sample was redissolved with 80 μL of methoxypyridine solution of 20 mg/mL concentration, and followed by 30 s of vortex vibration and 60 min of incubation at 60 ◦C. Sequentially, 100 μL of BSTFA-TMCS (99:1) derivatization reagent was added, vortexed for 30 s, incubated at 70 ◦C for 90 min, and centrifuged at 14,000 rpm for 3 min, and the supernatant of 100 μL was put into the detection bottle. The final samples were temporarily stored in a sealed cup for testing, and GC-TOF detection would be completed within 24 h.

Samples were analyzed by gas chromatography with a DB-5MS capillary column (Agilent, Santa Clara, CA, USA). In the analysis, the constant flow of helium was at 1 mL/min. A 1 μL sample was injected through the autosampler in a split ratio of 1:10, helium was used as the carrier gas at constant flow rate of 1 mL/min, the inlet temperature was 280 ◦C, and the transfer line and ion source temperatures were 320 ◦C and 230 ◦C, respectively. The heating program initially started at 50 ◦C for 0.5 min, increased to 320 ◦C at a rate of 15 ◦C/min, and remained at 320 ◦C for 9 min. Mass spectrometry was carried out with a full scope at a scan rate of 10 spec/s, an electron energy of −70 V, and a solvent delay of 3 min.

#### *2.4. Metabolomics Analysis*

The raw data of mass spectrometry were converted into net CDF format by Agilent MSD Chem Station for peak identification, peak filtering, and peak alignment, and the data matrix comprised the mass to charge ratio (*m*/*z*), retention time, and intensity. Metabolites were annotated in conjunction with the AMDIS program using the National Institute of Standards and Technology (NIST) commercial database and Wiley Registry metabolome database. Among them, the alkane retention index was used for further substance characterization according to the retention index provided by The Golm Metabolome Database (GMD) (http://gmd.mpimp-golm.mpg.de/ (12 October 2021). Then, the metabolomic data were analyzed using the MetaboAnalystR package in R language [30].

#### *2.5. cDNA Library Construction and Sequencing*

The leaves of each sample were fully ground into powder under liquid nitrogen cooling. Then, about 0.2 g of powder was weighed for total RNA extraction using the RNAprep Pure Plant Plus Kit (Tiangen, Beijing, China). The RNA purity and concentration were measured using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Approximately 1.5 μg of RNA per sample was used to construct the sequencing library according to the previously described method [31]. Then, the library preparations were sequenced to generate 150 bp paired-end reads using the Illumina HiSeq 4000 platform (Illumina, San Diego, CA, USA). The raw sequencing data were uploaded to the NCBI BioProject database, given a project number of PRJNA833588.

#### *2.6. Transcriptomic Analysis*

The raw data obtained by sequencing included a small number of reads with sequencing adapters or with low sequencing quality. In order to ensure the quality and reliability of data analysis, the original data should be filtered, and the adaptor contamination and low-quality reads were removed using Trimmomatic software (version 0.38). After removing unclear reads (unknown nucleotides > 5%) and low-quality reads (Q-value < 20), a clean sequence was obtained. The clean reads were assembled de novo using Trinity software, and the redundancies in the assembled transcript sequences were removed using CD-HIT-EST software with the global sequence identity threshold cut-off setting at 90%.

Uni genes were functionally annotated with alignments to the KOG, NT, R, Uniprotand Swiss-Prot databases. The DESeq R package (1.10.1) was used to analyze the differentially expressed genes (DEGs) with an adjusted *p*-value of <0.05. GO annotation and KEGG annotation were completed by a Python script and the GhostKOALA website (https: //www.kegg.jp/ghostkoala/ (accessed on 9 November 2021), respectively. GO and KEGG enrichment analysis was performed by TBtools and Mapman [32].

#### *2.7. Quantitative Real-Time PCR (qRT-PCR) Analysis*

First-strand cDNA was synthesized from 1 μg of RNA using the PrimeScript RT reagent Kit with gDNA Eraser (Takara, Dalian, China). qRT-PCR analyses were conducted using SYBR® Premix Ex Taq™ II (Takara, Dalian, China) on an ABI Prism 7700 (Applied Biosystems, Foster City, CA, USA). The PCR procedures went through at 95 ◦C for 10 s, followed by 40 cycles at 95 ◦C for 5 s, and then at 60 ◦C for 30 s. The gene-specific primers were obtained in Primer Premier 5.0, and all primer sequences are listed in Table S1. The comparative CT method (2−ΔΔCT method) was used to quantify the relative expression of genes, and the actin gene was selected as an internal control to normalize the data. Each experiment was performed in three technical replicates.

#### **3. Results**

#### *3.1. Physiological Responses to Heat Stress*

In order to explore the effect of heat stress on basil growth, after treating at high temperature (42/36 ◦C, light/dark, 12 h/12 h) with different times, the phenotypes of basil were analyzed. As shown in Figure 1, the plants as the control group (CK) grew well and the leaves unfurled (Figure 1a). The leaves appeared slightly curly and saggy after 1 d of heat treatment (Figure 1b), became more curly and drooping after 3 d of heat treatment (Figure 1c), and the leaves as well as the plants were severely withered to approach death after 5 d of heat treatment (Figure 1d). In addition, the heat-induced oxidation damage to basil growth was investigated. It was found that compared with the CK, the MDA content as an indicator of cell membrane damage significantly increased by 1.3-fold after 1 d of heat treatment, and increased continuously with the extension of treatment time (Figure 2a). Similarly, the ROS content including H2O2 and O2 − was increased dramatically with increasing heat treatment time (Figure 2b,c), and subsequent experiments only analyzed plants that had been heat-treated for 1 d and 3 d, as the plants were close to death after 5 d of heat treatment.

**Figure 1.** Phenotype of basil under heat stress. (**a**) Phenotype of control group (CK), (**b**) phenotype of high-temperature treatment for 1 d, (**c**) phenotype of high-temperature treatment for 3 d, (**d**) phenotype of high temperature treatment for 5 d.

**Figure 2.** Determination of oxidizing substance content. (**a**) MDA content, (**b**) H2O2 content, (**c**) O2 content. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d; H5, group with heat stress treatment for 5 d.

Moreover, after being exposed to heat stress for 1 d and 3 d, the photosynthetic capacity and photochemical activity of basil seedlings were also measured. The net photosynthetic rate (Pn) and the maximum photochemical efficiency of PSII (Fv/Fm) showed no significant difference compared with the CK after 1 d of treatment, but dramatically declined after 3 d of treatment (Figure 3a,b). Furthermore, under heat stress, the value of nonphotochemical quenching (NPQ) indicated that the change in plant heat dissipation capacity increased gradually with treatment time, while the value of photochemical quenching coefficient (qP) indicated that the openness of the PSII active center gradually decreased with time (Figure 3c,d).

**Figure 3.** Measurement of photosynthetic parameters and chlorophyll fluorescence parameters under heat stress. (**a**) Pn, the net photosynthetic rate, (**b**) Fv/Fm, the maximum photochemical efficiency of PSII, (**c**) NPQ, the nonphotochemical quenching, (**d**) qP, the photochemical quenching coefficient. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d.

#### *3.2. Metabolomics Analysis of the Leaves under Heat Stress*

As an aromatic plant, basil comprises various volatiles with important medicinal activities. To further understand the changes in volatiles under heat stress, the metabolites in basil leaves after 1 d and 3 d of heat treatment were analyzed by GC–MS technology, and a total of 222 metabolites were identified, including various nucleic acids, organic acids, peptides, and steroids (Table S2, Figure 4a). Principal component analysis (PCA) showed the internal structures of several variables related to different principal components. In the PCA score plot (Figure 4b), all biological replicates were clustered together, indicating a good correlation among core samples within groups, and also verifying the data reliability. In addition, there was an obvious distinction between the CK and the heat-stressed groups, where the latter almost clustered together and thus possessed similar metabolic profiles and phenotypes, which were significantly different from the former. Overall, the above results revealed that heat stress could significantly influence the metabolism of volatile substances in basil.

#### *3.3. Differentially Accumulated Metabolites (DAMs) Induced by Heat Stress*

Orthogonal projections to latent structures discriminant analysis (OPLS-DA) as a method of multivariate statistical analysis with monitoring pattern recognition can effectively eliminate irrelevant effects on the research, and scientifically screen different metabolites. As shown in Figure 5, distinct metabolic differentiations were observed between the CK and treatment 1 (H1), between the CK and treatment 2 (H3), and between treatment 1 (H1) and treatment 2 (H3). All Q2 values in the three models were greater than 0.5, suggesting that all models were adequately reliable. Subsequently, the fold change (FC) and variable importance in project (VIP) values of the OPLS-DA model were combined to further screen the DAMs. By setting FC ≥ 2 or ≤ 0.5 together with VIP ≥ 1.0 as thresholds for significant differences, the significantly DAMs of 29 between CK and H1 (H1 had 16 up-regulated and 13 down-regulated), 37 between CK and H3 (H3 had 17 up-regulated and 20 down-regulated), and 9 between H1 and H3 (H1 had 5 up-regulated and 4 downregulated) were determined, respectively, as presented in Table S3. Furthermore, Venn diagram analysis revealed that although no common DAMs were found in each comparison group, 18 DAMs were observed in CK-H1 and CK-H3 (Figure S1, Table S3), indicating that these metabolites may play crucial functions in response to heat stress.

**Figure 5.** The score plots of orthogonal projections to latent structures discriminant analysis (OPLS-DA) pairwise comparisons of differentially accumulated metabolites (DAMs). (**a**) Comparison groups of CK vs. H1, (**b**) comparison groups of CK vs. H3, (**c**) comparison groups of H1 vs. H3. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d.

#### *3.4. Pathway Enrichment Analysis of DAMs under Heat Stress*

Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, DAMs from each comparison group were mapped to metabolic pathways. As displayed in Figure 6, KEGG pathway enrichment analysis showed that DAMs in the CK-H1 pair were mainly enriched in alanine, aspartate, and glutamate metabolism; glycolysis or gluconeogenesis; tyrosine metabolism; amino sugar and nucleotide sugar metabolism; and another 3 pathways. DAMs in the CK-H3 pair were mainly enriched in 7 pathways including aminoacyl-tRNA biosynthesis; alanine, aspartate, and glutamate metabolism; phenylalanine, tyrosine, and tryptophan biosynthesis; and nitrogen metabolism. DAMs in the H1-H3 comparative pair were enriched in the citrate cycle (TCA cycle); pantothenate and CoA biosynthesis; and valine, leucine, and isoleucine biosynthesis. All these enriched pathways in the above were associated with the heat tolerance of basil.

**Figure 6.** Koto Encyclopedia of Genes and Genomes (KEGG) pathway classification of comparison groups. (**a**) Comparison groups of CK vs. H1, (**b**) comparison groups of CK vs. H3, (**c**) comparison groups of H1 vs. H3. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d.

#### *3.5. Transcriptomic Analysis of Basil Responses to Heat Stress*

To dissect the molecular response of basil to heat stress, the global gene expression profiles were detected by RNA-seq. A total of 9 cDNA libraries were constructed from CK, H1, and H3 groups. Sequencing analysis revealed that about 145.2 Gb of reads were generated. The clean reads were assembled de novo using Trinity software, and the redundancies in the assembled transcript sequences were removed using CD-HIT-EST software. In total, 205,848 unigenes were obtained with an average length of 419 bp and an N50 length of 977 bp. Whereafter, all sequences of the assembled unigenes were searched in the KOG, NT, NR, Uniprot, and Swiss-Prot databases and then annotated using the BLASTX program. In order to obtain the reliable profiles of gene expressions, the transcripts per million (TPM) values were used to quantify gene level. Differentially expressed genes (DEGs) were selected to annotate the differential expression between the CK and the heat-stressed groups based on the filter criteria of FC (a ratio of stress/control) values ≥ 2 and *p*-value < 0.05. In the result, a total of 19,400 DEGs were identified from the three comparison groups (WT-H1, WT-H3, and H1-H3), including 15,066, 15,445, and 192 DEGs in WT-H1, WT-H3, and H1-H3 pairs, respectively (Table S4). Among them, 11,183 DEGs were common response genes after heat stress of 1 d and 3 d, including 5437 down-regulated DEGs and 5746 up-regulated DEGs (Figure 7a,b). Interestingly, only 17 of 11,183 common DEGs showed a significant difference in the H1-H3 pair (Figure 7a). In addition, all DEGs were classified into 6 sub-groups with different numbers using Kmeans cluster analysis (Figure S2), indicating that DEGs could respond to heat stress in specific ways. All above results showed significant differences in gene expression between the CK and a heat-stressed group, while the differences among the heat-stressed groups were indistinctive.

**Figure 7.** Venn diagram analysis of differentially expressed genes (DEGs). (**a**) Venn diagram showing the overlapping and unique DEGs amongst the comparison groups, (**b**) Venn diagram showing the up-regulated and down-regulated DEGs in CK vs. H1 and CK vs. H3. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d.

#### *3.6. GO and KEGG Analysis of DEGs under Heat Stress*

The DEGs were further classified and characterized according to the functional category defined by GO and KEGG. To identify various GO terms from the significantly enriched DEGs, the function of DEGs was analyzed using TBtools software after heat stress. All 19,400 DEGs were classified into three main GO groups including molecular function, cellular component, and biological process. The 30 most enriched GO terms were shown in Figure 8, where the DEGs annotated in the molecular function group included glyceraldehyde-3-phosphate dehydrogenase (NAD+) (phosphorylating) activity, glyceraldehyde-3-phosphate dehydrogenase (NAD(P) +) (phosphorylating) activity, disordered domain specific binding, catalytic activity and hydrolase activity, and hydrolyzing O-glycosyl compounds; the DEGs in the cellular component consisted of 20 GO terms such as stromule, plastid membrane, plastid thylakoid, cell wall, chloroplast, and another 15 GO terms; and the DEGs in the biological process involved the small-molecule biosynthetic process, carbohydrate metabolic process, carbohydrate catabolic process, glucan catabolic process, and polysaccharide metabolic process.

**Figure 8.** GO enrichment analysis of all differentially expressed genes (DEGs) in basil under heat stress.

To further uncover the effect of heat stress on the biological pathway and signal transduction pathway, all DEGs were mapped to the KEGG database. As a result, the 30 most enriched KEGG pathways were identified and are shown in Figure 9, of which glyoxylate and dicarboxylate metabolism was the most dominant pathway, followed by starch and sucrose metabolism, the biosynthesis of other secondary metabolites, and metabolism. In addition, based on the Mapman results (Figure 10), most of the DEGs were significantly down-regulated in photosynthetic pathways relating to lipids, ascorbate, glutathione, starch, and sucrose, indicating that these metabolic pathways were suppressed by heat stress, while most of the DEGs were significantly up-regulated in the pathways relating to terpenes, phenylpropanoids, phenolics, and nucleotides, suggesting that these metabolic pathways were enhanced by heat stress.

**Figure 9.** KEGG enrichment analysis of all differentially expressed genes (DEGs) in basil under heat stress.

**Figure 10.** Comprehensive view of DEGs involved in multiple metabolic pathways analyzed by MapMan.

#### *3.7. The DEGs Participating in Alanine, Aspartate, and Glutamate Metabolism*

By combining transcriptome analysis with metabolome analysis, it was found that the metabolism of alanine, aspartate, and glutamate was significantly affected under heat stress. As presented in Table S3, after 1 d and 3 d of heat treatment, the L-alanine content dramatically reduced to 4.79% and 5.9% of its original content in the CK, respectively, and the oxalacetic acid content was similar to the original content after 1 d of heat treatment but decreased to 17% of the original content after 3 d of heat treatment (Table S3). A total of 84 DEGs were found to be involved in this pathway through KEGG analysis (Table S5). Alanine-glyoxylate aminotransferase (AGT) could catalyze the conversion of glyoxylate to glycine, and owing to the largest proportion of all DEGs in this pathway, 12 genes encoding AGT were significantly decreased after heat stress. Moreover, the expression levels of glutamate-glyoxylate aminotransferase (GGT1) and glutamine synthetase 2 (GS2) were also down-regulated after heat stress. On the contrary, many other genes such as glutamine-dependent asparagine synthase 1 (AS1), asparagine synthetase 3 (AS3), and carbamoyl phosphate synthetase A (carbamoyl phosphate synthetase A) were up-regulated after heat stress (Figure 11a). In addition, K-means cluster analysis showed that all the 84 DEGs could be classified into 4 sub-groups (Figure 11b), implying that they had specific expression patterns under heat stress.

**Figure 11.** DEGs related to alanine, aspartate, and glutamate metabolism. (**a**) Hierarchical clustering of expression profile of DEGs, (**b**) K-means cluster analysis of DEGs.

#### *3.8. Differentially Expressed Transcription Factors in Response to Heat Stress*

Transcription factors (TFs) play important roles in plant growth and development, hormone production, as well as the response to abiotic and biotic stress [33]. In this study, 206 differentially expressed TFs were found in the DEGs, including 139 genes being up-regulated and 67 genes being down-regulated under heat treatment (Table S6). As shown in Figure 12, all these differentially expressed TFs could be divided into 17 families, including bZIP, bHLH, WRKY, and TCP. With holding the most members in all TF families, 26 DEGs were identified to belong to the heat-shock transcription factor (HSF) gene family. Interestingly, DEGs in E2F/DP, ERF, and Trihelix families only showed an upregulated expression trend under heat treatment, while DEGs in C2H2 and SBP families only showed a down-regulated expression trend. It was also worth mentioning that most of the differentially expressed WRKYs were up-regulated except for one gene (Cluster-9139.92087). All the above results suggested that basil would adapt to heat stress by altering the expression of its own TFs.

**Figure 12.** Differentially expressed transcription factors involved in heat stress.

#### *3.9. Validation of RNA-seq Data by qRT-PCR*

To validate the reliability of the transcriptome analysis results, 20 genes related to hormone metabolism, transcription factors, and photosynthesis were selected for qRT-PCR analysis to detect their expression levels after heat stress treatments. Overall, the two analysis methods demonstrated basically consistent trends of gene expression (Figure 13), confirming that the RNA-seq data were reliable.

**Figure 13.** Validation of RNA-seq by qRT-PCR. (**a**) qRT-PCR analysis of DEGs in comparison groups of CK vs. H1, (**b**) qRT-PCR analysis of DEGs in comparison groups of CK vs. H3. CK, control group; H1, group with heat stress treatment for 1 d; H3, group with heat stress treatment for 3 d.

#### **4. Discussion**

#### *4.1. The Ability of Basil to Resist Heat Stress*

With the intensification of global warming, heat stress is becoming a severe abiotic stress to restrict the quality and yield of plants. Uncovering the mechanisms of plants in response to heat stress will facilitate the breeding of heat-resistant vegetable and crop cultivars. In this work, physiological changes combined with multi-omics analysis were used to investigate basil's response to heat stress of different times. Generally, photosynthesis is often considered as one of the most heat-sensitive physiological processes [34–36]. Heat stress commonly causes serious thermal damages to thylakoid structures, dramatically inhibits carbon metabolism and photochemical reactions, and reduces the rate of photosynthesis [11]. For instance, after being stressed at 43 ◦C for 2 h, tobacco would show an obvious decrease in net photosynthetic rate (Pn), stomatal conductance, and carboxylation efficiency (CE) of photosynthesis [37]. Heat stress would damage chloroplast membranes and thylakoid membranes, thus reducing the capacity of photosynthesis in soybean [38].

Consistent with previous studies, photosynthesis in basil was also suppressed under heat stress, as the values of Pn and Fv/Fm were severely decreased after 3 d of heat treatment, and most of the photosynthetic pathway-related genes were significantly downregulated under heat stress (Figure 3a,b and Figure 10). In addition, it is worth pointing out that the temperatures in the treatment conditions used by this study (42/36 ◦C, light/dark, 12 h/12 h) were abnormal in South China, and it was rare to have such a high temperature continuously for several days. Interestingly, the values of Pn and Fv/Fm showed no significant difference compared with CK after 1 d of heat treatment (Figure 3a,b), and the leaves just showed slightly curling (Figure 1b). Meanwhile, under heat stress, NPQ and qP displayed the trend to increase and decrease, respectively (Figure 3c,d). These results revealed the photosynthetic thermal adaptation of basil to short-term heat stress by regulating heat dissipation and the openness of the PSII active center. Previous studies have shown that heat stress often caused oxidative stress in plants [8]. The contents of H2O2, O2 −, and MDA were increased dramatically with increasing heat treatment time (Figure 2), indicating that basil was suffering from severe oxidative stress. It is worth pointing out that in *Ocimum basilicum* L. cv. 'Genovese', the content of H2O2 did not change significantly and the MDA content decreased under high temperature at elevated CO2 [9]. The difference in oxidative stress between the two varieties under high temperature may be caused by the difference in treatment conditions and genetic background, which needs further study. Despite suffering from a reduced photosynthetic capacity and severe oxidative stress, basil could survive for five days under heat stress (Figure 1d), and might live longer if watered again. Overall, basil (*Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth.) exhibits a certain ability to resist heat stress.

#### *4.2. Involvement of Key Metabolites in Response to Heat Stress*

Secondary metabolites play a variety of functions in plant physiological processes such as the response to environmental changes, innate immunity generation, growth, and development [39,40]. Metabolomics has been widely used in abiotic stress research due to its high-throughput detection of chemical substances, and has been proposed as an effective means of enhancing heat resistance [41]. In soybean, for example, an increase in sugar and nitrogen metabolism provided an obvious help in dealing with heat stress [42], and metabolomics analysis of wheat grain at the filling stage also showed that the content of carbohydrates decreased significantly under heat stress, while the contents of amino acid and biosynthetic starch increased [43]. In this study, 29 significantly DAMs were found between CK and H1, while 37 DAMs found between CK and H3 (Table S2). Venn diagram analysis revealed that 10 DAMs were specifically found in the CK-H1 comparison pair, and 14 DAMs were found in the CK-H3 pair, suggesting that these DAMs were specifically involved in response to early or long-term heat stress (Figure S1). Moreover, 18 common DAMs were observed in CK-H1 and CK-H3 comparison pairs, including 10 up-regulated DAMs such as pecazine, 4-hydroxylphenyllactic acid, beta-sitosterol, conduritol-betaepoxide, and quinic acid, and 8 DAMs including L-serine, D-glucosone, D-arabinopyranose, and L-alanine. All these key DAMs may play crucial roles in regulating heat resistance, and the detailed mechanism needs further study.

#### *4.3. Molecular Basis of Basil Resistance to Heat Stress*

Transcriptome analysis showed that a total of 19,400 DEGs were identified in the three groups (WT-H1, WT-H3, and H1-H3), and 15,066, 15,445, and 192 DEGs were found in WT-H1, WT-H3, and H1-H3, respectively. Venn analysis revealed that 11,166 common DEGs exhibiting similar expression trends were observed in the CK-H1 and CK-H3 comparison group, accounting for about 58% of all DEGs, and their expression levels altered immediately after 1 d of heat treatment, being maintained at these levels after 3 d of heat treatment; 3830 DEGs were distinctively observed in WT-H1, 4212 DEGs were observed in WT-H3, and their expressions responded specifically to heat stress time without causing significant differences in H1-H3 (Figure 7). Interestingly, it should be noted that 17 common DEGs were observed in all comparation groups, such as the heat shock protein 70 family (HSP70) and the gibberellin regulated protein and multicopper oxidase family, and their expressions levels were more closely correlated with heat treatment time. In sum, there

was a significant difference at the transcriptional level between the heat-stressed groups and the CK, while a slight difference was caused by the time of heat treatment.

All DEGs were highly enriched in 30 KEGG pathways including glyoxylate and dicarboxylate metabolism, starch and sucrose metabolism, and the biosynthesis of other secondary metabolites (Figure 9). According to the Mapman analysis results, the majority of DEGs were significantly down-regulated in the pathways related to the photosynthesis of lipids, ascorbate, glutathione, starch, and sucrose, indicating that these metabolic pathways were suppressed by heat stress; on the contrary, some metabolic pathways were strengthened by heat stress, because most DEGs involved in terpenes, phenylopropanoids, phenolics, and nucleotides were significantly up-regulated (Figure 10).

When plants are subjected to heat stress, the expression and accumulation of many heat shock proteins (HSPs) are always induced to enhance their adaptability to extreme environments [44]. In the current study, a total of 37 HSPs were found in all DEGs (Table S7). Among them, 35 HSPs were significantly down-regulated after heat treatment of 1 d and 3 d, probably because the expressions of HSPs are always activated immediately after a short-term heat treatment and would decrease at a later time [45]. In comparison with the CK, however, a small heat shock protein HSP18.2 (Cluster-9139.239298), respectively, increased by 3.36-fold and 5.72-fold after 1 d and 3 d of heat treatment, and the mitochondrial localized HSP60-3A (Cluster-9139.246542) also increased by 1.96-and 3.51-fold after 1 d and 3 d of heat stress, respectively, implying that the two HSPs could be involved in the regulation of long-term heat stress. In addition, some research studies have demonstrated that transcription factors could play crucial roles in heat stress. For example, a heat stress transcription factor *LlHsfA4* enhanced thermotolerance via regulating ROS metabolism in Lilies [46], and *AtMYB30* regulated the responses of *Arabidopsis* to oxidative and heat stress through calcium signaling, partially mediated by ANNEXIN (ANN) genes [47]. Li et al. (2011) reported that the constitutive expression of *WRKY25*, *WRKY26*, or *WRKY33* enhanced the heat stress resistance of *Arabidopsis* [48]. In this work, a total of 206 differentially expressed transcription factors (TFs) were identified in all DEGs, and they belonged to various families (Table S6). Among them, 58 TFs were significant differentially expressed after 1 d of heat stress, and 59 TFs were significant differentially expressed after 3 d of heat stress, indicating that these TFs were specifically involved in the early-stage or long-term response to high temperature. The other common 128 differentially expressed TFs were observed in the CK-H1 and CK-H3 comparison group, but not in the H1-H3 group, meaning that these TFs may be rapidly induced by heat stress and maintained at a certain level of expression to regulate other genes expression. Overall, all these DEGs could contribute to heat resistance, and their potential regulatory mechanisms should be explored in the future.

#### **5. Conclusions**

In conclusion, the potential mechanisms of basil (*Ocimum basilicum* Linn. var. *pilosum* (Willd.) Benth.) in response to heat stress were explored by combining physiological changes with multi-omics analysis. Heat stress disturbed the development of basil by causing severe oxidative damage and inhibiting photosynthesis. Metabolome and transcriptome analyses identified a total of 51 DAMs and 19,400 DEGs in basil after heat stress, and they could be enriched in diverse metabolic pathways such as glyoxylate and dicarboxylate metabolism, glycolysis or gluconeogenesis, and aminoacyl-tRNA biosynthesis. This paper has provided meaningful information to further understand the mechanism of the plant response to heat stress, which may also be useful in crop breeding.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12061434/s1, Figure S1: Venn diagram analysis of DAMs amongst the comparison groups; Figure S2: K-means cluster analysis of all DEGs; Table S1: Primer sequences in this study; Table S2: All metabolites identified in samples: Table S3: DAMs in each comparison group: Table S4: DEGs in each comparison group; Table S5: DEGs in alanine, aspartate, and glutamate metabolism; Table S6: Differentially expressed transcription factors under heat stress; Table S7: Differentially expressed HSPs under heat stress.

**Author Contributions:** Q.W., A.Q. and L.Y. conceived the ideas and designed the methodology; L.Q., C.L. and J.W. collected the data; L.Q., C.L. and D.L. analyzed the data; Q.W. and A.Q. led the writing of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Department of Science and Technology of Ningbo (2021IJCGY020198), and the Natural Science Foundation of Zhejiang (LQ19C140002).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Any data and codes used in this study are available upon request from the corresponding author.

**Acknowledgments:** We are especially grateful to Yingxian Zhao for his guidance on the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

