Transcriptome Sequencing Analysis Provides Insights Into the Response to Fusarium oxysporum in Lilium pumilum

Introduction

Lily is a perennial bulbous flower in the monocotyledon subclass, Liliaceae family and Lilium genus, with an extremely high ornamental value. China is a large producer of lily cut flowers. As the planting of lily continues and the planting area of lily expands, lily diseases are becoming more serious. Lily basal rot (also known as root rot or stem rot) has become one of the most important diseases endangering lily. It is a soil-borne disease that mainly harms the underground part of the plant, resulting in wilting, yellowing, and rotting of the bulbs. This disease seriously affecting the yield and quality of lily and resulting in huge economic losses. ¹ Fusarium oxysporum f. sp. lilii is the main pathogen that causes lily basal rot. ² F. oxysporum can survive for a long time as mycelium in the bulbs or as mycelium, chlamydospore, or sclerotium with diseased residue in the soil. At present, chemical control methods are the main measures to prevent and control lily basal rot. However, due to the soil-borne nature and the development of pathogenic resistance of this disease, chemical control methods have shown dwindling effects. ³ At the same time, long-term, large-scale, and repeated use of chemical pesticides causes environmental pollution. Therefore, selective breeding of resistant varieties is the preferred choice to control lily basal rot.

The germplasm resources of Lilium species highly resistant to F. oxysporum have an important role in cultivating lily varieties against this disease. To date, no wild Lilium species or cultivars showing complete resistance to lily basal rot have been found, but several wild Lilium species and cultivars are highly resistant to it. ⁴ Among lily cultivars, Asian lily hybrid is highly resistant and Oriental lily is the least resistant. ⁴ The lily germplasm resources in China are abundant. Some wild species with strong disease resistance, such as Lilium henryi, Lilium pumilum, and Lilium regale, have been widely used in disease resistance breeding. ⁴

Although the germplasm resources of Lilium species resistant to F. oxysporum have been used in disease-resistant breeding, the molecular mechanism of disease resistance is rarely reported. In 2013, a cDNA library of the root of F. oxysporum–infested Lilium regale was constructed by suppression subtractive hybridization and 180 sequences homologous to known proteins in the model plant by sequencing analysis, most of which belonged to pathogenesis-related (PR) 3, PR10, glutathione S-transferase, cytochrome P450, antioxidant enzymes, and peroxidases, were obtained. ⁵ In 2012, Lilium leichtlinii var. maximowiczii Baker was used to construct suppression subtractive hybridization library of lily after induced by F. oxysporum. They identified 6 types of disease-resistant expressed sequence tags (ESTs), including serine/threonine protein kinase, glutathione S-transferase, peroxidase, and cyclophilin homologs. ⁶ The above 2 studies identified some genes of lily in response to F. oxysporum infection, but the mechanism of this resistance has not been investigated in depth. Some studies have used the molecular marker system to construct the genetic map of Asian lily backcross populations. To date, a total of 10 potential quantitative trait loci (QTLs) for F. oxysporum resistance have been mapped. ⁷ , ⁸ However, none of these studies identified the specific genes involved in F. oxysporum resistance. Therefore, the molecular mechanism of resistance to F. oxysporum in lily needs to be further studied.

In this study, L. pumilum that was highly resistant to F. oxysporum was used as the experimental material. High-throughput transcriptome sequencing technology combined with bioinformatics method was applied to identify genes related to F. oxysporum resistance in L. pumilum, understand the function of these genes, and fully integrate the regulatory network, thereby clarifying the molecular mechanism of lily resistance-associated genes in response to pathogen invasion. Our findings will lay a theoretical foundation for the cloning of disease-resistance-associated genes and for disease-resistance breeding.

Materials and Methods

The tissue culture seedlings of wild L. pumilum with the same genetic background were used in the experiment. The pathogen F. oxysporum was isolated from the plants with the symptom of lily basal rot in the Greenhouse of Beijing University of Agriculture. L. pumilum and F. oxysporum were stored in the Lily Breeding Laboratory, College of Landscape Architecture, Beijing University of Agriculture.

