Comparative transcriptome analysis of the resistance mechanism of Hemerocallis citrina Baroni to Puccinia hemerocallidis infection

ABSTRACT

Hemerocallis citrina Baroni (daylily) is an important cash crop in many regions of China. Daylily rust caused by Puccinia hemerocallidis (Ph) is a major disease in H. citrina. Rust infection on the leaves of two resistant and susceptible varieties of H. citrina was examined. Transcriptome sequencing and defense enzyme activity were used to analyze the differential expression of related genes after 0, 72, 144, 192, and 288 h, to reveal the response mechanism of H. citrina to rust infection and screen-related resistance genes. The results suggested that 72 h after Ph infection was the key time point for determining differences between resistant and susceptible varieties. DEGs were significantly enriched in phenylpropanoid biosynthesis, plant–pathogen interaction, plant hormone signal transduction, and the MAPK signaling pathway. TFs WRKY33 and CYP98A are key genes closely related to the disease-resistance response of H. citrina. This study is the first large-scale transcriptomic analysis of H. citrina.

KEYWORDS:

• Daylily
• enzymatic activities
• RNA-seq
• transcription factors
• resistance genes

1. Introduction

Hemerocallis citrina Baroni is a perennial herbaceous plant of Hemerocallis (Asphodelaceae), which is known as the traditional mother's flower in Chinese culture and has a history of more than 1000 years (Qing et al. 2021; Ti et al. 2022). It is an important cash crop in many regions of China, mainly produced in Hunan, Gansu, Shanxi, Shaanxi, Henan, Jiangsu, Sichuan, and other provinces. Hemerocallis citrina has a variety of uses, including in gardening, and has a high ornamental value; its flower buds after steaming, drying, and processing into dried vegetables, namely huang hu cai in Chinese, are sold at home and abroad due to their rich nutritional value. Hemerocallis citrina is also a traditional Chinese herbal medicine that can improve lactation, sleep, depression, stomach pain, diuresis, and swelling (Xu et al. 2016; Li et al. 2017; Liang et al. 2021; Zhong et al. 2021; Liu et al. 2022).

However, daylily rust is an important disease in H. citrina caused by Puccinia hemerocallidis Thuem (Smith 2009; Buck 2013). This disease was originally found in Asia, but it gradually spread to North America in 2000 (Williams-Woodward et al. 2001; Hernández et al. 2002), South America in 2001 (Inokuti et al. 2012), Oceania in 2002 (Buck 2013), Central America and the Caribbean in 2006 (Buck 2013), South Africa in 2008 (Mostert et al. 2008), and Europe in 2016 (Silva et al. 2016).

Daylily rust is a serious hazard, infecting the leaves, flowers, and stems and seriously affecting the yield and quality of H. citrina. The incidence of diseased leaves is more than 50%, often causing yield reductions of 10–50%, thus restricting the further development of the H. citrina industry in China. Improving rust resistance is therefore a prime objective in H. citrina breeding. To avoid yield and quality losses, screening, breeding, and utilizing resistant varieties of H. citrina are the most economical and effective measures to prevent and control daylily rust, which is an important problem to be solved in current production. However, the research on the molecular mechanism of the interaction between H. citrina and P. hemerocallidis has not yet been carried out, which seriously hinders the development of H. citrina disease-resistance breeding projects.

RNA-seq technology has been widely used to study the mechanisms of plant disease resistance. RNA-seq technology can be used to analyze the interactions of plants in response to pathogen infection at a deeper level, discover more genes and metabolic pathways related to disease resistance, and understand the immune response of plants more systematically. Using RNA-seq technology, Dorostkar et al. (2022) explained the molecular mechanism of leaf rust resistance in Aegilops tauschii. Zhao et al. (2022) revealed the candidate genes and related metabolic pathways of powdery mildew resistance in melon (Cucumis melo L.) responding to Podosphaera xanthii infection. Meng et al. (2021) identified genes conferring resistance and susceptibility to Phytophthora nicotianae in tobacco. Cui et al. (2017) identified DEGs and long non-coding (lnc) RNAs between resistant and susceptible tomato, and found that lncRNA16397 conferred resistance to Phytophthora infestans by co-expressing glutaredoxin. Cremer et al. (2013) found that genes involved in the phenylpropanol pathway and terpenoid synthesis were significantly upregulated, while genes involved in photosynthesis were significantly downregulated after lettuce was infected with Botrytis cinerea.

In this study, RNA-seq technology was used to sequence the transcriptome of two resistant and susceptible varieties of H. citrina. This is the first large-scale transcriptomic analysis of H. citrina. The resistance mechanism of H. citrina response to rust infection was studied, and the differentially expressed gene information of H. citrina to daylily rust infection was obtained. Bioinformatics analysis was used to screen candidate genes and metabolic pathways related to rust resistance and provide new insights into the molecular mechanism of the interaction between H. citrina and the rust fungus. These results will bring more effective and long-term measures for the prevention and control of daylily rust and for improving the yield and quality of H. citrina and will provide a scientific theoretical basis for the breeding of varieties resistant to daylily rust.

