Integrated transcriptome and metabolome revealed the drought responsive metabolic pathways in Oriental Lily (Lilium L.)

Test materials and drought treatment

Oriental Lily was planted in the garden laboratory of Fujian Forestry Vocational Technical College. The bulbs (average mass 65 g) were planted in plastic pots (15 cm (H) × 10 cm (inner diameter)), containing three small holes in the pelvic floor to ensure ventilation of the underlying soil. Each pot contained only one seed ball. The volume ratio of cultivated substrate was 3:1 (coconut brick and nutrient soil, pH = 6.15). Each pot contained approximately 350 g substrate.

A total of 60 lily seedlings with basically the same growth and no pests and diseases were selected. After 7 days of culture, they were placed in the artificial climate box to slow the seedlings. The ambient conditions were set as: day temperature (25 ± 1) °C, night temperature (18 ± 1) °C, light intensity 12,000 Lx, photocycle 12 h/12 h (day/night), relative humidity 60%.

The experiment set up two water treatments, each treatment had five replicates: the control group had normal water supply (ZC), and the soil moisture content was 70–80% of the field water holding capacity; in the drought stress group (GH), the water control was 40–50% of the field water holding capacity. Water loss was measured daily at 5 p.m. by weighing and water control, while portable TDR (TZS-PHW, China) was used to measure soil volumetric water content.

Paraffin section preparation

The fresh tissue was fixed with fixed liquid for more than 24 h. We removed the tissue from the fixed liquid and trimmed the target tissue with a scalpel in the ventilation cupboard and then transferred to dehydration box. The dehydration box was shifted to the dehydrator in order to dehydrate with gradient alcohol. The rest of the procedure was performed following the method of Yuan et al. (2023).

Observation of morphological indicators and observation of leaf anatomy

The growth indexes of each lily plant were recorded at 10, 15, 20, 25 and 30 days after treatment, and the leaves of the same leaf position were measured. The specific morphological indexes include leaf length, leaf width and leaf area using a centimeter-scale ruler. The thickness of blade, palisade tissue and spongy tissue were measured by Motic Image Plus 2.0 measuring software, and each structural parameter was the average value measured within 20 field of view.

(1) $$CTR={\displaystyle \frac{\mathrm{T}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}}{Thicknessofleaf}\times 100\mathrm{\%}}$$

(2) $SR={\displaystyle \frac{\mathrm{T}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{g}\mathrm{y}\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}}{Thicknessofleaf}\times 100\mathrm{\%}.}$

Sample preparation for metabolome analysis

The vacuum freeze-dried method was used for biological samples to grind the samples (MM 400, Retsch); 100 mg of powder was weighed and dissolved in 1.0 mL of extract and refrigerated overnight at 4 °C with vertexing 3 times. The sample underwent centrifugation for 10 min at 12,000 rpm, the supernatant was aspirated and filtered with a 0.22 μm microporous membrane for liquid chromatography tandem mass spectrometry (LC-MS/MS) analysis.

Data processing and multivariate analysis

The Proteowizard software package was used to convert the original mass spectrometry file into mzXML file format (Smith et al., 2006), and the RXCMS software package was used to analyze the chromatographic peaks. The obtained qualitative and quantitative results of metabolites. Dos Santos & Canuto (2023) were saved for further analysis. The identified metabolites were annotated using HMDB (Wishart et al., 2022), massbank (Horai et al., 2010), LipidMaps (Sud et al., 2007), mzcloud (Abdelrazig et al., 2020) and KEGG databases, (Kanehisa et al., 2022). The parameter was set to ppm <30 ppm (Asif et al., 2023).

The R software package Ropls was used for principal component analysis of normalized data (Thevenot et al., 2015). On the basis of orthogonal partial least squares discriminant analysis (OPLS-DA), score plot, load plot and S-plot were plotted to show the difference of metabolite composition among samples. The p value was calculated according to the statistical test, the projected importance of variable (VIP) was calculated by OPLS-DA dimensionality reduction method, and the component difference multiple was calculated by fold change. The differential metabolites were identified with p < 0.05 and VIP > 1 as the screening conditions for differential compounds.

Pathway analysis

The Metabolomic Analyst software package was used to conduct functional pathway enrichment and topological analysis of the screened differential metabolites.

The enrichment results were taken as the unit of KEGG pathway, and hypergeometric test was applied to obtain the p-value of pathway enrichment. With p-value ≤ 0.05 as the threshold.

(3) $$p-\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}={\sum }_{i=m}^M{\displaystyle \frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}=1-{\sum }_{i=0}^{m-1}{\displaystyle \frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}.}}$$

Library construction and sequencing for transcriptome analysis