F. oxysporum was inoculated on potato dextrose agar (PDA) medium and cultured at 27°C for 7 days. The spores were washed with sterile water and adjusted to 1 × 10⁶ spores/mL for the subsequent inoculation. The tissue culture seedlings of L. pumilum grew at 25°C under a 16/8 h light/dark period. The tissue culture seedlings of L. pumilum with a bulb diameter 1 to 1.5 cm and strong root tissue were selected for F. oxysporum inoculation. For each treatment, 3 bottles with 3 plant tissue culture seedlings in each bottle were used for inoculation, and each seedling was inoculated with 600 µL of spore suspension. L. pumilum tissue culture seedlings were inoculated with sterile water at the same time as the control group. The roots of the tissue culture seedlings in the treatment group and the control group were collected at 6, 12, or 24 h after inoculation and immediately frozen in liquid nitrogen. The samples were stored at −80°C for later use. The 3 time-point samples for F. oxysporum inoculation were named F6h, F12h, and F24h, respectively, and the corresponding samples for the control group were named M6h, M12h, and M24h.

The total RNA of the samples was extracted using the RNAprep Pure Plant Kit (TIANGEN Ltd; Beijing, China). After the obtained sample was tested, the eukaryotic mRNA was enriched using magnetic beads with Oligo (dT). Subsequently, fragmentation buffer was added to break the mRNA into short fragments. The mRNA was used as the template to synthesize a single strand of cDNA, followed by the addition of buffer, dNTPs, DNA polymerase I, and RNase H to synthesize the double-stranded cDNA, which was purified with AMPure XP beads. The purified double-stranded cDNA was first subjected to end repair by adding poly-A tails and linking to the sequencing adaptor, and then AMPure XP beads were used for the selection based on fragment size. Finally, polymerase chain reaction (PCR) amplification was performed, and the PCR product was purified with AMPure XP beads to obtain the final library. After the library was constructed, Qubit 2.0 was used for the preliminary quantification, and the library was diluted to 1.5 ng/µL. The insert size of the library was detected using an Agilent 2100 Bioanalyzer. When the expected insert size was detected, the effective concentration of the library was accurately quantified (effective concentration of the library >2 nM) to ensure the quality of the library.

The constructed cDNA library was sequenced by Illumina HiSeq 2500 sequencing platform (Illumina, Inc, San Diego, CA, USA) from Novogene. The raw reads from sequencing were filtered. Clean reads were collected after removing reads with the adaptor and those of low quality. The clean reads were spliced by Trinity. ⁹ The resulting transcript sequence was used as a reference for the subsequent analysis. The longest transcript in each gene was used as a Unigene. Functional annotation was carried out for the obtained Unigenes with 7 major databases: National Center for Biotechnology Information (NCBI), non-redundant protein (Nr), NCBI nucleotide sequences (Nt), Protein family (Pfam), euKaryotic Ortholog Groups (KOG), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO).

Using Bowtie2 software, the clean reads were aligned to the assembled Unigene library, and the results were statistically analyzed using RNA-Seq by Expectation Maximization (RSEM) software to further obtain the number of read counts for each gene sample corresponding to each gene. ¹⁰ , ¹¹ The expression of the Unigene was calculated using the fragments per kb per million fragments (FPKM) method. ¹² The read-count data were first standardized using trimmed mean of M values (TMM), followed by the differential analysis using DEGSeq, with the screening threshold q < 0.005 and |log2FoldChange| > 1. ¹³ The obtained differentially expressed genes were subjected to GO functional enrichment analysis and KEGG pathway analysis.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis was used to analyze the expression of genes responsive to F. oxysporum identified by RNA-seq. Total RNA of roots of the tissue culture seedlings were extracted according to the above methods. Total RNA was used for synthesizing reverse transcripts using the PrimeScript RT reagent Kit (TaKaRa, Dalian, China) according to the manufacturer’s instructions. cDNA was amplified using the following primer pairs: UTP-glucose-1-phosphate uridylyltransferase, 5P (5′-AAAGGTTTCGGGTGGTTGTG-3′) and 3P (5′-TCCGACGACGGATACAACTG-3′); alpha-1,4 glucan phosphorylase, 5P (5′-CTCCGGCGTAGAAGGTATCG-3′) and 3P (5′-CGAGGATCCAAGGCAAAGAC-3′); phosphoglucomutase, 5P (5′-CGTCAACGATCCCTTCATCA-3′) and 3P (5′-AGGGTTCCTGCATCCAAAGA-3′); NADP-dependent glyceraldehyde-3-phosphate dehydrogenase (GAPDH), 5P (5′-CCCCTCCTCATCGACATAGC-3′) and 3P (5′-GCTCCTTTGTTCCCCAACAC-3′); and 18 rRNA, 5P (5′-GCCTGAGAAACGGCTACCAC-3′) and 3P (5′-ACCAGACTTGCCCTCCAATG-3′). Quantitative reverse transcription polymerase chain reaction was performed using SYBR Premix Ex Taq II (TaKaRa) under the following conditions: 40 cycles at 95°C for 15 s, 60°C for 30 s, and 72°C for 30 s. The qRT-PCR reactions were performed on the BIO-RAD iQ5 (Applied Biosystems, Foster City, CA, USA). Each sample was processed in triplicate, and the relative expression was calculated using 2^−ΔΔCT. ¹⁴ The 18S rRNA was used as references to normalize the expression level of genes.