2. Materials and methods

2.1. Plant culture

Two H. citrina varieties were used for inoculation with Puccinia hemerocallidis (Ph). One was ‘Qiezihua' (QZH) (Ph resistant) and the other one was ‘Malanhua' (MLH) (Ph susceptible), they were obtained from the germplasm resource nursery of the Institute of Suqian, Jiangsu Academy of Agricultural Sciences. Sixty healthy QZH and MLH were selected. The leaves and excess fleshy roots were cut off and transplanted into the nutrient bowl, and water was poured into the plant’s nutrient bowl until the substrate was saturated. The nutrient bowl was then placed in a light incubator at 25°C with a humidity of 65–70% and a light/dark cycle of 14 h/10 h. Two weeks later, when H. citrina had grown new leaves, they were ready for inoculation. Before inoculation, H. citrina needed to be sprayed with a watering can, and then the lid and vent were closed to moisturize the daylily for 4–8 h.

2.2. Pathogen inoculation

An isolate of Ph was collected from an infected H. citrina plant (‘Mengzihua’). Urediniospores were suspended in 0.05% (v/v) Tween 20 (Mueller et al. 2003). The dewaxed leaf was repeatedly rubbed with a sponge to slightly damage the leaf and then sprayed with the spore suspension until runoff, and the hood was covered and the vent closed. Every day, the plants were cultured in light for 14 h at 28°C with a humidity of 90% and in the dark for 10 h at 22°C with a humidity of 85%. The plants were moisturized for 48 h, and the smother cap was opened for ventilation.

2.3. Enzymatic determination

To clarify the changes in defense enzyme activities in resistant and susceptible H. citrina cultivars infected with rust, leaf samples (about 6 g) treated for 0 (control), 72, 144, 192, or 288 h were collected. The samples were wrapped with aluminum foil, quickly frozen in liquid nitrogen, they were stored at −80°C. The activities of SOD (superoxide dismutase), CAT (catalase), POD (peroxidase), and PPO (polyphenol oxidase)-related defense enzymes were determined. All enzyme extraction procedures were conducted at 4°C. SOD activity was measured using nitroblue tetrazolium (NBT) with a SOD kit (Suzhou Keming Biotechnology Co., Ltd., Suzhou, China). CAT, POD, and PPO activities were measured using CAT, POD, PPO kit (Nanjing Cluster Biotechnology Co., Ltd., Nanjing, China).

SOD activity was measured based on the reduction of NBT by the superoxide anion generated by the riboflavin system under illumination. One unit of SOD (U) was defined as the amount of crude enzyme extract required to inhibit the reduction rate of NBT by 50% (Chen et al. 2009). H₂O₂ has a characteristic absorption peak at 240 nm, and CAT can decompose H₂O₂, resulting in a decrease in the reaction solution with increasing reaction time at 240 nm. CAT activity can be calculated according to the rate of absorbance change. The catalytic degradation of 1 μmol of H₂O₂ per g tissue per minute was defined as one unit of CAT (U). POD catalyzes H₂O₂ oxidation of specific substrate and has characteristic light absorption at 470 nm. One unit of POD activity was defined as a change of 0.01 in absorbance at 470 nm. PPO can catalyze catechol to produce quinone, which has characteristic light absorption at 525 nm. One unit of PPO (U) was defined as a change of 0.005 in absorbance at 525 nm. Three repetitions were completed (Mohammadreza and Ali 2016).

2.4. RNA extraction, library preparation, and sequencing

To determine the changes in the gene expression levels of MLH and QZH at the stage of interaction with Ph, about 1 g of leaf samples treated for 0 (control), 72, 144, 192, and 288 h was collected, sealed in plastic bags, and immediately frozen in liquid nitrogen. Three biological replicates were set for each condition. All leaves were stored at −80°C before RNA extraction. Total RNA was extracted with a DP441 RNA extraction kit (Tiangen Biochemical Technology Co. Ltd., Beijing, China) and sent to Genepioneer Biotechnologies (Nanjing, Jiangsu, China) for transcriptome sequencing. RNA was checked for purity (OD260/280), concentration, and nucleic acid absorption peak using Nanodrop, and for integrity using an Agilent 2100 Bioanalyzer. After qualifying the total RNA samples, oligo (dT) magnetic beads were used to enrich mRNA. Fragmentation buffer was added to the obtained mRNA to generate short fragments. Then, using mRNA as a template, the first cDNA strand was synthesized with six base random hexamers, and the second strand of cDNA was synthesized by adding buffer, dNTPs, RNase H, and DNA polymerase I. Next, the cDNA was purified by Ampure XP beads. The purified double-stranded cDNA was end-repaired, A-tailing was added, and sequencing connectors were connected. AMPure XP beads were used for segment size selection. Finally, the cDNA library was enriched via PCR. After library construction, a NovaSeq 6000 platform was used for sequencing (Jiangtao et al. 2021; Dorostkar et al. 2022; Gurjar et al. 2022).