Total RNA was extracted, and the Oligo (dT) magnetic beads were used to enrich the mRNA with polyA structure. Using RNA as a template, the first strand of cDNA was synthesized using 6-base random primers and reverse transcriptase, and the second strand cDNA was synthesized using the first strand cDNA as template. Then, the constructed library was inspected by Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), and the total concentration and effective concentration were detected. The libraries were put forward for paired-end (PE) sequencing using next-generation sequencing (NGS) based on Illumina sequencing platform (Eid et al., 2009). The longest transcript was selected as Unigene, and further used for GO, KEGG, eggNOG, SwissProt, Pfam annotation, ORF prediction, SSR prediction (Su et al., 2021).

Quantitative real-time PCR

Total RNA was extracted using the RNA Easy Fast Plant Tissue Kit (Tiangen, Beijing, China) according to the manufacturer’s instructions. The expression were analyzed by qRT-PCR using the SYBR green PCR mix (Vazyme Biotech Co., Ltd, Beijing, China). The primer details are available for TRINITY_DN2608 (F:TGAACACCTCTTGTCGAGCCR:CCGTCAACCTGACCGTAACA). Actin expression was used as an internal control. (F:ATTCAGATGCCCAGAAGTCTTGTTR:GAAACATTTTCTGTGAACGATTCCT). The relative gene expression was calculated by 2^−ΔCT (Livak & Schmittgen, 2001).

Arabidopsis transformation

The open reading frame of the cloned TRINITY_ DN2608 gene was inserted into the plant expression vector pCambia 1305. The restriction enzymes used in this experiment are EcoRI and NotI. Agrobacterium vectors were inoculated into LB liquid medium containing 25 mg·L^-1Rif and 50 mg·L^-1Kana resistance for activation culture for 2 days. The culture conditions were set as 28 °C and 200 rpm. According to the ratio of bacteria solution: LB medium = 1:50 (V/V), that is, 1 mL of the above bacteria solution was absorbed and added into LB liquid medium containing the corresponding resistance at 28 °C, overnight culture at 200 rpm until the OD600 is between 0.6–1.0. On the second day, the mixed transformation vector was placed in a 50 mL centrifuge tube, and the inflorescences of the white flower buds were immersed in the bacterial solution for about 30 S. After transformation, the plants were shaded with black plastic bags and placed in a article box for 24 h of moisture. At the end of shading, the material was placed under normal conditions for growth, and the inflorescences were re-treated with the above transformation method after 1 week. Finally, the material was placed under normal conditions to grow until the seeds matured (Jia et al., 2023).

Phenotypic analysis of lily under drought stress

The drought-stressed lily plants significantly restricted the growth compared to that of control (CK) (Figs. 1A and 1B). Under drought stress, the length, width and area of lily leaves were measured (Figs. 1C–1F). A significant decrease in the length and width was observed following drought stress. Noticeable difference in the length and width of CK and drought-stressed lily plants was recorded at the 15^th day. For the leaf area, the significant difference between CK and drought-stressed lily plants was observed at the 20^th day. Paraffin section experiment was used to investigate the effect of drought stress on the anatomical structure of lily leaf (Figs. 1G–1J). That thickness of fence tissue was about 54.49 µm under drought stress, which was approximately 12.87% lower than the control. The spongy tissue with a thickness value of 104.55 µm under drought stress was 11.41% lower than the control. Meanwhile, the drought-stressed leaf thickness was recorded at a value of 161.14 µm, which was 10.05% lower than the control. Following the morphological response of lily to drought stress, samples were taken from the control and drought treated plants for transcriptomic and metabolomic analysis.

Metabolomic analysis of Oriental Lily

The primary screening yielded a total of 37,285 metabolites, among which 3,057 were down-regulated and only 688 were up-regulated, while the remaining were insignificant (Fig. 2A). We also used volcano plot to view and analyze the differential expression of metabolites (Fig. 2B). The results showed metabolites M130T140 (mz: 130.1224), M538T961 (mz: 538.1734), M160T144 (mz: 160.1335), M134T49 (mz: 134.045) and M292T977 (mz: 291.9267) were discovered with most credible P value. According to the mass/charge ratio and p value of metabolites, the scatter plot was drawn to clearly see the distribution of different substances in the samples (Fig. 2C). It revealed that M130T140, M134T49, M160T144, M292T977 and M215T564 were significantly down-regulated whereas M433T331, M538T961, M537T961, M523T332 and M539T96 showed significant up-regulated expression.