Results and Discussion

High-throughput transcriptome sequencing technology is an effective means of studying the molecular mechanisms of the plant without whole-genome information. To date, transcriptome sequencing has been used in studies of genetic divergence, ¹⁵ flavonoid biosynthesis, ¹⁶ carbohydrate metabolism, ¹⁷ cold response, ¹⁸ and vernalization ¹⁹ in lily. In this study, transcriptome sequencing was used to identify the genes of L. pumilum that were involved in responses to F. oxysporum infection. Transcriptomes of the root tissues of the L. pumilum at 6, 12, and 24 h after inoculation with F. oxysporum or sterile water (control) were sequenced. A total of more than 50 GB of raw reads were obtained in this study. After removing adaptors, low-quality reads, and contaminants, over 7 GB of clean reads were obtained from each sample (Supplementary Table 1). The raw data (Accession Number: SRA633315) used in the study can be obtained from the sequence read archive (SRA) database. A total of 217 098 Unigenes with the average length of 517 bp and the N50 of 664 bp were obtained after the de novo assembly of the obtained clean reads (Supplementary Table 2). The obtained Unigenes were aligned in 7 large databases, and 83 289 of the Unigenes were annotated in more than 1 database, accounting for 38.36% of the total Unigenes (Supplementary Table 3). The number of Unigenes with successful annotation was the highest in Nr, accounting for 23.06% of the total Unigenes, and lowest in KOG in the KOG database, only 7.89% of Unigenes annotated.

The analysis of differential gene expression at 6, 12, or 24 h after inoculation with F. oxysporum and sterile water was performed. There were 111 differentially expressed genes in F6h vs M6h, of which 61 were upregulated and 50 were downregulated. There were 254 differentially expressed genes in F12h vs M12h, of which 188 were upregulated and 66 were downregulated. There were 2500 differentially expressed genes in F24h vs M24h, of which 2009 were upregulated and 491 were downregulated (Figure 1).

OPEN IN VIEWER

Figure 1.

Volcano plot indicated upregulated and downregulated Unigenes after inoculation using F. oxysporum. L. pumilum was inoculated with F. oxysporum and sterile water, and cultured for (A) 6 h, (B) 12 h, and (C) 24 h.

The results of GO enrichment analysis of the differentially expressed genes showed that the differentially expressed genes of F6h vs M6h were enriched in the category of molecular function (MF) with only 1 subset annotated as catalytic activity. The differentially expressed genes of F12h vs M12h were enriched in the categories of MF and BP (biological process), including 19 subsets with distribution mainly in catalytic activity, single-organism metabolic process, oxidation-reduction process, oxidoreductase activity, carbohydrate metabolic process, and riboflavin metabolic process. The differentially expressed genes of F24h vs M24h were enriched in the categories of BP, CC (cellular component), and MF, including 196 subsets, which were mostly annotated in the BP categories of metabolic process, single-organism metabolic process, and biosynthetic process; the CC categories of macromolecular complex, cytoplasm, and non-membrane-bound organelle; and the MF categories of catalytic activity, oxidoreductase activity, and structural molecule activity (Figure 2).

OPEN IN VIEWER

Figure 2.