2.5. Bioinformatics analysis of RNA-seq data

The initial results of transcriptome sequencing were in the form of original images, which could be employed for base recognition and transformation via base calling to obtain the original sequence data (raw data or raw reads). After obtaining high-quality sequencing data, Trinity software (V2.11.0) (Grabherr et al. 2011) was used for sequence assembly of the obtained high-quality clean data to generate the transcript sequence of each gene (transcript). The longest transcript in each gene was selected as a unigene. There is no reference genome for H. citrina, so the Unigene library of the species was constructed using the assembled transcripts, based on which bioinformatics analysis, including gene expression analysis and gene functional annotation, was conducted.

2.6. Unigene annotation

Unigenes were annotated according to Blastx comparison results from different databases. BLAST (Altschul et al. 1997) software was used to compare unigene sequences with those of the non-redundant (NR) protein sequence (Deng et al. 2006), Swissprot (Apweiler et al. 2004), Gene Ontology (GO) (Ashburner et al. 2000), Cluster of Orthologous Groups of proteins (COGs) (Tatusov et al. 2000), Clusters of Orthologous Groups for eukaryotic complete genomes (KOG) (Koonin et al. 2004), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2004) databases. After predicting the amino acid sequences of unigenes, HMMER software was used to compare them with the Pfam database (Finn et al. 2014). Unigene annotation information was then obtained.

2.7. Principal component analysis (PCA)

To demonstrate the differences between samples, PCA was performed using the princomp function in R package (V3.6.2), namely dimension reduction, and a PCA diagram was drawn.

2.8. Transcription factor analysis

iTAK (V1.7a) software was used to predict plant transcription factors (TFs). DEG-associated protein sequences in each group were used to predict the TFs, and the TFs regulating the rust resistance mechanism of H. citrina cultivars were predicted by comparing the differences in TFs in rust-resistant and rust-susceptible H. citrina varieties at four time points.

2.9. Identification of differentially expressed genes (DEGs)

DEGs among samples were analyzed by calculating differential multiples using the DESeq2 (1.26.0) program package in R (v3.6.2). The screening criteria for differential genes are |log₂FC|≧1 and p-value≦0.05 (p-value obtained by Benjamini–Hochberg correction methods). log₂FC≧1 indicated an upregulated gene, and log₂FC≦1 indicated a downregulated gene. A Venn diagram was used to analyze the expression of DEGs 72, 144, 192, and 288 h after inoculation with the rust fungus to identify the DEGs at the four time points and to analyze the functions of these genes.

2.10. GO and KEGG enrichment analyses of DEGs

GO function enrichment analysis was conducted for MLH and QZH at important treatment time points after rust inoculation. By classifying DEGs and all expressed genes (background genes) at the secondary GO term, the distribution of DEGs and background genes in each secondary item was obtained. The significant KEGG pathway enrichment of DEGs was used to predict the most important biochemical metabolic pathways and signal transduction pathways of DEGs at key time points after rust inoculation, and key candidate genes of the anti-rust mechanism of daylily were screened from key metabolic pathways. A heatmap of key candidate genes before and after rust infection in MLH and QZH was generated based on FRKM values using TBtools v1.0 software.

2.11. Verification of RNA-seq results by qRT-PCR

To verify the accuracy of the obtained differential gene expression patterns, we screened 19 DEGs for verification analysis by qRT-PCR. qRT-PCR using SYBR Premix Ex TapTMII (TaKaRa, Japan) was conducted with a real-time PCR detection system (Bio-Rad, Hercules, California, USA). The primers used in the reactions are listed in Table S1. The total reaction mixture volume of 20 μL consisted of 0.5 μL of each primer (10 μM), 2 μL of template (10× diluted cDNA), 10 μL of SYBR Premix Ex Tap™II, and 7 μL of ddH₂O along with distilled water as a template for blank control. The qRT-PCR program followed at 95°C for 5 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. The specificity of the reaction was verified using melting curve analysis, and the relative gene expression was quantified using the Pfaffl 2^−ΔΔCt method (Pfaffl 2001; Falcão et al. 2022). All samples were assayed in triplicate.

3. Results

3.1. Phenotypic response of H. citrina to rust infection

The phenotypic response of H. citrina infected with rust in the field is shown in Figure 1(a). After 14 days of rust infection, the susceptible cultivar ‘Malanhua' (MLH) showed high susceptibility (the rust spot area accounted for more than 10% of the leaf area) (Figure 1(b)). Meanwhile, the resistant cultivar ‘Qiezihua' (QZH) showed high resistance (no disease spots on the leaf surface) (Figure 1(c)). Based on the results of phenotypic analysis, a comparative transcriptomic analysis of the mechanism of rust resistance of QZH with high resistance and MLH with high sensitivity was carried out.

