Transcriptome analysis revealed the molecular mechanism by which Chrysanthemum naktongense roots respond to low-concentration lanthanum treatment

Lanthanum (La) is the second most abundant rare earth element. Owing to improper mining of rare earth minerals and improper disposal of rare earth waste, lanthanum is overused and overenriched in soils^[1]. The Inner Mongolia region of China is one of the primary sources of light rare earth elements^[2]. In recent years, extensive mining activities have led to a significant enrichment of rare earth elements, including lanthanum^[3]. Many studies have shown that lanthanum at low concentrations can promote the growth of plants^[4]. A study showed that lanthanum at low concentrations increased the root activity of bahia grass (Paspalum notatum) and improved respiration and photosynthesis in bahia grass, thus promoting root growth^[5]. However, when the concentration of lanthanum in the soil is excessively high, growth and development in plants is inhibited^[6], and an increase in soil lanthanum concentration decreases the chlorophyll level and photosynthetic capacity of rice^[7]. Another study showed that the stem length of Medicago sativa L. plants treated with 20 mM La was significantly shorter than that of control plants^[8].

Wild Chrysanthemum is a perennial herb with both ornamental and medicinal value and strong stress resistance^[9]. Chrysanthemum naktongense, a type of wild chrysanthemum with high biomass, is widely distributed in Northeast China and has strong stress resistance. Previous studies have shown that the content of rare earth elements in plant roots is greater than that in the aboveground parts. For example, Hordeum vulgare^[10], and tea trees^[11] have higher contents of rare earth elements in their roots than in their leaves and stems. In addition, our previous work confirmed that the roots of Chrysanthemum naktongense are strongly enriched with La. However, few studies have investigated the molecular mechanism by which plant roots respond to lanthanum treatment. An in-depth exploration of the differential gene expression patterns and response pathways of Chrysanthemum naktongense roots under lanthanum treatment can help in the cultivation of lanthanum-enriching varieties of this plant. Lanthanum at excessively high concentrations affects the growth and development of Chrysanthemum naktongense. Plant hormones are involved in the intrinsic mechanisms by which plants respond to abiotic stress, and they usually play an important role in regulating plant physiological and biochemical mechanisms by stimulating signal transduction pathways^[12,13]. Studies have shown that endogenous auxin can increase the resistance of AUX1 mutant seedlings to arsenic and that auxin translocation and signaling mechanisms play a role in mediating plant resistance to heavy metal stress and in root development^[14]. Dai et al. reported that exogenous jasmonic acid decreased Cu accumulation in alfalfa and improved antioxidant system activity^[15]. Kocaman reported that the application of exogenous leaf ABA (μmol·L⁻¹ 40) limited the transport of Cd, Cr, Hg, and Sn to strawberry leaves^[16].

Heat shock protein 70 (Hsp70) exhibits a high degree of sequence conservation and is one of the most extensively studied heat shock proteins^[17]. It plays a crucial role in plants' resistance to abiotic stresses such as heat, salinity, and drought.^[18] Studies have indicated that plant hormones can regulate the expression of Hsp70. The promoter regions of the apple (Malus pumila) Hsp70 gene family members contain a plethora of plant hormone and stress-responsive elements; under treatment with methyl jasmonate (MeJA), gibberellin (GA₃), and abscisic acid (ABA), the expression levels of MdHsp70-23, MdHsp70-26, and MdHsp70-45 are the highest^[19]. In potato (Solanum tuberosum), StHsp70-7 exhibits increased expression under ABA and gibberellin (GA₃) treatment, while StHsp70-17 shows upregulated expression in response to ABA, GA₃, and salicylic acid (SA) treatments^[18]. In tobacco, ABA treatment induces NtHsp70-1, enhancing the plant's stress tolerance^[20]. However, the functions of Hsp70 in many plants, including Chrysanthemum naktongense, remain to be determined.

In recent years, high-throughput sequencing has been widely used in molecular biological research to identify key differentially expressed genes (DEGs) related to plant response to abiotic stress^[21]. The stressors include salt^[22], drought^[23], low temperatures^[24], high temperatures^[25], and heavy metals^[26]. Lanthanum is a heavy metal^[27], and studies have explored the molecular response mechanisms of plants treated with lanthanum through transcriptome analysis. For example, lanthanum at a low concentration improved the root activities of bahia grass^[5]. A combination of transcriptome and ionomic analyses revealed a fundamental pattern in Arabidopsis in response to the uptake and toxicity of the rare earth elements lanthanum and ytterbium^[28].