Gene Ontology classification of upregulated and downregulated Unigenes after inoculation using F. oxysporum. L. pumilum was inoculated with F. oxysporum and sterile water, and cultured for (A) 6 h, (B) 12 h, and (C) 24 h. BP indicates biological process; CC, cellular component; MF, molecular function.

The results of KEGG enrichment analysis of the differentially expressed genes showed that there were 7 KEGG pathways enriched among the differentially expressed genes of F6h vs M6h. The upregulated genes were significantly enriched in the metabolic pathway of flavonoid biosynthesis (Supplementary Table 4). There were 18 KEGG pathways enriched among the differentially expressed genes of F12h vs M12h. The upregulated genes were significantly enriched in 13 metabolic pathways, mainly including vitamin B6 (VB6) metabolism, flavonoid biosynthesis, phenylalanine metabolism, phenylpropanoid biosynthesis, plant hormone signal transduction, and tyrosine metabolism (Supplementary Table 4). There were 19 KEGG pathways enriched among the differentially expressed genes of F24h vs M24h. The upregulated genes were significantly enriched in 3 metabolic pathways, including biosynthesis of amino acids, proteasome, and ribosome (Supplementary Table 4).

Previous studies have shown that starch and sucrose metabolism changes when pathogens infect plants. ²⁰ , ²¹ Gómez-Ariza reported that defense-related genes were upregulated systemically in rice leaves after adding sucrose to the roots. ²² Sucrose can promote the lupine defense responses against Fusarium. ²³ In this study, starch and sucrose metabolism was the pathway with the most upregulated genes at 12 h after inoculation, and a total of 9 genes were upregulated, including sucrose-phosphate synthase (c117443_g1); UTP-glucose-1-phosphate uridylyltransferase (c95325_g1); aglucan phosphorylase (c112961_g2); alpha-1,4 glucan phosphorylase (c121101_g1); 1,4-alpha-glucan branching enzyme (c108888_g1); phosphoglucomutase (c116359_g1); glucose-1-phosphate adenylyltransferase (c108264_g1 and c117572_g1); and fructokinase-2-like (c94420_g1) (Figure 3). The upregulation of these genes is conducive to transforming starch and glycogen into sucrose. At 24 h after inoculation, 9 genes in the starch and sucrose metabolism pathway were downregulated, including sucrose-phosphate synthase (c113968_g1), pectinesterase (c74217_g1), glucose-1-phosphate adenylyltransferase (c108264_g1, c117572_g1, and c115394_g1), granule-bound starch synthase 2 (c114915_g1), phosphoglucomutase (c116359_g1), alpha-amylase (c117598_g1), and beta-fructofuranosidase (c119079_g2) (Figure 3).

OPEN IN VIEWER

Figure 3.

A heatmap showing the expression profiles of resistance-associated genes. This graph presents the log base 2 of the change fold of gene expressions at 6, 12, and 24 h after inoculation using F. oxysporum. Each row corresponds to 1 gene, and each column corresponds to the 1 time point.

Glyceraldehyde-3-phosphate dehydrogenase, the key enzyme in the glycolysis/gluconeogenesis pathway, is a multifunctional enzyme that is involved in the regulation of reactive oxygen species (ROS), autophagy, and plant immune responses. ²⁴ The knockout of GAPDH in Arabidopsis enhances the resistance to the incompatible pathogen Pseudomonas syringae pv. Tomato. ²⁴ In the present study, 5 genes in the glycolysis/gluconeogenesis pathway were upregulated at 12 h after inoculation with F. oxysporum, including GAPDH (c103781_g1), NADP-dependent GAPDH (c105111_g1), phosphoglucomutase (c116359_g1), glucose-6-phosphate 1-epimerase (c111964_g1), and phosphoglycerate kinase (c109005_g2). The expression levels of these genes were increased by 2.5-fold and 4.9-fold. The expression level of GAPDH (c103781_g1) was 2.5 times the levels of the control, but there was no significant difference between the expression of these genes and the control at 24 h after inoculation (Figure 3).