Figure 1. Response of Hemerocallis citrina (H. citrina) to rust infection. (a) Phenotypic response of H. citrina infected with rust in the field. (b) Phenotypic response of rust fungus after 14 days of infection with susceptible cultivar ‘Malanhua’ (MLH). (c) Phenotypic response of rust fungus after 14 days of infection with resistant cultivar ‘Qiezihua’ (QZH).

3.2. Defense enzyme activity

Before rust inoculation of H. citrina leaves (0 h), the enzymatic activities of SOD, POD, CAT, and PPO in the susceptible variety MLH were higher than those of the resistant variety QZH (62.68, 231.68, 6.98, and 9.74%, respectively). The enzyme activity obviously changed after rust infection of H. citrina. SOD activity in the leaves of the two varieties decreased after infection with Ph. After 72 h, the SOD activity of QZH decreased to the lowest level and then began to increase in a wavelike pattern. After 288 h, the SOD activity of QZH increased by 22.2% compared with that before inoculation; in contrast, the SOD activity of MLH decreased by 4.73% after 288 h compared with that before inoculation. The fluctuation in SOD activity throughout the process of infection was very obvious in QZH but not in MLH, which may indicate that QZH had a faster SOD response than MLH and that it could more effectively remove reactive oxygen species (ROS) produced by fungus infection in the cells (Figure 2(a)). POD enzyme activities in the leaves of the two varieties were contradictory after being infected with Ph. After 72 h, the POD activity of QZH rapidly increased to the highest level and then decreased, while the POD activity of MLH decreased to the lowest level. However, the POD activity increased in MLH, indicating that resistance in QZH was enhanced by rapidly increasing POD activity after rust infection. However, the POD activity of MLH was not effectively activated but decreased after pathogen infection (Figure 2(b)). The changes in CAT enzyme activities and POD enzyme activities in the leaves of the two varieties were consistent after rust infection. The CAT enzyme activity of QZH increased rapidly after pathogen infection to enhance plant resistance, while the CAT activity in MLH failed to be effectively increased, and its enzyme activity gradually decreased after pathogen infection (Figure 2(c)). The changes in PPO enzyme activity were somewhat complicated. After 72 h of infection with daylily rust, the enzyme activities of both varieties increased rapidly, and MLH reached its highest level and then began to decrease. The PPO enzyme activity in QZH increased after rust infection, decreased after 72 h, began to increase after 144 h, and was higher than that in MLH after 288 h. These results indicate that the PPO enzyme activity of daylily-infected plants was activated rapidly after rust infection, but the PPO enzyme activity of the resistant variety QZH was higher than that of the susceptible variety MLH (Figure 2(d)).

Figure 2. Changes in the enzyme activity of Hemerocallis citrina (H. citrina) before and after Puccinia hemerocallidis (Ph) infection. (a) SOD activity before and after Ph infection. (b) POD activity before and after Ph infection. (c) CAT activity before and after Ph infection. (d) PPO activity before and after Ph infection.

The changes in enzyme activities in daylily-infected plants showed that after 72 h of infection, the enzyme activities of the two varieties were at the peak, either rising to the highest level or falling to the lowest level. Serba et al. (2015) found a massive upregulation of transcripts at 60 hpi in switchgrass, which suggests the widely distributed regulation of transcriptional activity 2 d after rust inoculation. Dhariwal et al. (2011) reported differential expression between resistant and susceptible cultivars after 48 hpi, with the highest abundance appearing at 96 hpi. In this study, 72 hpi was the key time point to determine the difference between resistant and susceptible varieties, and the enzyme activities at this time point reflected the gene expression differences between the two cultivars. The resistance of resistant varieties was enhanced by increasing the activity of defense enzymes, while susceptible varieties were highly susceptible to disease and consumed a large amount of defense enzymes when defense enzymes were activated in a short time to resist pathogen infection; therefore, the activity of defense enzymes of susceptible varieties was reduced compared with that before inoculation.

3.3. Transcriptome characteristics of rust-infected H. citrina

To analyze the transcriptome changes of rust-resistant and -susceptible H. citrina leaves and identify the related rust-resistant genes, 30 samples were sequenced. A total of 438,198 transcripts and 191,439 unigenes were obtained by assembly with high assembly integrity. There were 143,394 transcripts with a length of 300–500 nt (32.72%) and 140,619 transcripts with a length of 500–1000 nt (32.09%). There were 80,566 unigenes between 300 and 500 nt, accounting for 42.08% of the total, and 64,335 unigenes between 500 and 1000 nt, accounting for 33.61% of the total (Table S2). In this project, 79,271 unigenes with annotated information were obtained by unigene-selecting BLAST, with a parameter e-value ≤ 10⁻⁵. Of these, 50,956 unigenes were between 300 and 1000 bp in length, and 28,315 unigenes were longer than 1000 bp (Table S3). PCA was used to explain the grouping and variability of transcriptomes between QZH and MLH cultivars of H. citrina. The intra-group repeatability of the combination of QZH and MLH varieties was better, and the difference between groups was significant, and principal components 1 and 2 explained 11.72% and 7.00% of the variation in gene expression, respectively (Choudhury et al. 2017) (Figure 3).