To study the gene expression of the root system of Chrysanthemum naktongense treated with lanthanum, identify the regulatory genes associated with lanthanum treatment, and propose a response mechanism, we conducted transcriptome sequencing via an Illumina high-throughput sequencing platform. This is the first report on the transcriptome sequencing of Chrysanthemum naktongense, and the results provide insights into the gene expression patterns and pathways of Chrysanthemum naktongense in response to lanthanum treatment.

The seeds of Chrysanthemum naktongense used in the test were collected from Chunkun Mountain (110° E, 41° N), Guyang County, Baotou city, Inner Mongolia, China. The seeds were soaked and disinfected with 5% sodium hypochlorite for 10 min, washed with distilled water 5−6 times until clean, and then placed evenly on two layers of infiltrated filter paper in a glass Petri dish. After germination, the plants were transferred to nutrient-rich soil that was not contaminated with heavy metals for culture. After two months (seedling to the 7 true leaf stage), plants exhibiting consistent and healthy growth were selected for lanthanum treatment. LaCl₃, a pure reagent, was employed as the lanthanum source. The samples were divided into three groups according to treatment duration: 0 h (control group CK-0), 24 h (TR-24), and 72 h (TR-72), and three biological replicates were included in both the treatment groups and the control group. The concentration of La³⁺ was set to 100 mg·L⁻¹; this concentration was selected based on early screening by our research group (under the premise of maintaining normal growth in Chrysanthemum naktongense and yielding the highest antioxidant system index). At the same time, the whole root system was sampled, frozen in liquid nitrogen, and stored at −80 °C for later use.

RNA extraction and cDNA library construction

A total of about 1 g of Chrysanthemum naktongense root was preserved in liquid nitrogen and ground into a powder. Total RNA was extracted with a TRIzol reagent. The concentration and purity of total RNA were determined with a NanoDrop2000 instrument. The integrity of total RNA was determined via agarose gel electrophoresis. RIN values were obtained via an Agilent RNA 6000 Nano Kit. The total RNA that met standards was used to construct a library of nine samples. High-throughput sequencing was conducted at Beijing Kaitai Ming Mirror Technology Co., Ltd (Beijing, China).

Transcriptome assembly

First, the original reads sequenced via Illumina (Illumina Novaseq 6000) were treated to remove reads with contamination, high unknown base N content (> 5%), and low quality. High-quality clean reads were subsequently obtained. Trinity^[29] was then used for de novo assembly, Kallisto^[30] software was used to estimate the expression levels of the assembled transcripts, and an expression level FPKM below 1 was selected as the threshold. Low-expression transcripts were filtered out, and only transcripts with the highest expression levels of a gene were retained, representing the genetic signature. This reduced false positive results due to assembly. Then, CD-hit-est^[31] was used for cluster analysis of the above results, and only one gene with high similarity was retained. The main parameters for CD-hit-est were -c0.98-p1-d0-b3-T10. The GC content and the Q20 and Q30 data were also evaluated.

DEGs screening and qRT‒PCR verification