In the present study, 6 genes in the phenylpropanoid biosynthesis pathway were highly expressed at 12 h after inoculation with F. oxysporum, including 4-coumarate-CoA ligase 2 (c115266_g1), caffeoyl-CoA O-methyltransferase (c105247_g1), cinnamyl-alcohol dehydrogenase (c119082_g1), and peroxidase (c118468_g3, c115650_g1, and c120335_g1). At 24 h after inoculation, the expression levels of 13 genes in this pathway were downregulated, including phenylalanine ammonia-lyase (c103771_g1), cinnamyl-alcohol dehydrogenase (c117898_g1, c118292_g2, and c119633_g1), peroxidase (c115552_g2, c101935_g1, c116520_g1, c101806_g1, c96040_g1, c116501_g1, c120735_g1, and c103286_g1), and cationic peroxidase 1 (c80387_g1) (Figure 3).

Flavonoids are involved in the resistance of wheat against F. graminearum ²⁵ and the resistance of Brassica napus against Xanthomonas campestris pv. Campestris. ²⁶ In the present study, the expression levels of 3 chalcone synthase genes (c113690_g1, c113916_g1, and c113916_g2) in the flavonoid biosynthesis pathway were upregulated at 6 and 12 h after inoculation with F. oxysporum. Caffeoyl-CoA O-methyltransferase (c105247_g1 and c105247_g2), which is involved in the synthesis of feruloyl-CoA, was upregulated at 12 h after inoculation (Figure 3). Caffeoyl-CoA O-methyltransferase is also involved in the phenylpropanoid pathway and lignin production, playing an important role in the resistance process of maize against Cochliobolus heterostrophus, Cercospora zeae-maydis, and Setosphaeria turcica. ²⁷

At 12 h after inoculation with F. oxysporum, some genes in the jasmonic acid (JA), abscisic acid (ABA), and auxin signaling pathways were induced. Jasmonic acid is an endogenous growth regulator in higher plants and a stress signal molecule that will accumulate rapidly when the plant tissues are attacked by pathogens or insects. The transcription factor MYC4 (c114016_g1 and c113516_g3) of the JA signaling pathway was upregulated in the present study. Abscisic acid can negatively regulate the disease resistance in plants by inhibiting salicylic acid (SA) signaling.^28–30 In the present study, the expression of Protein phosphatase 2C 51 (c115633_g1), which negatively regulates ABA signaling, was significantly upregulated. ³¹ GH3 in the auxin signaling pathway can positively regulate SA-mediated plant disease resistance. ³² , ³³ The gene of auxin-responsive GH3 gene family (c117107_g1) was significantly upregulated in this study (Figure 3).

Vitamin B6 participates in the resistance process of Arabidopsis thaliana against Pseudomonas syringae pv. Tomato DC3000 and Botrytis cinerea. ³⁴ Pyridoxal phosphate is the enzyme in the last step of the process of VB6 biosynthesis. In the present study, the expression levels of pyridoxal phosphate (c121490_g2, c121490_g1, and c108372_g2) at 12 h after inoculation with F. oxysporum were 14.9, 13, and 119.4 times those of the control, respectively (Figure 3).

Inoculation of pathogens can result in the accumulation of different amino acids in the host plant. ³⁵ The inoculation of pathogens in Arabidopsis can result in changes in the expression of some enzymes in the biosynthesis of amino acids. ³⁶ , ³⁷ Homoserine is a precursor of threonine, isoleucine, and methionine. The mutation of downy mildew resistant (DMR1)–encoded homoserine kinase leads to the accumulation of homoserine. Mutations in dihydrodipicolinate synthase 2 (DHDPS2) and aspartate kinase 2 (AK2) genes can lead to the accumulation of the amino acids threonine, methionine, and isoleucine. Mutations in the DMR1, DHDPS2, and AK2 genes are reported to enhance resistance to Hyaloperonospora in Arabidopsis. ³⁸ , ³⁹ In this study, a total of 69 synthase genes were upregulated in acid biosynthesis pathways at 24 h after inoculation with F. oxysporum. These enzymes are involved in the synthesis of multiple amino acids, including threonine, isoleucine, tryptophan, tyrosine, lysine, valine, leucine, phenylalanine, arginine, and proline (Figure 4).

OPEN IN VIEWER

Figure 4.

Upregulated genes in amino acids biosynthetic pathway at 24 h after inoculation using F. oxysporum.