Figure 3. Principal component analysis showing the genetic variation between rust-infected QZH and MLH H. citrina varieties based on transcriptome data.

3.4. Analysis of the expression of differential transcription factors (TFs)

The expression of TFs in QZH was higher than that in MLH, and 72 h was the key defense period for both cultivars to be infected by pathogens. TFs, such as WRKY, MYB, NAC, bHLH, bZIP, AP2/ERF, C2C2, GARP, GNAT, HB, NE-Y, and B3, were expressed at the four time points of rust infection in the two H. citrina varieties, and expression in resistant varieties was higher than that in susceptible varieties. Therefore, these differential TFs were the key TFs regulating the anti-rust defense mechanism against daylily, and their expression was beneficial in improving the resistance of H. citrina to rust pathogen infection (Figure 4).

Figure 4. TF analysis showing the expression of rust-infected QZH and MLH H. citrina varieties. (a) Expression of TFs after infection of QZH by rust. (b) Expression of TFs after infection of MLH by rust.

3.5. Analysis of differentially expressed genes (DEGs)

DEGs were detected in both cultivars at different time points after rust inoculation. In the four treatment groups, a total of 27,978 DEGs were detected in the resistant variety QZH, and 24,583 DEGs were detected in the susceptible variety MLH (Figure 5). At 72, 144, 192, and 288 h after inoculation of MLH leaves, 3252, 1867, 2911, and 4129 DEGs, respectively, were upregulated, and 2608, 2166, 3692, and 3958 DEGs, respectively, were downregulated. In QZH, 2691, 3365, 2585, and 4226 DEGs were upregulated at 72, 144, 192, and 288 h, respectively, while 2814, 4703, 3331, and 4263 DEGs, respectively, were downregulated (Figure 5(a–c)). During the sampling period of 72–288 h, the number of DEGs in QZH was higher than that in MLH. More QZH genes were downregulated than upregulated (Figure 5(c)). A total of 663 DEGs were detected at the four time points in QZH (Figure 4(a)). These 663 genes were mainly enzymes, proteins, and TFs, such as peroxidase, oxalate – CoA ligase, caffeoyl-CoA acyltransferase, glutamyl-tRNA (Gln) amidotransferase, tryptophan synthase, Tryptamine hydroxycinnamoyltransferase, UDP-glycosyltransferas, cytochrome P450 94A1, chlorophyll a-b binding protein, photosystem I subunit, leucine-rich repeat protein, pentatricopeptide repeat-containing protein, bHLH, ERF, bZIP, and NAC. A total of 606 DEGs were detected at the four time points in MLH (Figure 5(b)). These 606 genes mainly include pathogenesis-related protein, peroxidase, pectate lyase, pectinesterase inhibitor, UDP-glucuronate, L-ascorbate oxidase homolog, MYB, and GATA.

Figure 5. Differentially expressed genes (DEGs) across four time points after Ph infection in resistant (QZH) and susceptible cultivars (MLH). (a) Venn diagram of DEGs in the resistant H. citrina cultivar at each of the four time points after Ph infection. (b) Venn diagram of DEGs in the susceptible H. citrina cultivar at each of the four time points after Ph infection. (c) Distribution of upregulated and downregulated gene expression patterns across the four time points after Ph infection in resistant and susceptible cultivars.

3.6. Enrichment analysis of DEGs

After 72 h of daylily rust infection, significant enrichment analysis of GO function showed that 2816 unigene sequences in the DEGs of MLH and QZH were classified into three GOs. The differences in DEGs in the cellular component were mainly manifested in the membrane, membrane part, extracellular region, extracellular region part, and supermolecular complex. The differences in molecular function were mainly manifested in transcription regulator activity, antioxidant activity, nutrient reservoir activity, and molecular carrier activity. The differences in biological processes were related to stimulus, multi-organism processes, immune system processes, detoxification, and rhythmic processes, among others (Figure 6(a)). As shown in the enrichment pathway diagram at 72 hpi, the upregulated DEGs of MLH were mainly enriched in photosynthesis-antenna proteins (ko00196), photosynthesis (ko00195), circadian rhythm-plant (ko04712), flavonoid biosynthesis (ko00941), stilbenoid, diarylheptanoid, and gingerol biosynthesis (ko00945), and plant–pathogen interaction (ko04626) (Figure 6(b)). The upregulated DEGs of QZH were mainly enriched in the MAPK signaling pathway (ko04016), plant hormone signal transduction (ko04075), and phenylpropanoid biosynthesis (ko00940), photosynthesis-antenna proteins (ko00196), photosynthesis (ko00195), circadian rhythm-plant (ko04712), and plant–pathogen interaction (ko04626) (Figure 6(c)). The downregulated MLH DEGs were enriched in phenylpropanoid biosynthesis (ko00940), flavonoid biosynthesis (ko00941), and stilbenoid, diarylheptanoid, and gingerol biosynthesis (ko00945) (Figure 6(d)). The downregulated QZH DEGs were enriched in phenylpropanoid biosynthesis (ko00940), flavonoid biosynthesis (ko00941), stilbenoid, diarylheptanoid, and gingerol biosynthesis (ko00945), and ribosome (ko03010) (Figure 6(e)).