The transcriptome data of the root system of plants treated with lanthanum for 24 h (TR-24), and 72 h (TR-72) were compared with those of the control group (CK-0). DESeq2 software was used for differential expression gene analysis, and |logFC| (FC: fold change) > 1 and p < 0.05 were selected as the thresholds for differential gene screening^[32]. The R cluster Profiler package^[33] was used for functional annotation of Clusters of Orthologous Groups (COG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genes and Genomes (KEGG) enrichment analyses were carried out for differentially expressed unigenes. The real-time fluorescence-based quantitative differential gene primers used are shown in Supplementary Table S1. The annotation results were visualized on the online cloud platform of Beijing Kaitai Spiegel Biological Co., Ltd. (https://kaitai.cloud/all).

The identification and bioinformatics analysis of Hsp70 gene family members in Chrysanthemum naktongense

Identification and bioinformatics analysis of Hsp70 family members in Chrysanthemum naktongense were conducted using the online software InterPro (www.ebi.ac.uk/interpro) to search for and download the conserved Hsp70 protein domain (accession number: PF00012). Local software HMMER 3.0 was employed to search for potential Hsp70 protein sequences within the Chrysanthemum naktongense transcriptome protein sequences. Additionally, online tools SMART (http://smart.embl-heidelberg.de), Pfam (http://pfam.xfam.org), and NCBI-CDD (www.ncbi.nlm.nih.gov/cdd) were utilized to identify conserved domains in the Chrysanthemum naktongense Hsp70 protein sequences and to eliminate redundant sequences. The physicochemical properties of the aforementioned protein sequences were predicted using the online software Expasy (https://web.expasy.org/protparam). Subcellular localization of the proteins was predicted using the online software WoLF PSORT (https://psort.hgc.jp). Hsp70 gene sequences from Arabidopsis and rice were downloaded from the Arabidopsis Genome Initiative (www.arabidopsis.org), and the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu), respectively. Phylogenetic trees were constructed using MEGA 7.0 with the Neighbor-Joining (NJ) method, and a bootstrap analysis was performed with 1,000 replicates using pairwise gap deletion mode. The phylogenetic tree results were visualized using the online software iTOL (https://itol.embl.de). Finally, the conserved motifs of the Hsp70 family members were analyzed using the online software MEME (http://meme-suite.org), and the results were visualized using the TBtools software (version: 2.1.19).

Quantitative RT-PCR verification

The E.Z.N.A Plant RNA Kit (OMEGA, USA) was employed to extract total root RNA from whole roots subjected to treatments for 0, 24, and 72 h. The extracted RNA was subsequently reverse transcribed into cDNA using the ReverTra Ace qPCR RT Master Mix (TOYOBO, Japan). SYBR Green Realtime PCR Master Mix (TOYOBO, Japan) was analyzed, using 7 μL ddH₂O, 1 μL cDNA, 2 μL of upstream and downstream primers, and 10 μL of SYBR enzyme as the reaction system. The relative gene expression was analyzed by the 2^−ΔΔCᴛ method. Gene-specific primers were designed utilizing the online NCBI Primer-BLAST tool (www.ncbi.nlm.nih.gov/tools/primer-blast, Supplementary Table S1).

Data analysis

Microsoft Excel (2019) software was used to conduct a one-way analysis of variance and produce line charts. Duncan's multiple intervals were used to analyse significant differences, and differences were considered statistically significant at p < 0.05.

As shown in Table 1, a total of 56.48 Gb of clean bases were obtained from nine libraries of roots exposed to lanthanum for different durations. The data availability rates ranged from 97.86% to 98.52%, the percentage of Q20 bases was 98.6% or greater, the percentage of Q30 bases was 95.65% or greater, and the GC content was 43.20% or greater, which suggested that the sequencing data were reliable and could be used for subsequent analysis. A total of 62,125 unigenes were obtained after sequence assembly. By comparing the sequences with those in the COG, GO, and KEGG databases, functional annotation of unigenes was carried out, and 43,841 annotated unigenes were ultimately obtained.

COG annotation

To preliminarily explore the response mechanism of the Chrysanthemum naktongense root system to lanthanum treatment, we performed COG annotation of the sequencing results. As shown in Fig. 1a, a total of 44,284 unigenes were annotated in the COG database, but most of them had unknown functions although they ranked first in the COG prediction, indicating that many sequences in the root system of Chrysanthemum naktongense unigenes have not been effectively verified and studied. Among the unigenes, signal transduction mechanisms (T), posttranslational modification, protein turnover, chaperone posttranslational modification, protein turnover, chaperone (O), translation, ribosomal structure and biogenesis (J), and transcription (K) were found at positions 2-5, suggesting that the intrinsic mechanism of root response to lanthanum treatment may involve the basic processes of growth and development and signal transduction.

Figure 1. 

COG annotation results of DEGs and DEGs in the roots of Chrysanthemum naktongense under lanthanum stress for different durations. (a) COG annotation results of differentially expressed genes. All the predicted genes were consistent with the COG database and were functionally grouped into 24 molecular families. The X-axis represents the number of differentially expressed genes assigned to the corresponding COG class, and the Y-axis represents the COG class to which the differentially expressed genes are annotated. (b)−(d) shows the volcano map of differentially expressed genes in TR-24 vs CK-0, TR-72 vs CK-0, and TR-72 vs TR-24. (e), (f) show Venn diagrams of upregulated and downregulated DEGs in TR-24 vs CK-0, TR-72 vs CK-0, and TR-72 vs TR-24, respectively. The threshold for differentially expressed gene screening was as follows: (|logFC| > 1, p < 0.05).

Screening and identification of differentially expressed genes

After 24 h of lanthanum treatment (TR-24), we identified 3,266 differentially expressed genes (DEGs) in the root system of prickly lettuce compared to the control group without lanthanum treatment (CK-0). Among these, 1,551 genes were upregulated and 1,715 were downregulated (Fig. 1b). Following 72 h of lanthanum treatment (TR-72), we observed 5,993 DEGs, with 3,146 upregulated and 2,847 downregulated genes (Fig. 1c). A comparison between TR-72 and TR-24 revealed a total of 1,456 DEGs; among them, 866 were upregulated and 590 were downregulated (Fig. 1d). Utilizing Venn diagrams for further analysis (Fig. 1e & f), we found that there were: 977 genes upregulated in both TR-24 vs CK-0 and TR-72 vs CK-0, 1,261 genes downregulated in both comparisons, the comparisons between TR-72 vs CK-0 and TR-72 vs TR-24 revealed a total of 596 co-upregulated genes and 358 co-downregulated genes. These results indicated an increase in specific differentially expressed genes as well as co-expressed ones following lanthanum treatment over time.

DEGs GO enrichment analysis and KEGG enrichment analysis

To further understand the functional pathways related to the above DEGs, GO enrichment analysis of the top 20 upregulated DEGs in each treatment group was performed (Fig. 2a−c). Compared with those in the 0 h treatment group (CK-0), 125 DEGs related to BPs were upregulated in the 24 h treatment group (TR-24 vs CK-0, Fig. 2a). The most significantly enriched terms (p < 0.05) were related to secondary metabolic processes (GO:0019748, 19 DEGs), response to jasmonic acid (GO:0009753, 16 DEGs), secondary metabolite biosynthetic processes (GO:0044550, 14 DEGs), and hormone metabolic processes (GO:0042445, 13 DEGs). Twenty-five upregulated DEGs related to molecular function (MF), including O-methyltransferase activity (GO:0008171, seven DEGs), palmitoyl-(protein) hydrolase activity (GO:0008474, three DEGs), were enriched. Compared with CK-0, 272 upregulated DEGs related to the biological process (BP) were enriched under 72 h lanthanum treatment (TR-72 vs CK-0, Fig. 2b). The most abundant terms were related to hormone metabolic processes (GO:0042445, 31 DEGs), response to nitrate (GO:0010167, 26 DEGs), response to chitin (GO:0010200, 25 DEGs), etc. There were 123 annotations related to molecular function (MF), molecular transducer activity (GO:0060089, 41 DEGs), signaling receptor activity (GO:0038023, 38 DEGs), anion transmembrane transporter activity (GO:0008509, 33 DEGs), etc. In addition, 33 upregulated DEGs related to cellular components (CC) and apoplasts (GO:0048046) were significantly enriched. In addition, in TR-72 vs TR-24 (Fig. 2c), 123 upregulated DEGs were annotated to the biological process (BP) category. The most prominent enrichment was found in external encapsulating structure organization (GO:0045229, 22 DEGs). In addition, 71 DEGs were annotated to the cellular component (CC) category, and the most significant enrichment was found in cell wall organization (GO:0071555, 21 DEGs), etc. Twenty upregulated DEGs were significantly enriched in molecular functions (MF) and were related to the plastid thylakoid lumen (GO:0031978, four DEGs), transferase activity, and transfer of hexosyl groups (GO:0016758, 16 DEGs).

Figure 2. 

Column diagram of the top 20 enriched GO terms and bubble diagram of the enriched KEGG pathways of the upregulated DEGs in each treatment group. (a)−(c) Show the top 20 GO enrichment terms in the sequences of TR-24 vs CK-0, TR-72 vs CK-0, and TR-72 vs TR-24. BP refers to biological process, CC refers to cellular component, and MF refers to molecular function. The vertical axis shows the name of the GO enrichment term, and the horizontal axis displays the number of genes enriched in each GO term. The length of the rectangle represents the number of DEGs associated with each enrichment term, and the colour depth represents the significance of each enrichment term. (d)−(f) Shows the bubbles of the TR-24 vs CK-0, TR-72 vs CK-0, and TR-72 vs TR-24 upregulated DEGs in the top 20 enriched KEGG terms. The vertical axis represents the name of the enriched KEGG pathway, and the horizontal axis represents the GeneRatio of each enriched pathway. The size of the points represents the number of DEGs in each enriched pathway, and the depth of the colour of the points represents the significance of each enriched term.

Through GO enrichment analysis, it was found that the GO enrichment terms of the top 20 DEGs in each treatment group that were upregulated (Fig. 2a−c) were mainly related to biological processes, in which the secondary metabolite generation pathway, hormone metabolism, and plant cell wall-related biological processes played crucial roles in the response of the C. naktongense root system to lanthanum treatment.

To further explore which biological pathways played important roles in the response to lanthanum treatment, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs under the two treatment durations. The results revealed that (Supplementary Tables S2, S3) in TR-24 vs CK-0, 1,243 DEGs were annotated to 282 pathways. In TR-72 vs CK-0, 2,880 DEGs were annotated to 332 pathways. TR-24 vs CK-0 showed significant enrichment in 37 KEGG pathways (p < 0.05), TR-72 vs CK-0 exhibited significant enrichment in 58 KEGG pathways (p < 0.05), and all of them were dominated by upregulated DEGs.

Compared with CK-0, TR-24 had 208 significantly upregulated DEGs (p < 0.05) enriched in 29 KEGG pathways. The top 20 enriched pathways of the upregulated DEGs (Fig. 2d) included the oestrogen signaling pathway (ko04915, 19 DEGs). The multiple DEGs were related to the longevity-regulating pathways (ko04213, 14 DEGs), phenylpropanoid biosynthesis (ko00940, 13 DEGs), and plant hormone signal transduction (ko04075, 13 DEGs). In TR-72 vs CK-0, 556 upregulated DEGs were significantly enriched (p < 0.05) in 42 KEGG pathways, and the top 20 (Fig. 2e) included. Plant hormone signal transduction (ko04075, 38 DEGs), longevity regulating pathway - multiple species (ko04213, 24 DEGs), phenylpropanoid biosynthesis (ko00940, 19 DEGs), etc. TR-72 vs TR-24 had 134 upregulated enriched DEG pathways related to 21 KEGG pathways, the top 20 (Fig. 2f) of which included phenylpropanoid biosynthesis (ko00940, nine DEGs), plant hormone signal transduction (ko04075, nine DEGs). These results indicated that, under lanthanum treatment, the main pathways that significantly changed in the root system were the oestrogen signaling pathway, phenylpropanoid biosynthesis, plant hormone signal transduction, and other pathways. Importantly, the GO and KEGG enrichment results revealed that upregulated DEGs were involved in multiple pathways related to the response to lanthanum stress, among which plant hormone-related pathways were the most significant and highly enriched pathways.

Effects of lanthanum treatment on plant hormone signal transduction pathways in the roots of Chrysanthemum naktongense

After lanthanum treatment for 24 and 72 h, the number of upregulated DEGs enriched in the plant hormone signal transduction pathway was greater than that enriched in the other KEGG pathways, and the difference was significant. To investigate the effects of lanthanum treatment on the hormone signaling pathway in the root system of Chrysanthemum naktongense, we analyzed the expression patterns of genes involved in this pathway. The results (Supplementary Tables S2, S3) revealed that in this pathway, 13 and 38 samples from TR-24 vs CK-0 and TR-72 vs CK-0, respectively, presented upregulated DEGs and 19 samples from TR-72 vs CK-0 presented downregulated DEGs and were enriched in this pathway. With increasing treatment duration, the number of DEGs enriched in this biological pathway also increased. Among them, a total of 19 auxin signal transduction-related DEGs were predicted (Fig. 3b: 13 upregulated, six downregulated): five encoded transport inhibitor response proteins (TIR1: one downregulated, four upregulated), seven encoded AUX/IAA proteins (one downregulated, six upregulated), one encoded auxin response factor (ARF: acetic acid-amido synthetase), three encoded acetic acid-3-amido synthetase proteins (GH3: acetic acid-induced protein) and arg7-like protein (ARG7:1 acetic acid-induced protein, two acetic acid-induced proteins), which were downregulated, and one acetic acid-induced protein was upregulated. In the jasmonic acid signaling pathway (Fig. 3a), seven DEGs (Fig. 3b) were predicted, including two that encoded coronatine-insensitive protein (COI: one upregulated, one downregulated) and five that encoded the jasmonate-zim domain protein (JAZ: all upregulated). In the abscisic acid signal transduction pathway (Fig. 3a), a total of eight DEGs were predicted (Fig. 3b), among which three encoded abscisic acid receptor (PYR/PYL: one downregulated, two upregulated) and three encoded phosphatase 2C protein (PP2C: PP2C). Both genes encoded ABSCISIC ACID-INSENSITIVE 5-like protein (ABF: both downregulated); three DEGs were related to gibberellin signal transduction (Fig. 3b), two encoded transcription factors (TFs: upregulated), one encoded gibberellin receptor (GID1: downregulated); three were ethylene signal transduction DEGs (Fig. 3b: these were upregulated, but their expression decreased with increasing treatment duration); and one encoded the ethylene-responsive transcription factor (ERF). ETHYLENE INSENSITIVE 3-like 3 (EIL3) and the DNA-binding domain in plant proteins such as APETALA2 and EREBPs (AP2/ERF); a salicylic acid signal transduction DEG (Fig. 3b) encodes the basic form of pathogenesis-related protein 1-like (PR-1: downregulated). Taken together, these results suggest that lanthanum treatment affects plant hormone signal transduction pathways, thereby altering plant growth and development. In addition, the auxin signal transduction pathway and jasmonic acid signal transduction pathway are significantly enriched in response to lanthanum treatment.

Figure 3. 

Plant hormone signal transduction pathway and category of DEGs in this pathway. (a) Expression heatmap of genes associated with the auxin, jasmonic acid, and abscisic acid signal transduction pathways and related enriched DEGs. (b) Number of DEGs annotated to different plant hormones.

DEGs with significantly different expression levels under lanthanum treatment

To further explore the gene expression patterns in the roots of Chrysanthemum naktongense in response to lanthanum treatment, we analyzed the KEGG pathways enriched with upregulated DEGs under the two treatment durations and the DEGs in each pathway. The results revealed that there were seven KEGG pathways in which both TR-24 vs CK-0 and TR-72 vs CK-0 showed enrichment of upregulated DEGs (Supplementary Table S4). These pathways included stilbenoid, diarylheptanoid, and gingerol biosynthesis (a total of 13 DEGs); longevity-regulating pathways, including multiple species (a total of 26 DEGs); phenylpropanoid biosynthesis (a total of 27 DEGs); flavonoid biosynthesis (a total of 14 DEGs); plant hormone signal transduction (a total of 41 DEGs); plant‒pathogen interactions (a total of 30 DEGs); and sulphur metabolism (a total of 11 DEGs). Among them, the DEGs significantly upregulated in TR-24 vs CK-0 and TR-72 vs CK-0 in each pathway were as follows (a total of 39 DEGs, Supplementary Table S5): stilbenoid, diarylheptanoid, and gingerol biosynthesis, three DEGs; longevity-regulating pathway, 12 DEGs in multiple species; phenylpropanoid biosynthesis, two DEGs in flavonoid biosynthesis; and 10 DEGs in plant hormone signal transduction. There were seven DEGs related to plant‒pathogen interactions and four DEGs related to sulphur metabolism. A total of 21 metal transport-related genes were predicted among the 39 upregulated DEGs mentioned above. These genes included those for heat-activated protein 70 (Hsp70), mitogen-activated protein kinase (MAPK), JAZ, F-box, WRKY, ERF, and GH3. Among them, the largest number of genes were predicted for Hsp70, at 12. These findings suggest that the Hsp70 gene may play a crucial role in the intrinsic mechanism by which Chrysanthemum roots respond to lanthanum treatment.

Identification and bioinformatics analysis of members of the Hsp70 family of Chrysanthemum naktongense based on transcriptome data

To deeper understand the Hsp70s in Chrysanthemum naktongense after Pfam, SMART, and NCBI CDD-search analysis, a total of 16 Hsp70 members were identified from the transcriptome data in this study, and they were named CnHsp70-1 to 16 according to the order of their appearance in the transcriptome data. The physicochemical characteristics of the 16 CnHsp70 family members were comprehensively analyzed, revealing that the protein sequence lengths varied from 265 amino acids in CnHsp70-7 to 864 amino acids in CnHsp70-12 (Table 2). The isoelectric points (pI) spanned a range from 5.11 for CnHsp70-8 to 8.71 for CnHsp70-9, classifying CnHsp70-9, CnHsp70-7, and CnHsp70-15 as basic proteins, while the remaining 13 members, with pI values below 7, were characterized as acidic proteins. Molecular weights ranged from 29.47 kDa for CnHsp70-7 to 96.54 kDa for CnHsp70-12. Instability indices varied from 30.19 for CnHsp70-5 to 50.64 for CnHsp70-9, with CnHsp70-9, CnHsp70-13, CnHsp70-16, and CnHsp70-7 identified as unstable proteins due to their instability indices exceeding 40, while the others were considered stable. The Grand average of hydropathicity (GRAVY) values, which were all negative, ranged from −0.461 for CnHsp70-13 to −0.061 for CnHsp70-6, indicating that all 16 CnHsp70 proteins are hydrophilic in nature. Subcellular localization prediction results showed that nine CnHsp70 members were localized in cytoplasm, four in chloroplast, and one member in extracellular space, mitochondrial, and plasma membrane.

Furthermore, we conducted a phylogenetic analysis of 18 Hsp70 gene family members from Arabidopsis thaliana, 32 from rice (Oryza sativa), and 16 from Chrysanthemum naktongense using an evolutionary tree (Fig. 4a). The resulting 66 Hsp70 genes were categorized into seven distinct clusters, designated as Group I to Group VIII. Specifically, the 16 CnHsp70 members were found distributed across Groups II, IV, V, VI, and VII. Additionally, we utilized the online MEME software to analyze the 16 CnHsp70 protein sequences, and the results (Fig. 4b) revealed the presence of 10 distinct motifs, designated as Motif 1 through Motif 10. Each member contained between 5 to 10 motifs, with Motif 1, Motif 2, and Motif 8 being notably conserved across all CnHsp70 member sequences.

Figure 4. 

Phylogenetic tree of Hsp70 family members of Arabidopsis thaliana (At), Oryza sativa (Os), and Chrysanthemum naktongense (Cn) and structure of 16 Hsp70 motifs of Chrysanthemum naktongense. (a) Phylogenetic tree of Hsp70 family members of AtHsp70-1 to 18, OsHsp70-1 to 32, and CnHsp70-1 to 16. Different colors are used to distinguish different groups. The asterisk indicates CnHsp70-1 to 16. (b) Phylogenetic tree of CnHsp70-1 to 16 and their motifs, diverse colors symbolize distinct motifs.

DEGs determined by qRT-PCR and expression analysis of CnHsp70s

To verify the accuracy of the transcriptome results, in this study, such as those related to the Longevity regulating pathway - multiple species, photosynthesis, and minerals, were selected, and real-time fluorescence quantitative verification was carried out on genes related to heavy metal transport and absorption pathways (Fig. 5). Linear regression analysis revealed that the correlation coefficients were similar, R² = 0.89, indicating that the results of qRT-PCR and RNA-Seq were consistent in terms of gene expression trends. Therefore, the transcriptome data is reliable.

Figure 5. 

Six genes were verified via real-time quantitative PCR.

In addition, to further understand lanthanum expression in CnHsp70s in response to treatment at different times, we analyzed lanthanum expression in 16 Hsp70 members identified above for 24 h and 72 h (Fig. 6). The analysis revealed that the majority of the 16 CnHsp70 family members exhibited upregulation in response to lanthanum treatment at both 24-h and 72-h intervals, with expression patterns that were largely congruent with the transcriptomic data. This concordance validates the authenticity of the transcriptomic data. Furthermore, a subset of 10 CnHsp70 members, including CnHsp70-2, -3, -5, -6, -8, -10, -12, -13, -15, and -16, displayed upregulation at both time points, with their expression levels escalating in tandem with the duration of lanthanum exposure (Supplementary Fig. S1).

Figure 6. 

Expression levels of 10 Hsp70 members of Chrysanthemum naktongense in RNA-seq and qRT-PCR. The blue line represents the results of RNA-seq and the yellow represents qRT-PCR.

REEs are nonessential heavy metals for plants, and long-term plant exposure to REEs negatively affects plant growth and development^[34]. Plant hormones play important roles in the response of plants to abiotic stress. In addition, many studies^[35−37] have suggested that plant hormones, as signal transduction molecules, can reduce the toxicity of heavy metals in plants, thus promoting plant growth and development^[38]. Compared with those in the control group, a total of 41 plant hormone signal transduction pathways related to upregulated DEGs were identified after lanthanum treatment for 24 h and 72 h (Fig. 3).

Auxin plays an important role in plant resistance to abiotic stress. Some studies have shown that the upregulation of auxin-regulatory genes under heavy metal treatment drives auxin biosynthesis and signaling^[39,40]. Previous studies have shown that genes encoding AUX/IAA proteins involved in plant hormone signaling in maize are significantly upregulated under lanthanum treatment^[41]. Our results are essentially consistent with these findings. In this study, we predicted a total of 19 DEGs related to auxin signal transduction (Fig 3a−c), among which 13 DEGs were upregulated, and seven genes encoding AUX/IAA proteins were upregulated, accounting for 85.71% of the above 13 DEGs.

Moreover, with increasing treatment duration, the expression of DEGs related to the auxin signal transduction pathway increased. Six DEGs were downregulated, and their expression levels decreased with increasing treatment duration. These results indicated that low concentrations of lanthanum promoted the expression of DEGs involved in the auxin signal transduction pathway in the roots of Chrysanthemum naktongense. Moreover, the response of the auxin signal transduction pathway to treatment with low concentrations of lanthanum increased with increasing treatment duration. Therefore, lanthanum at a low concentration may promote the auxin signal transduction pathway in the root system of Chrysanthemum naktongense. In addition, jasmonic acid can effectively mitigate the effects of heavy metals on plants and help plants adapt to their environment under adverse conditions^[42]. Zhao et al. reported that endogenous jasmonic acid could effectively limit the absorption and upward transport of Cd by tomato roots^[43]. Jasmonic acid alleviates Cd-induced toxicity in Arabidopsis through signal transduction pathways^[44]. In this study, we predicted a total of seven DEGs related to jasmonic acid signaling pathways (Fig. 3a, b & d), six of which were upregulated, including five genes encoding JAZ proteins and one gene encoding a COI. These results suggest that lanthanum treatment may promote jasmonic acid signal transduction pathways in the roots of Chrysanthemum naktongense. Abscisic acid, jasmonic acid, ethylene, and salicylic acid are generally considered stress hormones and play highly important roles in plant response to Cd stress^[45,46]. Abscisic acid (ABA) is a key plant hormone involved in regulating plant response to various abiotic stresses, such as drought, heat, cold, and heavy metals^[47]. Wang et al. demonstrated that arsenic toxicity upregulated the expression of genes involved in ABA signal transduction pathways in tomato roots^[48]. Other studies have shown that arsenic toxicity upregulates the expression of the ABA signal transduction regulatory genes OsPP2C4, OsPP2C5, OsbZIP10, and OsbZIP12 in rice seedlings^[49]. Similar to the above results, in this study, in the abscisic acid signal transduction pathway (Fig. 3a, b & e), low-concentration lanthanum treatment downregulated one gene encoding the PYR/PYL protein and two genes encoding ABF. These results suggest that lanthanum at low concentrations may inhibit the expression of genes involved in the ABA signal transduction pathway. However, the underlying mechanism by which plant hormone pathways respond to lanthanum treatment needs to be further investigated.

DEGs that regulate the response of the Chrysanthemum naktongense root system to lanthanum treatment

Heavy metal treatment induced differential expression in genes encoding stress-responsive proteins. In this study, we found that 21 genes encoding stress response proteins, including Hsp70, MAPK4, JAZ, F-box, WRKY, ERF, and GH3, were consistently and significantly (p < 0.05) upregulated at both treatment durations. Twelve of the genes encoded the Hsp70 protein. Moreover, the expression heatmap (Fig. 4a) revealed that the expression of the genes encoding these 12 protein members of Hsp70 essentially increased with increasing lanthanum treatment duration. Under abiotic stress such as high temperature, drought, high salt, and heavy metals, Hsp70 can be rapidly synthesized to improve plant stress resistance^[50]. Some studies have shown that Hsp70 is related to the internal mechanism of plant response to heavy metals. Shaheen et al. reported that the expression level of the Hsp70 protein in giant reeds increased with increasing metal concentration and concluded that the Hsp70 protein may significantly improve plant tolerance to metal stress^[51]. The overexpressing transgenic strain of RcHsp70 in Chinese rose increased the resistance of the plant to heavy metal stress^[52]. Genome-wide analysis and expression profile analysis of the SlHsp70 gene family in Nightshade by Abbas et al. revealed that the expression level of SlHsp70-11 in roots was greater under Cd²⁺ stress^[53]. These findings are essentially consistent with our findings. To gain a deeper comprehension of the Hsp70s expression profile in Chrysanthemum naktongense treated with lanthanum, in this study, we identified 16 CnHsp70 members in the transcriptome data by bioinformatics. The expression analysis found that a subset of 10 CnHsp70 members, including CnHsp70-2, -3, -5, -6, -8, -10, -12, -13, -15, and -16, after treatment with 100 mg·L⁻¹ lanthanum for 24 h and 72 h, the expression level was on the rise (Supplementary Fig. S1). These findings indicated that Hsp70s may enhance the ability of plants to withstand heavy metal stress. Furthermore, previous research has demonstrated that plant hormones can influence the expression levels of Hsp70^[18]. In this study, genes involved in plant hormone signal transduction pathways and Hsp70s were found to be up-regulated following lanthanum treatment (Figs 3, 5 & 6). Consequently, both plant hormones and Hsp70s may play a beneficial role in C. naktongense's response to lanthanum exposure. Nevertheless, the potential mechanism remains to be further explored.

In this study, we preliminary revealed the molecular mechanism by which the root system of Chrysanthemum naktongense responds to lanthanum treatment. Based on TR-24 vs CK-0 (3,266) and TR-72 vs CK-0 (5,933) analysis, it was found that, with increasing lanthanum treatment duration, the number of DEGs increased. In addition, we found that DEGs were significantly (p < 0.05) enriched in multiple KEGG pathways under the two treatment durations, among which DEGs in the plant hormone signaling pathway were most enriched. In this pathway, low-concentration lanthanum treatment promoted the signal transduction pathways of auxin and jasmonic acid and inhibited the signal transduction pathways of abscisic acid, thus effectively enhancing the resistance of the root system to lanthanum treatment. In addition, 16 members of CnHsp70s were identified, and a subset of 10 CnHsp70 members, including CnHsp70-2, -3, -5, -6, -8, -10, -12, -13, -15, and -16, after treatment with 100 mg·L⁻¹ lanthanum for 24 and 72 h, the expression level was on the rise. In future studies, researchers should further elucidate the plant hormone-mediated root response to lanthanum treatment, investigate the role of the key genes screened herein, and identify genes that have the potential to improve the accumulation of the rare earth element lanthanum in Chrysanthemum naktongense.