The ubiquitin-26S proteasome system is an important post-translational protein regulation system in eukaryotes. It plays an important role in pathogen defense, ⁴⁰ , ⁴¹ abiotic stress tolerance, ⁴² hormone signaling, ⁴³ morphogenesis, ⁴⁴ and chromatin modification ⁴⁵ in plants. In the present study, the expression levels of 24 subunits of the 26S proteasome in lily root were significantly increased to more than 32 times the levels of the control, and the expression levels of all 14 subunits of the 20S proteasome were significantly increased at 24 h after inoculation with F. oxysporum (Figure 5).

OPEN IN VIEWER

Figure 5.

Upregulated genes in proteasome biosynthetic pathway at 24 h after inoculation using F. oxysporum.

Ribosomal proteins can maintain the stability of the ribosomal complex and play an important role in the synthesis of proteins. Several ribosomal proteins are involved in the disease resistance of plants. ⁴⁶ , ⁴⁷ In this study, 86 ribosome protein-coded genes were significantly upregulated 24 h after inoculation, including 37 genes coding for small-subunit proteins and 39 genes coding for large-subunit proteins. All 86 ribosome protein genes were upregulated by more than 32-fold, and most of the genes were upregulated by approximately 1000-fold (Figure 6).

OPEN IN VIEWER

Figure 6.

Upregulated genes in ribosome protein biosynthetic pathway at 24 h after inoculation using F. oxysporum.

Transcription factors can regulate target genes by binding to specific cis-regulatory elements in their promoters. ⁴⁸ The transcription factors of the WRKY, AP2/ERF, MYC, bZIP, and MYB families are involved in the regulation of defense gene expression when attacked by pathogens. ⁴⁹ , ⁵⁰ MYC transcription factors are important regulators of JA-responsive genes, negatively regulating the expression of the defense genes in Arabidopsis. ⁴⁹ , ⁵¹ ERF is a transcription factor of the AP2/ERF subfamily and is involved in the defense responses of plants to pathogenic bacteria via the SA and Ethylene (ET)/JA-dependent signal transduction pathways. ⁵² , ⁵³ The transcription factors of the WRKY family are involved in the protection response of plants against various pathogens. ⁵⁴ In this study, 20 differentially expressed transcription factors at 3 time points were identified (Figure 3). The number of differentially expressed transcription factors was significantly increased at 24 h after inoculation, and all 12 transcription factors in the WRKY and ERF families had upregulated expression. At 6 h after inoculation, 2 transcription factors MYC4 (c113516_g3 and c114016_g1) were upregulated. The c113516_g3 was upregulated 97-fold at 6 h after inoculation, but no change in the gene expression was observed at 12 and 24 h after inoculation. At 24 h after inoculation, 7 ERF members were upregulated, including Ethylene-responsive transcription factor 1A (c118561_g1) by 24.3-fold. WRKY transcription factor 33 (c96521_g1) was upregulated by 2.7-fold and 12.1-fold at 12 h and 24 h after inoculation respectively. WRKY transcription factor 33 (c113118_g1) and WRKY transcription factors 2, 70, and 40 were upregulated by 3.7-fold and 8.6-fold at 24 h after inoculation (Figure 3). Therefore, we speculate that WRKY and ERF proteins play important roles in the resistance of L. pumilum to F. oxysporum. The role of MYC in the responses of L. pumilum to F. oxysporum needs to be further investigated.

To validate the results obtained by high-throughput sequencing, 4 genes were selected for further confirmation by qRT-PCR. Our results showed that the 4 genes including UTP-glucose-1-phosphate uridylyltransferase (c95325_g1), alpha-1,4 glucan phosphorylase (c121101_g1), phosphoglucomutase (c116359_g1), and NADP-dependent GAPDH (c105111_g1) showed significantly higher expression in 12 h after inoculation with F. oxysporum (Figure 7). Thus, our qRT-PCR results were consistent with those using Illumina sequencing method.

OPEN IN VIEWER

Figure 7.

Transcriptome sequencing differentially expressed genes data validation by quantitative reverse transcription PCR (qRT-PCR).

In this study, transcriptome sequencing was used to study resistance responses during early stages of L. pumilum infected by lily pathogen F. oxysporum. Biological pathways and transcription factors related to the responses of L. pumilum to F. oxysporum were identified using bioinformatics analysis. Our research greatly improves the current understanding of molecular mechanism in Lilium to resist F. oxysporum.

Supplemental Material