Figure 6. Enrichment analysis of DEGs 72 h post-inoculation (hpi). (a) Enrichment analysis of GO function at 72 hpi. (b) Enrichment analysis of the upregulated DEGs in MLH at 72 hpi. (c) Enrichment analysis of the upregulated DEGs in QZH at 72 hpi. (d) Enrichment analysis of the downregulated DEGs in MLH at 72 hpi. (e) Enrichment analysis of the downregulated DEGs in QZH at 72 hpi.

3.7. Screening of candidate genes related to rust resistance

The relevant differential genes involved in these metabolic pathways were screened. The criteria for screening genes were as follows: large expression difference between samples (log₂(FC)), and high gene expression (RPKM/FPKM greater than 50 in at least one sample). Nineteen key candidate genes were screened (Table S4). A heatmap of 19 key candidate genes before and after rust infection in MLH and QZH showed that after 72 h of daylily rust infection, the PR1, 4CL, CYP73A, CYP98A, and Rboh genes were significantly upregulated, and the Psb Q, Lcha1, Lchb1, and Lcb2 genes involved in photosynthesis were downregulated. The WRKY33 transcription factor was significantly upregulated at 144 h (Figure 7).

Figure 7. Heatmap of key candidate genes. (a) Heatmap of key candidate genes before and after rust infection in MLH. (b) Heatmap of key candidate genes before and after rust infection in QZH.

3.8. Verification of RNA-seq results by qRT-PCR

To verify the reliability of the RNA-seq analysis, the expression of 12 candidate genes related to rust resistance in H. citrina was analyzed by qRT-PCR (Figure 8). The expression patterns of the genes determined by qRT-PCR were consistent with those identified by RNA-seq, indicating that the results of this study were reliable for transcriptional expression. The results showed that HCT, 4CL, WRKY33, PP2C, CDPK, CYP98A, CYP73A, Psb Q, and Lhcb were key genes that regulated the expression of resistance genes through transcription.

Figure 8. qRT-PCR validation of key candidate genes.

4. Discussion

When infected by fungal pathogens, plants initiate their own defense mechanisms to resist pathogens and produce corresponding physiological and biochemical reactions, including changes in the activities of defense enzymes, such as SOD, CAT, POD, and PPO. SOD, POD, and CAT are the three major enzymes in the ROS scavenging system, which can effectively inhibit the damage caused by ROS free radicals to the plant body. SOD catalyzes the disproportionation reaction of the superoxide radical (O₂⁻) to produce less toxic H₂O₂ and O₂ (Clarke et al. 2002; Jelali et al. 2014; Gao et al. 2021). POD (Saed-Moucheshi et al. 2014) and CAT (Camejo et al. 2016) degrade H₂O₂ to non-toxic H₂O and O₂. These three enzymes coordinate and interact with each other in the cell to maintain a normal level of free radicals in the cell, reduce the toxicity caused by free radicals to plants, improve the resilience of plants, and are closely related to the disease resistance of plants.

SOD, as the first enzyme to remove ROS (Zhang et al. 2008), plays an important regulatory role in plant resistance to pathogen infection. When the pathogen infects the host, it releases internal toxins and destroys the plant cell structure. In the process of defending against toxins, hosts produce a large amount of ROS beyond the normal level, which breaks the dynamic balance between ROS production and elimination in plants; so, the change in SOD levels reflects the stress resistance of plants to a certain extent. Aamir et al. (2019) showed that one of the regulatory manifestations of tomato defense after Trichoderma infection is the accumulation of ROS molecules, which leads to a gradual increase of SOD activity from 0 to 72 h with the extension of time. This indicates that the SOD level in tomato could be increased by Trichoderma infection to enhance its defense against Fusarium wilt. In our study, SOD activity in the rust-resistant variety QZH decreased to the lowest level 72 hpi with rust and then gradually increased. POD and CAT activities increased to the highest level 72 hpi and then decreased.

PPO is the main enzyme of polyphenol oxidase. It participates in the oxidation of phenolic substances in plants to produce quinones and in the synthesis of lignin to kill pathogens and inhibit their propagation and also plays a role in disease resistance. Correlations between plant resistance to fungi and PPO activities have been reported in many plant–pathogen interactions. Ballhorn et al. (2010) found that the susceptibility of high cyanogenic accessions to fungal pathogens was strongly correlated to reduced activity of resistance-associated PPO in leaves. Dorleku and Gbewonyo (2021) showed that increased levels of PPO correlated with increased resistance to Phytolyma lata-induced galling in Milicia spp. The PPO activity of resistant varieties increased rapidly, which was higher than that of susceptible varieties. PPO plays an important role in the process of plant disease resistance. The resistance mechanism of PPO is mainly to participate in the secondary metabolism of plants, produce some toxic substances, and directly kill pathogenic fungi.