The obtained metabolites were identified by using the Human Metabolome Database (HMDB) (http://www.hmdb.ca), massbank (http://www.massbank.jp/), LipidMaps (http://www.lipidmaps.org), and mzcloud (https://www.mzcloud.org) via MS/MS model. Finally, 39 metabolites were obtained, among which two displayed downregulated while the other 37 were upregulated. These prominent upregulated metabolites including picolinic acid, L-Gulose, D-sorbose, UDP-D-Xylose and alpha-D-glucose. Meanwhile, an array of downregulated metabolites were also recorded such as beta-carotene and glycerophosphocholine (Figs. 3A and 3B). Correlations between different metabolites were calculated to reveal the synergy (Fig. 3C). The correlation between glycerophosphocholine and ancymidol displayed significant negative correlation, whereas positive correlation was found amongst numerous others.

Transcriptome analysis of Oriental Lily

The regulation of metabolites are often governed by different transcription factors (TFs). In order to understand the molecular mechanism of these metabolites in response to drought stress, the transcriptome analysis of the control sample and drought-treated sample was conducted. All unigenes were functionally annotated, mainly in databases such as NR (NCBI non-redundant protein sequences), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genome), eggnog (evolutionary genealogy of genes: Non-supervised Orthologous Groups), Swiss-Prot and Pfam.

DESeq software was used to analyze the expression of differential genes, and the screening criteria was set as |log2FoldChange| > 1 and p-value < 0.05. Compared with the control, a number of 720 differentially expressed genes (DEGs) were mapped, among which 513 were down-regulated and 207 were up-regulated (Figs. 4A–4C). KEGG enrichment analysis of DEGs (Figs. 4D and 4E) showed that up-regulated genes were mainly enriched in plant hormone signal transduction, MAPK signaling pathway, glyoxylate and dicarboxylate metabolism. On the other hand, down-regulated genes were also enriched in plant hormone signal transduction, and other different pathways including tyrosine metabolism, ABC transporters, pentose phosphate pathway, glycolysis, gluconeogenesis, RNA polymerase, fructose and mannose metabolism and peroxisome.

Integrated analysis of transcriptome and metabolome data

To determine the correlation between DMs and DEGs, a correlation analysis was conducted based on the Pearson correlation coefficient (differential genes (p < 0.05) and differential metabolites (VIP > 1, p < 0.05)). Results with a Pearson correlation coefficient greater than 0.8 were selected. The following table shows the correlation analysis between DMs and DEGs, and it can be seen that DMs and DEGs are more negatively correlated (Fig. 5A).

A bar chart was drawn to display the significance of simultaneously enriched pathways. The results showed that the metabolic pathways potentially involved in lily drought stress response were largely associated to carbohydrate metabolism (galactose, starch and sucrose, glycolysis and gluconeogenesis) (Fig. 5B). At the same time, we selected transcripts with differential metabolites (VIP > 1, p < 0.05) and p < 0.05. Based on the metabolite information in the KEGG database, we extracted transcripts corresponding to the relevant enzymes to obtain the corresponding relationship between the two and combined them with multiple differences to display plotting (Fig. 5C). Finally, we acquired nine pairs of differentially expressed metabolites and the corresponding differential transcripts (p < 0.05). Interestingly, gene TRINITY_DN2608 showed higher association with a certain metabolite and was selected for further work. The online Blast platform analysis yielded that TRINTIY_DN2608 encode a galactose protein, mapped in the galactose pathway. Based on these findings, we put forward TRINITY_DN2608 as a candidate gene for functional characterization.

Identification and phenotypic observation of TRINITY_DN2608 gene in Arabidopsis thaliana

The TRINITY_DN2608 gene was heterogeneically transformed into Arabidopsis thaliana by using the floral dip method. The T3 generation seeds were obtained from the transgenic Arabidopsis lines carrying the TRINITY_DN2608 gene. The seeds were disinfected with alcohol and then spread on kanamycin plate for screening and cultivation. After 14 days, Arabidopsis thaliana with two to four true leaves was transplanted (Fig. 6B). Leaf DNA of positive strains was extracted, and gene specific PCR primer was used to identify whether the gene was transferred into Arabidopsis thaliana (Fig. 6A). The transgenic plants were allowed to grow in normal conditions for 3 weeks. Following that, drought stress treatment was carried out to analyze the response of TRINITY_DN2608 gene. In overexpressed Arabidopsis thaliana plants, the expression level of TRINITY_DN2608 was 3–4 times that of the wild type (WT) (Fig. 6C). No obvious phenotypic difference was recorded under normal condition. However, under drought treatment, the WT plants maintained green leaves, whereas the leaves of TRINITY_DN2608-OE showed obvious wilting phenotype (Fig. 6D and 6E). In addition, the water loss from TRINITY_DN2608-OE was higher than that of the wild type between 0 and 12 h (Fig. 6F). These results suggest that allogeneic expression of TRINITY_DN2608 gene can enhance the sensitivity of Arabidopsis thaliana to drought.