Plants activate the expression of certain disease-resistance genes through the regulation of related TFs to resist the infection by fungal pathogens. TFs are a class of proteins that regulate gene transcription in living organisms. They can directly bind or indirectly act on specific nucleotide sequences to activate or inhibit the expression of target genes. WRKY TFs play an important role in tuning the cellular responses triggered by pathogens in many plant species (Dong et al. 2003; Chen et al. 2010). Tao et al. (2022) revealed that WRKY33 directly regulates IGS biosynthesis, conferring resistance to Alternaria brassicicola. Some studies have found that Botrytis cinerea infection or FLG22 treatment can activate the transcriptional activity of WRKY33 and induce a plant disease resistance response (Verma et al. 2021). Zhibing Lai’s study showed that the Arabidopsis sp. WRKY33 transcription factor positively regulates plant resistance to necrotrophic fungal pathogens (2010), which is consistent with the results of this study. WRKY33 plays an important regulatory role in the resistance of H. citrina to rust infection. The expression of the WRKY33 gene decreased rapidly after 72 h of rust fungus infection and then increased rapidly after 144 h. The expression of the WRKY33 gene in the disease-resistant variety QZH was significantly higher than that in the susceptible variety MLH, indicating that the WRKY33 gene was positively correlated with rust resistance in H. citrina and that the high expression of the WRKY33 gene effectively inhibited rust fungus infection.

Phenylpropanoid biosynthesis (ko00940) is a key pathway for plant disease resistance (Lai 2010). The metabolic pathway of phenylpropane in plants is to synthesize lignin, flavonoids, and many other phenolic compounds, which are implicated in disease resistance in plants (Karamat et al. 2012; König et al. 2014; Zernova et al. 2014; Lee et al. 2019; Dong and Lin 2021). Phenylpropanoid biosynthesis utilizes very important enzymes, such as HCT, 4CL, CYP73A, and CYP98A, which are closely related to resistance to rust. CYP98A is in the cytochrome P450 family. Fazeelat Karamat et al. (2012) found that CYP98A22 overexpression positively impacted the concentration of furanocoumarins in Ruta graveolens. This is similar to our results, which showed that the CYP98A gene (TRINITY_DN2149_c0_g1) may be the key gene involved in rust resistance. CYP98A was highly expressed in the disease-resistant variety QZH but hardly expressed in the susceptible variety MLH, indicating that this gene is closely related to resistance to daylily rust. CYP98A is mainly involved in phenylpropanoid biosynthesis, styrene acrylic biosynthesis, stilbenoid, diarylheptanoid, and gingerol biosynthesis, and the flavonoid metabolic pathway, which is mainly involved in lignin synthesis. In this experiment, the CYP98A gene was expressed in disease-resistant varieties but not in the susceptible varieties, indicating that this gene improved the rust resistance of H. citrina by controlling lignin synthesis. If expressed in plants, the varieties were rust resistant; if not expressed, the varieties were not rust resistant. This gene may be a key gene for rust resistance.

Finally, according to previous studies and the results of this experiment, we predicted the possible resistance mechanism of H. citrina in response to rust fungus infection, as shown in Figure 9. When the pathogen of rust infects the leaves of H. citrina, its defense mechanism is mainly activated through the synergistic effect of several metabolic pathways to enhance resistance to the rust fungus. The first mechanism is the plant–pathogen interaction pathway. On the one hand, the high-sensitivity response triggered by pathogen-associated molecular pattern (PAMP) is formed after rust infection, and the expression of the Rboh gene generates a large amount of ROS to enhance the disease resistance of tissues. On the other hand, TF WRKY33 is activated to induce the production of rust resistance genes. The second mechanism is plant hormone signal transduction; the metabolites produced by phenylalanine metabolism regulate the expression of the PR-1 gene mediated by the first signal molecule, salicylic acid. The third mechanism is phenylpropanoid biosynthesis, in which the formation of lignin is catalyzed by the HCT enzyme, strengthening the cell wall, and flavonoids in CHS are under the catalysis of enzymes, antioxidants, and phytoalexins. Fourth, the MAPK signal pathway for rust infection can quickly activate WRKY33 transcription activity, inducing a plant disease resistance response. Camalexin synthesis can activate the resistance protein SUMM2 by disrupting MPK4 through pathogen effectors, mediating plant immunity and cell death. In addition, after the rust fungus infects H. citrina, plants affect photosynthesis by regulating the photosynthetic antenna protein, which leads to a change in the circadian rhythm.

Figure 9. Possible resistance mechanism of H. citrina in response to daylily rust fungus infection.

5. Conclusions

In this study, we reported the results of the first large-scale transcriptomic analysis of H. citrina in response to rust infection. This study screened the DEGs related to anti-rust biological pathways, and these gene expression differences were related to protein, protein kinases, TFs, photosynthesis, plant and pathogen interaction path, phenylpropanoid biosynthesis, and flavonoid synthesis and thus were closely related to rust resistance in H. citrina by being involved in the disease resistance reaction. The rust resistance genes of H. citrina were further mined, and the key DEGs were verified by RT–PCR. TF WRKY33 activated the disease resistance response of plants, thus participating in the regulation of rust resistance genes in H. citrina. CYP98A is a key gene related to rust resistance, and it was overexpressed in QZH, a disease-resistant variety. CYP98A is a key enzyme that regulates phenylpropanoid synthesis and flavonoid synthesis. The expression of this enzyme is key for disease-resistant varieties to resist pathogen infection, and it promotes lignin synthesis and increases lignin content in plant cell walls, thereby enhancing rust resistance in H. citrina. Our data has provided some new information on understanding the functionary of CYP98A gene and TF WRKY33, mealwhile, it perhaps facilitates rust resistance genetic improvement in the furture. Furthermore, We believe that our findings could be bring more effective and long-term measures for the prevention and control of daylily rust and for improving the yield and quality of H. citrina and will provide a scientific theoretical basis for the breeding of varieties resistant to daylily rust.

Acknowledgements

We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript. Conceptualization, methodology, validation, formal analysis, writing – original draft preparation, L.J.Z.; investigation, J.M.; resources, S.W.; data curation, S.L. and N.Y.; visualization, F.T.; writing – review and editing, supervision, project administration, funding acquisition, L.L.Z and X.Y. All authors have read and agreed to the published version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The raw data were deposited in the National Centre for Biotechnology Information(NCBI) Sequence Read Archive (SRA) database with accession number PRJNA891118.

Additional information

Funding

This research was funded by the Suqian Science and Technology Program [K201921], the Suqian Agricultural Science and Technology Independent Innovation Project [SQCX202209].

Notes on contributors

Lijie Zhang

Lijie Zhang holds a degree in Biological Sciences from Anhui Institute of Science and Technology, China (2006) and a master's in Microbiology from Fujian Agriculture and Forestry University, China (2009). She is currently an associate professor at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. She has experience in the area of microbiology, molecular biology, plant pathology and plant breeding, working mainly on the following topics: plant and microbial interactions, transcriptomics and Breeding of disease-resistant daylily varieties.

Lingling Zhou

Lingling Zhou holds a degree in Agronomy from Shenyang agricultural university, China (2007) and a master's in Agronomy from Nanjing Agricultural University, China (2010). She is currently an associate professor at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. She has experience in the area of microbiology, molecular biology and plant breeding, working mainly on the following topics: breeding and cultivation techniques of daylily varieties.

Jiali Meng

Jiali Meng holds a degree in Agronomy from Nanjing Agricultural University, China (2009) and a master's in Agronomy from Nanjing Agricultural University, China (2012). She is currently a research associate at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. She has experience in the area of molecular biology, proteomics and plant breeding, working mainly on the following topics: cold tolerance of vegetable, plant grafting technique, proteomic analysis and breeding of low temperature tolerant varieties.

Shaojun Wu

Shaojun Wu holds a degree in Agronomy from Yunnan Agricultural University, China (2008) and a master's in Agronomy from Nanjing Agricultural University, China (2011). He is currently a research associate at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. He has experience in the area of vegetable cultivation, plant physiology and vegetable breeding, working mainly on the following topics: vegetable breeding technology, vegetable continuous cropping obstacle prevention technology, vegetable and disease interaction.

Shuhua Liu

Shuhua Liu holds a degree in Agronomy from Anhui Institute of Science and Technology, China (2017) and a master's in Agronomy from Nanjing Agricultural University, China (2020). He is currently a research assistant at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. He has experience in the area of plant nutrition and plant physiology, working mainly on the following topics: crop growing, plant and disease interaction.

Nianfu Yang

Nianfu Yang holds a degree in Agronomy from Anhui Agricultural University, China (2017), and a master's in Agronomy from Nanjing Agricultural University, China (2020). He is currently a research assistant at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. He has experience in the area of Vegetable science, vegetable breeding for disease resistance, vegetable cultivation, working mainly on the following topics: breeding of vegetable resistance varieties, high efficiency and high quality cultivation technology.

Fufa Tian

Fufa Tian holds a degree in Vegetable Science from Qingdao Agricultural University, China (2000) and a master's in Vegetable Science from Huazhong Agricultural University, China (2003). He is currently an associate professor at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. He has experience in the area of Vegetable culture, vegetable breeding and plant pathology, working mainly on the following topics: vegetable breeding, Vegetable cultivation and disease control technology.

Xiang Yu

Xiang Yu holds a degree in Agronomy from Nanjing Agricultural University, China (2005) and a master's in Agronomy from Nanjing Agricultural University, China (2008). He is currently an associate professor at Institute of Suqian, Jiangsu Academy of Agricultural Sciences. He has experience in the area of vegetable culture, vegetable breeding science and molecular biology, plant pathology and plant breeding, working mainly on the following topics: breeding of high yield and high quality vegetable varieties, efficient fertilization technology and improvement of vegetable quality.