Transcriptome Analysis of Fusarium Root-Rot-Resistant and -Susceptible Alfalfa (Medicago sativa L.) Plants during Plant–Pathogen Interactions

Alfalfa (Medicago sativa L.) is a perennial leguminous forage cultivated globally. Fusarium spp.-induced root rot is a chronic and devastating disease affecting alfalfa that occurs in most production fields. Studying the disease resistance regulatory network and investigating the key genes involved in plant–pathogen resistance can provide vital information for breeding alfalfa that are resistant to Fusarium spp. In this study, a resistant and susceptible clonal line of alfalfa was inoculated with Fusarium proliferatum L1 and sampled at 24 h, 48 h, 72 h, and 7 d post-inoculation for RNA-seq analysis. Among the differentially expressed genes (DEGs) detected between the two clonal lines at the four time points after inoculation, approximately 81.8% were detected at 24 h and 7 d after inoculation. Many DEGs in the two inoculated clonal lines participated in PAMP-triggered immunity (PTI) and effector-triggered immunity (ETI) mechanisms. In addition, transcription factor families such as bHLH, SBP, AP2, WRKY, and MYB were detected in response to infection. These results are an important supplement to the few existing studies on the resistance regulatory network of alfalfa against Fusarium root rot and will help to understand the evolution of host–pathogen interactions.


Introduction
Alfalfa (Medicago sativa L.), known as the "queen of forage", has become one of the most important forage crops globally because of its high yield, high quality, satisfactory palatability, and wide adaptability. Alfalfa root rot caused by Fusarium spp. is a major root disease affecting alfalfa production. After infection, the root neck and main root xylem gradually decay and become hollow, deteriorating the ability of the roots to absorb nutrients and water, which causes the plants to gradually die [1].
Fusarium root rot is a soil-borne disease that commonly occurs in alfalfa-growing areas. The disease was first reported in the USA in 1937 and has since been reported in Canada, New Zealand, Australia, India, Egypt, and Japan, with incidence rates of over 60% in some areas [2][3][4][5][6]. Similarly, Fusarium root rot commonly affects alfalfa production in the northern region of China. Fusarium spp. can cause disease throughout the whole growth period of alfalfa. The disease reduces the ability of the plant to fix nitrogen and the quality of alfalfa, greatly shortening its utilization timeframe; it can also cause alfalfa to lose its processing value [7]. Alfalfa consumes stored organic matter during overwintering, and the soluble sugar content of Fusarium-infested alfalfa plants becomes significantly reduced, resulting in the number of tillers being reduced in the following year, which seriously affects yield and quality [7]. Additionally, Fusarium spp. can survive in the soil for a long time and

Experimental Design
The experimental design is shown in Figure 1. After two rounds of disease resistance screening, one resistant and one susceptible clonal line of alfalfa "AmeriGraze401 + Z" were used for transcriptome sequencing analyses.

Screening of Susceptible and Resistant Alfalfa Plants
Alfalfa is a highly heterozygous cross-pollinated plant, and individual alfalfa plants have different genotypes within the variety. To overcome the differences caused by these differing genotypes, clonal lines of resistant and susceptible individual plants were screened as plant materials for this study. Our previous study [19] indicated that "Amer-iGraze401 + Z"was a cultivar that was more resistant to Fusarium root rot based on an evaluation of disease resistance. The pathogen used in this experiment was F. proliferatum strain L1, which was collected and preserved in our laboratory [34]. Two rounds of screening were conducted to ensure the accuracy of screening resistant and susceptible individual plants. The first round entailed a preliminary soil culture screening method [35], followed by hydroponics for the second round of screening [36]. For the first round of screening, alfalfa seeds were sterilized in 25% bleach (v/v) for 15 min, followed by rinsing three times with sterilized water. Next, 1000 sterilized seeds were planted in a sterile soil mixture (soil:vermiculite:perlite = 3:2:1 (w/w/w)) with five seeds per pot. After the alfalfa seeds grew for 30 d, they were inoculated with sterile sorghum (Sorghum bicolor) seeds containing F. proliferatum and cultivated in a greenhouse (25 ± 2 °C day/20 ± 2 °C night, with 75-80% relative humidity, 16 h light/8 h dark). After inoculation 45 d, resistant plants with small root disease spots and susceptible plants with heavy root rot were screened and then propagated using the hydroponic cutting method [37]. After cutting, the plants were rooted for 20 d in a growth chamber under controlled conditions (25 ± 1 °C day/20 ± 1 °C night, with 80% relative humidity and 16 h light (200 mol/m 2 s)/8 h dark). Uniform clonal line plants were inoculated with 5 × 10 6 spores/mL of F. proliferatum L1 under hydroponic cultivation. After 14 d, one of the most resistant and one of the most susceptible plants

Screening of Susceptible and Resistant Alfalfa Plants
Alfalfa is a highly heterozygous cross-pollinated plant, and individual alfalfa plants have different genotypes within the variety. To overcome the differences caused by these differing genotypes, clonal lines of resistant and susceptible individual plants were screened as plant materials for this study. Our previous study [19] indicated that "AmeriGraze401 + Z" was a cultivar that was more resistant to Fusarium root rot based on an evaluation of disease resistance. The pathogen used in this experiment was F. proliferatum strain L1, which was collected and preserved in our laboratory [34]. Two rounds of screening were conducted to ensure the accuracy of screening resistant and susceptible individual plants. The first round entailed a preliminary soil culture screening method [35], followed by hydroponics for the second round of screening [36]. For the first round of screening, alfalfa seeds were sterilized in 25% bleach (v/v) for 15 min, followed by rinsing three times with sterilized water. Next, 1000 sterilized seeds were planted in a sterile soil mixture (soil:vermiculite:perlite = 3:2:1 (w/w/w)) with five seeds per pot. After the alfalfa seeds grew for 30 d, they were inoculated with sterile sorghum (Sorghum bicolor) seeds containing F. proliferatum and cultivated in a greenhouse (25 ± 2 • C day/20 ± 2 • C night, with 75-80% relative humidity, 16 h light/8 h dark). After inoculation 45 d, resistant plants with small root disease spots and susceptible plants with heavy root rot were screened and then propagated using the hydroponic cutting method [37]. After cutting, the plants were rooted for 20 d in a growth chamber under controlled conditions (25 ± 1 • C day/20 ± 1 • C night, with 80% relative humidity and 16 h light (200 mol/m 2 s)/8 h dark). Uniform clonal line plants were inoculated with 5 × 10 6 spores/mL of F. proliferatum L1 under hydroponic cultivation. After 14 d, one of the most resistant and one of the most susceptible plants were selected and propagated using the hydroponic cutting method [37]. The two clonal lines were used for transcriptome sequencing analyses.

Determination of Sampling Time and Sample Collection for RNA-Seq
In order to explore alfalfa early response genes during compatible interaction with F. proliferatum L1, we examined the expression of defense-related genes at different time points after inoculation using real-time quantitative reverse-transcription PCR (qRT-PCR).
NPR1 and NPR3 have been identified as important disease-resistant-related genes in Arabidopsis thaliana [38]. Additionally, NPR1 plays a significant role in the establishment of systemic acquired resistance (SAR) as well as induced systemic resistance (ISR). SPL15 is an IPA1 homologous gene and has been supposed as an upstream disease-resistant transcription factor in the disease resistance of rice [39]. In our previous study, we also verified that SPL15 participates in the process of alfalfa root rot resistance (unpublished data). Therefore, the relative expression levels (2 −∆∆t ) of the disease-resistant genes, NPR1 and NPR3, and their upstream transcription factor, SPL15, were determined by qRT-PCR at 0, 12, 24, 48, 72 h, and 7 d after inoculation with F. proliferatum L1. The inoculation method refers to the hydroponic inoculation previously established by our laboratory [34]. According to the expression level of NPR1, NPR3, and SPL15 genes, the sampling time was determined for transcriptome sequencing.
The two clonal lines prepared in Section 2.2 were inoculated with F. proliferatum L1 for transcriptome sequencing analyses with three replicates. The whole roots of the inoculated and uninoculated (control) groups were collected, frozen in liquid nitrogen, and stored at −80 • C for transcriptome sequencing.
Raw sequences were qualitatively controlled (short fragments and low-quality fragments were removed) to obtain high-quality clean sequences. The M. sativa cv. Xinjiangdaye genome was used as the reference sequence [40]. To ensure precise alignment, the raw RNA sequence data were removed with a 3-end adapter using Cutadapt (https://cutadapt.readthedocs.io/en/stable/ (accessed on 1 January 2022)), and then reads with average quality below Q20 and minimum read size (50 bp) were filtered out. The filtered reads were used for sequence alignment to the reference genome using HISAT2 software (https://daehwankimlab.github.io/hisat2/ (accessed on 3 April 2022)). The read count value was determined by HTSeq (https://htseq.readthedocs.io/en/master/ (accessed on 24 February 2022)).

Differential Gene Expression Analysis of Alfalfa
In order to better exploit resistance-related genes against Fusarium root rot in alfalfa, we analyzed four types of DEGs: (1) the DEGs between resistant (R) and susceptible (S) clonal lines before inoculation (control, C) at each time point (CR24 vs. CS24, CR48 vs. CS48, CR72 vs. CS72, CR7d vs. CS7d); (2) the DEGs between uninoculated (C) and inoculated (Treatment, T) of resistance clonal line at each time point (CR24 vs. TR24, CR48 vs. TR48, CR72 vs. TR72, CR7d vs. TR7d); (3) the DEGs between uninoculated and inoculated groups of susceptible clonal lines at each time point (CS24 vs. TS24, CS48 vs. TS48, CS72 vs. TS72, CS7d vs. TS7d); (4) the DEGs between resistance and susceptible clonal lines after inoculation at each time point (TR24 vs. TS24; TR48 vs. TS48; TR72 vs. TS72, TR7d vs. TS7d). DESeq 2 (DESeq2 R package 1.16.1) was used to screen the DEGs, and the screening conditions were p-adjusted < 0.05 and |log 2 FC| ≥ 1. The differential gene expression sequences were annotated by NR, IPR, TREMBL, and Swiss-Prot based on the reference genome. Clonal line type-specific DEGs were obtained by removing the genetically different genes DEGs generated in uninoculated group from all the DEGs induced by F. proliferatum L1 inoculation using FunRich software (http://www.funrich. org/ (accessed on 6 March 2022)). Gene ontology (GO) enrichment analysis of the DEGs was performed using topGO, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (p-value ≤ 0.05) were conducted using clusterProfiler according to biological process (BP), molecular function (MF), and cellular component (CC) among all up-and downregulated genes at each time point after inoculation.

Differential Gene Expression Analysis of F. proliferatum L1
The number and expression level of F. proliferatum L1 genes can be detected in the inoculation group of the two clonal lines. The F. proliferatum ET1 genome was used as the reference sequence for sequence alignment, transcriptome, and alignment. GO and KEGG enrichment analyses were performed on detected genes.

Real-Time Quantitative PCR (qRT-PCR) Analysis of the DEGs
In order to verify the expression patterns of genes observed in the RNA-seq analysis, we used the total RNA of all the samples subjected to transcriptome sequencing for qRT-PCR verification. A 1 µg sample of total RNA was reverse-transcribed into cDNAs using the PrimeScript TM RT reagent kit with gDNA Eraser (TaKaRa, Kusatsu, Japan). A NanoDrop One (Thermo Fisher Scientific Inc., Waltham, MA, USA) was used to quantify 100 ng of purified single-stranded cDNA for qRT-PCR. Relative quantitative analyses were performed in a CFX-96 Touch Real-time PCR detection system (Bio-Rad, Inc., Hercules, CA, USA) with the following cycling conditions: 95 • C for 30 s and 40 cycles at 95 • C for 5 s, 60 • C for 30 s.
Three technical replicates were performed for each sample. Gene-specific primers were designed using Primer Premier 5 software (PREMIER Biosoft, San Francisco, USA) (Table S1). Relative quantification was normalized to the housekeeping control gene (β-actin), and the fold change (FC) in gene expression was calculated using the 2 −∆∆t method.

Screening of Experimental Plants and Determination of Sampling Time
In the first screening, 58 susceptible and 42 resistant plants with obvious traits were selected for the second screening (Figure 2A,B). After propagation by the hydroponic cutting method and F. proliferatum L1 inoculation under hydroponic conditions, one of the most resistant and one of the most susceptible clonal lines ( Figure 2C,D) with obvious symptoms were screened after surface sterilization to generate uniform clonal lines by cutting propagation for transcriptome sequencing.
Finally, qRT-PCR was performed to determine the relative expression levels (2 −∆∆t ) of the disease-resistance-related genes, NPR1 and NPR3, and their upstream transcription factor SPL15 at six time points. The expression levels of NPR1 and SPL15 were significantly increased in the inoculated group compared to those in the uninoculated group, although the NPR3 gene was only increased slightly at 24 h ( Figure 2E). No obvious disease symptoms were found after inoculation until 7 d, when the plants exhibited yellowing or falling ( Figure 2F). Based on the qRT-PCR results and post-inoculation phenotype, we collected samples at 24 h post-inoculation (hpi), 48 hpi, 72 hpi, and 7 days post-inoculation (dpi) for transcriptome sequencing.

Alfalfa Sequence Analysis and Alignment with the Reference Genome
For each sample, approximately 8 Gb of clean sequences were obtained. The clean reads% and clean data% for the library transcriptome sequencing obtained from uninoculated or inoculated resistant and susceptible clonal lines were both higher than 90% (Table S2). Finally, qRT-PCR was performed to determine the relative expression levels (2 −∆∆t ) of the disease-resistance-related genes, NPR1 and NPR3, and their upstream transcription factor SPL15 at six time points. The expression levels of NPR1 and SPL15 were significantly increased in the inoculated group compared to those in the uninoculated group, although the NPR3 gene was only increased slightly at 24 h ( Figure 2E). No obvious disease symptoms were found after inoculation until 7 d, when the plants exhibited yellowing or falling ( Figure 2F). Based on the qRT-PCR results and post-inoculation phenotype, we collected samples at 24 h post-inoculation (hpi), 48 hpi, 72 hpi, and 7 days post-inoculation (dpi) for transcriptome sequencing.

Alfalfa Sequence Analysis and Alignment with the Reference Genome
For each sample, approximately 8 Gb of clean sequences were obtained. The clean reads% and clean data% for the library transcriptome sequencing obtained from uninoculated or inoculated resistant and susceptible clonal lines were both higher than 90% (Table S2).
On average, 85% of the total sequences were mapped to the M. sativa cv. Xinjiangdaye reference genome sequence, with most of the sequences being mapped to exons-approximately 96%-and only a small portion of them being mapped to intergenic regions, approximately 17% (Table S2). To evaluate the similarity of RNA-Seq data from different samples, the mapped RNA-Seq data from different samples were used for a principal component analysis (PCA) with DEseq2 ( Figure S1). In addition, Pearson's correlation among biological replicates ( Figure S2) for all the samples was analyzed and showed a On average, 85% of the total sequences were mapped to the M. sativa cv. Xinjiangdaye reference genome sequence, with most of the sequences being mapped to exons-approximately 96%-and only a small portion of them being mapped to intergenic regions, approximately 17% (Table S2). To evaluate the similarity of RNA-Seq data from different samples, the mapped RNA-Seq data from different samples were used for a principal component analysis (PCA) with DEseq2 ( Figure S1). In addition, Pearson's correlation among biological replicates ( Figure S2) for all the samples was analyzed and showed a high correlation between sequencing replicates (approximately 0.95). The results of PCA and Pearson's correlation indicated there is a high correlation between sequencing replicates.

DEGs between Resistant and Susceptible Clonal Lines in Uninoculated or Inoculated Conditions
The numbers of DEGs between resistant and susceptible clonal lines in uninoculated or inoculated conditions were analyzed at four time points ( Figure 3A).
To eliminate the effects of genetic differences, DEGs between resistant and susceptible clonal lines at the same time points were identified before inoculation (CR 24 vs. CS24, CR48 vs. CS48, CR72 vs. CS72, and CR7d vs. CS7d) (|log 2 FC| ≥ 1 and FDR < 0.05). There were 1057 overlapped DEGs detected at four time points ( Figure 3B). Among the 1057 DEGs, the gene expression level of 592 DEGs at 24 hpi, 748 DEGs at 48 hpi, 741 DEGs at 72 hpi, and 607 DEGs at 7 dpi in resistant plants were higher than that in susceptible clonal line. Additionally, the expression levels of 380 DEGs in the resistant clonal line were always higher than that in the susceptible clonal line at four time points, but the expression levels of only 237 DEGs in the susceptible clonal line were higher than that in the resistant clonal line. In addition, more DEGs were detected between the resistant and susceptible clonal lines at 24 h (CR24 vs. CS24). GO enrichment analysis and KEGG enrichment analysis were conducted for these 1057 DEGs, and these genes were enriched in 712 biological processes and 86 KEGG pathways (Table S3). Among the 712 biological processes, the top five biological processes with more DEGs were biological, metabolic, cellular, organic substance metabolic, and primary metabolic processes. Among the 86 KEGG pathways, the top five KEGG pathways with more DEGs were metabolic, biosynthesis of secondary metabolites, ribosome, carbon metabolism, and biosynthesis of amino acids pathways. We hypothesized that there might be differences in the existence of structural, physical, or chemical barriers between resistant and susceptible clonal lines based on DEG analysis before inoculation (CS vs. CR).
es 2022, 13, x FOR PEER REVIEW 7 of

DEGs between Resistant and Susceptible Clonal Lines in Uninoculated or Inoculated Conditions
The numbers of DEGs between resistant and susceptible clonal lines in uninoculat or inoculated conditions were analyzed at four time points ( Figure 3A). To eliminate the effects of genetic differences, DEGs between resistant and suscep ble clonal lines at the same time points were identified before inoculation (CR 24 vs. CS CR48 vs. CS48, CR72 vs. CS72, and CR7d vs. CS7d) (|log2FC| ≥ 1 and FDR < 0.05). Th were 1057 overlapped DEGs detected at four time points ( Figure 3B). Among the 10 DEGs, the gene expression level of 592 DEGs at 24 hpi, 748 DEGs at 48 hpi, 741 DEGs 72 hpi, and 607 DEGs at 7 dpi in resistant plants were higher than that in susceptible clo line. Additionally, the expression levels of 380 DEGs in the resistant clonal line were ways higher than that in the susceptible clonal line at four time points, but the expressi levels of only 237 DEGs in the susceptible clonal line were higher than that in the resista clonal line. In addition, more DEGs were detected between the resistant and suscepti clonal lines at 24 h (CR24 vs. CS24). GO enrichment analysis and KEGG enrichment an ysis were conducted for these 1057 DEGs, and these genes were enriched in 712 biologi processes and 86 KEGG pathways (Table S3). Among the 712 biological processes, the t Among the 1057 DEGs, we focus on the genes enriched in biological processes related to defense, such as the response to stress, cell wall composition, and the genes enriched in defense-related pathways, such as the biosynthesis of secondary metabolites, plant-pathogen interaction, plant hormone signal transduction, and the MAPK signaling pathway.
A total of 27 DEGs among the 1057 DEGs were enriched in biological processes related to defense. We  Figure 3C). More DEGs were detected between the resistant and susceptible clonal lines at 48 h (TR48 vs. TS48), and the top five biological processes with more DEGs were structural molecule activity, structural constituent of ribosome, small molecule binding, nucleotide binding, and nucleoside phosphate binding. Similarly, it was found that GO terms and KEGG of the top 20 observed in the TR48 vs. TS48 were almost the same as in the CR24 vs. CS24 pair ( Figure S3).

Transcriptional Changes in Response to F. proliferatum L1 Inoculation
Resistant and susceptible clonal lines inoculated with F. proliferatum L1 were compared with the uninoculated group, and the DEGs that responded to the pathogen infection in each clonal line at four time points were obtained (|log 2 FC| ≥ 1 and FDR < 0.05) (Table S4, Figure 4). The number of DEGs between the inoculated and uninoculated resistant clonal lines was much less than that of the susceptible clonal line at the early stage of inoculation (24,48, and 72 hpi) ( Figure 4A,C). In particular, there were more DEGs in the susceptible line in the early stage (24 h), but there were more DEGs in the resistant line at the later stage (7 dpi). We analyzed the top 30 genes with a high fold change between uninoculated and inoculated groups of each line at four time points (Table S5). The top 30 genes with a high fold change compared with those of the control in the two clonal lines at 24 hpi and 7 dpi are shown in Figure 4. At 24 h, we found that 19 genes of the 30 genes were upregulated in the resistant line, while all the genes were downregulated in the susceptible line. Additionally, there was no overlap in the DEGs between the two lines. GO enrichment analyses of these DEGs indicated enrichment in small molecule binding, catalytic activity, nucleotide binding, nucleoside phosphate binding, and anion binding. KEGG enrichment analyses of these DEGs indicated enrichment in glycolysis/gluconeogenesis, citrate cycle (TCA cycle), fatty acid degradation, pyruvate metabolism, alanine, aspartate and glutamate metabolism, phenylpropanoid biosynthesis, caffeine metabolism, and the pentose phosphate pathway ( Figure S4)

Common Transcriptional Changes in Response to the Inoculation of Resistant and Susceptible Clonal Lines
A total of 2078 overlapping DEGs between the uninoculated and inoculated groups were identified in both the clonal lines ( Figure 5, Table S6). We performed a functional enrichment analysis for GO terms and found that most of these genes were enriched in the polysaccharide catabolic process, polysaccharide biosynthetic process, reproduction, and the response to acid chemical processes ( Figure S5). KEGG enrichment analyses were performed for these DEGs, and after excluding genes with unspecified functions, common genes related to metabolic pathways (50.4%) accounted for the largest proportion, followed by genes related to plant-pathogen interactions (43%), including the categories: biosynthesis of secondary metabolites (29.7%), plant hormone signal transduction (5.9%), phenylpropanoid biosynthesis (5.1%), plant-pathogen interaction (4.3%), and the MAPK signaling pathway (3.1%) ( Figure S6).

Common Transcriptional Changes in Response to the Inoculation of Resistant and Susceptible Clonal Lines
A total of 2078 overlapping DEGs between the uninoculated and inoculated gr were identified in both the clonal lines ( Figure 5, Table S6). We performed a functi enrichment analysis for GO terms and found that most of these genes were enriche the polysaccharide catabolic process, polysaccharide biosynthetic process, reproduc and the response to acid chemical processes ( Figure S5). KEGG enrichment analyses performed for these DEGs, and after excluding genes with unspecified functions, com genes related to metabolic pathways (50.4%) accounted for the largest proportion lowed by genes related to plant-pathogen interactions (43%), including the catego biosynthesis of secondary metabolites (29.7%), plant hormone signal transduction (5 phenylpropanoid biosynthesis (5.1%), plant-pathogen interaction (4.3%), and the M signaling pathway (3.1%) ( Figure S6).  There were many transcription factors (TFs) in these DEGs that were shown to be involved in resistance to pathogens. They included bHLH (BHLHW, BHLH123, BHLH30, BBD1, BHLH25, BHLH80, and BHLH126), AP2 (BBM2, ANT, and WRI1), SBP (ALPL and SPL6), WRKY (WRKY23, WRKY75, WRKY6, WRKY51, WRKY48, WRKY70, WRKY40, WRKY1, and WRKY42), and MYB (MYB87, MYB78, MYB108, MYB88, and MYB4) TF families. These plant-pathogen interaction genes and disease-resistant-related transcription factors showed roughly the same trend at the four time points in both clonal lines after inoculation compared to the uninoculated plants ( Figure 6).

Clonal Line Type-Specific Transcriptional Changes in Response to Inoculation
After removing the DEGs of the genetically different genes that were generated in the uninoculated group from all the DEGs induced by F. proliferatum L1 inoculation (Fig-Figure 6.

Clonal Line Type-Specific Transcriptional Changes in Response to Inoculation
After removing the DEGs of the genetically different genes that were generated in the uninoculated group from all the DEGs induced by F. proliferatum L1 inoculation ( Figure 4A,C) using FunRich software, the remaining DEGs were treated as transcriptional changes in response to F. proliferatum L1 inoculation at four time points ( Figure 5, Table S6). Specific transcriptional changes induced by F. proliferatum L1 included 3512 and 4081 genes found to be differentially regulated in the resistant and susceptible clonal lines, respectively ( Figure 5). The most prevalent functional categories among the modulated genes were those of metabolic process (41.5% and 43.1%, respectively) and the biosynthesis of secondary metabolites (27.4% and 24.5%, respectively) for both clonal lines. Other categories of plant-pathogen interaction genes included plant-pathogen interactions (5.6% and 6.9%, respectively) and phenylpropanoid biosynthesis (4.7% and 3.3%, respectively). Additionally, plant hormone signal transduction (3.7% and 5.6%, respectively) and the MAPK signaling pathway (2.6% and 4.4%, respectively) were prevalent in both clonal lines. These percentages were approximately consistent with the percentage of DEGs associated with plant-pathogen interactions found in the common transcriptional changes in both clonal lines ( Figure S7).
In the resistant clonal line, 485 clonal line type-specific DEGs were detected at 24 hpi ( Figure 7). Among these genes, 212 were downregulated, and 273 were upregulated. Six genes (MS.gene65772, MS.gene38785, MS.gene067433, MS.gene011203, MS.gene00935, and MS.gene 055695) of the top ten DEGs with a high fold change were upregulated. These six genes were annotated as a transmembrane protein, methyltransferase DDB, protein nuclear fusion defective 4 isoform X1, cytochrome P450 88D2, epidermal patterning factor protein 2, and receptor-like protein 18, respectively. Four genes (MS.gene98474, MS.gene013794, MS.gene 038488, and MS.gene71722) of the top ten DEGs were downregulated. These four genes were annotated as LRR and NB-ARC domain disease-resistance protein, metal transporter Nramp5, P-loop nucleoside triphosphate hydrolase superfamily protein, and transmembrane protein, respectively. In susceptible clonal line, five genes, including MS.gene26169 (protein reveille 1), MS.gene43364 (polygalacturonase inhibitor), MS.gene37066 (MYB family transcription factor PHL5), MS.gene08790 (plasmodesmata callose-binding protein 3), and MS.gene66665 (aluminum-activated malate transporter 2 isoform X1), were upregulated at all four time points. Additionally, TF families, including bHLH, AP2, MYB, WRKY, and SBP, were analyzed in all the DEGs induced by F. proliferatum L1 in the resistant and susceptible clonal lines, respectively. These TF families exhibited different expression trends in specific and common DEGs, but these DEGs changed significantly during the early stage in the susceptible clonal line and changed significantly during the late stage in the resistant clonal line (Table S7).

F. proliferatum L1 Genes Were Enriched in Resistant and Susceptible Clonal Lines after Inoculation
Using F. proliferatum ET1 as the reference genome, the number of fungal genes detected in the inoculation groups was counted. GO and KEGG enrichment analyses were performed on these genes, and the results indicated that they were mostly low-molecular-weight secondary metabolites such as lipids, glycosides, polysaccharides, peptides, proteins, and glycoproteins, which are typical pathogenic factors of host-selective fungi. The patterns of expression of fungal genes were distinct during the interactions between the two lines. We found that the number of fungal genes detected in the susceptible clonal line was significantly higher (one-way analysis of variance (ANOVA), p < 0.001) than that detected in the resistant clonal line during the early stage. Although the number of fungal genes detected in the susceptible line at 24 h was slightly higher than that of the resistant line, they were consistent at the other time points, and the top 30 genes expressed at the highest levels detected in the two lines at the four time points were basically consistent (Table S8, Figure 8).
MS.gene37066 (MYB family transcription factor PHL5), MS.gene08790 (plasmodesmata callose-binding protein 3), and MS.gene66665 (aluminum-activated malate transporter 2 isoform X1), were upregulated at all four time points. Additionally, TF families, including bHLH, AP2, MYB, WRKY, and SBP, were analyzed in all the DEGs induced by F. proliferatum L1 in the resistant and susceptible clonal lines, respectively. These TF families exhibited different expression trends in specific and common DEGs, but these DEGs changed significantly during the early stage in the susceptible clonal line and changed significantly during the late stage in the resistant clonal line (Table S7).

F. proliferatum L1 Genes Were Enriched in Resistant and Susceptible Clonal Lines after Inoculation
Using F. proliferatum ET1 as the reference genome, the number of fungal genes detected in the inoculation groups was counted. GO and KEGG enrichment analyses were performed on these genes, and the results indicated that they were mostly low-molecularweight secondary metabolites such as lipids, glycosides, polysaccharides, peptides, proteins, and glycoproteins, which are typical pathogenic factors of host-selective fungi. The patterns of expression of fungal genes were distinct during the interactions between the two lines. We found that the number of fungal genes detected in the susceptible clonal line was significantly higher (one-way analysis of variance (ANOVA), p < 0.001) than that detected in the resistant clonal line during the early stage. Although the number of fungal genes detected in the susceptible line at 24 h was slightly higher than that of the resistant line, they were consistent at the other time points, and the top 30 genes expressed at the highest levels detected in the two lines at the four time points were basically consistent (Table S8, Figure 8).

Validation of RNA-Seq Data by qRT-PCR
To validate the RNA-Seq expression profiles of DEGs, qRT-PCR was performed on three disease-associated genes (NPR1, PR1, and PR4A), eight disease-associated transcription factors (WRKY22, WRKY29, WRKY33, SPL6, SPL7, SPL9, SPL15, and SPL16), and MS.gene03877 (NBS-LRRtype disease resistance protein). For the relative expression levels of these 12 genes in the resistant and susceptible clonal lines and inoculated relative to the uninoculated group at four time points, the qRT-PCR results were approximately consistent with those of the RNA-Seq data. However, the RNA-Seq analysis displayed a higher dynamic range, showing larger differences between the fold change compared with qRT-PCR (Figure 9).

Validation of RNA-Seq Data by qRT-PCR
To validate the RNA-Seq expression profiles of DEGs, qRT-PCR was performed on three disease-associated genes (NPR1, PR1, and PR4A), eight disease-associated transcription factors (WRKY22, WRKY29, WRKY33, SPL6, SPL7, SPL9, SPL15, and SPL16), and MS.gene03877 (NBS-LRRtype disease resistance protein). For the relative expression levels of these 12 genes in the resistant and susceptible clonal lines and inoculated relative to the uninoculated group at four time points, the qRT-PCR results were approximately consistent with those of the RNA-Seq data. However, the RNA-Seq analysis displayed a higher dynamic range, showing larger differences between the fold change compared with qRT-PCR ( Figure 9).

Discussion
In this study, we generated a unique set of DEGs induced by F. proliferatum L1 in alfalfa. Overall, 85% of the total sequences were mapped to the M. sativa cv. Xinjiangdaye reference genome sequence. This finding is significant for understanding the temporal expression of genes in alfalfa after inoculation with F. proliferatum L1. It is also can be used as a candidate gene pool against Fusarium spp.
The number of DEGs in the susceptible clonal line between the inoculated and uninoculated groups was twice that of the resistant clonal line between inoculated and uninoculated groups at 24 hpi. Many of these DEGs overlapped between the two uninoculated clonal lines. We hypothesize that this may be due to the genetic differences between the two clonal lines in their physical structure and chemical barrier against F. proliferatum. In our study, 17 respiratory burst oxidase homolog protein (RBOH) genes were detected. After inoculation, the change trends of these 17 genes in the two clonal lines were basically the same, but the differential expression in the susceptible clonal line occurred earlier than that in the resistant clonal line, possibly due to the pathogen invading the susceptible clonal line earlier and initiating oxidative activity. Notably, the miraculin gene was the most highly expressed gene at 24 h after inoculation in the resistant line. Previous researchers have focused on the benefits of miraculin to the human and have not reported on the function of miraculin on the plant itself [41,42]; however, a new study showed that miraculin might primarily play a role in regulating seed germination and maturation, resisting pathogen infection and environmental stress, and regulating plant growth [43]. These results indicate that the peculiar property of miraculin that modifies sour tastes to sweet tastes may be secondary, and the main meaning of its existence is to benefit itself. A total of eight HSP90 genes were detected, among which three DEGs (MS.gene30654, MS.gene38210, and MS.gene49154) were shared in two clonal lines after inoculation but were expressed earlier in the resistant line than in the susceptible line. Three DEGs (MS.gene24136, MS.gene58998, and MS.gene81285) were unique in the resistant line, and two DEGs (MS.gene80872 and MS.gene81288) were unique to the susceptible line after inoculation (Table S9). These results suggest that the early expression of HSP90 may help plants to remove ROS in the resistant line.
The difference in the plant immunity between the two clonal lines may also be the reason for their different defense response levels to the pathogens. Plants have evolved immunity to pathogens via either PTI or ETI [17,18]. Plants possess PRRs on their surface membrane that trigger PTI by PAMPs and belong to families of receptor-like kinases (RLKs) [44,45]. According to KEGG analyses, the DEGs of LRR receptor-like serine/threonine-protein kinase FLS2, LRR receptor-like serine/threonine-protein kinase EFR, receptor kinase-like protein Xa21, probable LRR receptor-like serine/threonineprotein kinase, and putative receptor-like protein kinase were involved in the alfalfa immunity. Additionally, some leucine-rich repeats (LRR) receptor-like serine/threonineprotein kinases were up-or downregulated after inoculation in one or both clonal lines. After inoculation, the expression of these genes in the resistant clonal line was significantly upregulated compared to that in the uninoculated group, whereas those in the susceptible clonal line were almost unchanged (Table S10).
Downstream signaling networks are triggered by pathogen recognition and mediated by protein kinases (PK), particularly calcium-activated protein kinases (CIPKs or CDPK) and mitogen-activated protein kinase (MAPK), which control defense responses [46][47][48]. Calcium-dependent protein kinases (CDPKs) comprise a large family of serine/threonine kinases in plants and protozoans [48]. In our data, no CDPK genes were detected in the common transcriptional changes in response to inoculation and in the DEGs specific to the resistant clonal line, whereas three CDPK genes were detected in the DEGs specific to the susceptible clonal line, compared with the inoculation group. The three CDPK genes (MS. gene41297, MS.gene50367, and MS.gene60795) increased slightly during the early stage and changed significantly at 72 h. The expression of MS. gene41297 decreased, while those of MS.gene50367 and MS.gene60795 increased (Table S11). MAPK cascades are highly conserved signaling modules downstream of receptors/sensors that transduce extracellular stimuli into intracellular responses in eukaryotes. Plant MAPK cascades play pivotal roles in signaling plant defenses against pathogen attacks. MAPK cascades have also emerged as battlegrounds for plant-pathogen interactions. In Arabidopsis, pathogen perception by plants leads to the activation of two independent MAPK cascades [49,50]: the MEKK1-MKK4/MKK5-MPK3MPK6 cascade and the MEKK1-MKK1/MKK2-MPK4 cascade. Two MEKK1 genes (MS.gene73935 and MS.gene94111) were detected in the DEGs of the resistant clonal line after inoculation. Both genes were slightly upregulated or downregulated in the early stages. However, both were significantly upregulated at 7 d.
In our study, we found that 18 CERK1 genes were expressed, and the expression levels of these genes were basically the same in the two clonal lines without inoculation, but expression levels in the resistant clonal line were slightly higher than those in the susceptible clonal line. The expression level of these genes slightly increased in the inoculated resistant clonal line compared to the uninoculated resistant clonal line, whereas there was a decrease in the susceptible clonal line during the early (24 hpi) and late (7 dpi) stages after inoculation and an increase in the middle stages (48 hpi and 72 hpi) (Table S12).
Downstream defense-responsive genes are normally positively or negatively regulated by TFs and are direct or indirect targets of various signal-transduction pathways. The TF families found in our study are widely reported to be involved in plant defense responses, including MYB, WRKY, ethylene-responsive factors (ERFs), squamosa promoter-binding protein-like (SPL), zinc finger domain proteins, and basic helix-loophelix (bHLH). The basic helix-loop-helix (bHLH) transcription factor family is one of the largest transcription factor gene families in Arabidopsis thaliana, and these transcription factors have the pleiotropic regulatory roles in plant growth and development, stress response, biochemical functions, and the web of signaling networks [51]. These common disease-resistancerelated transcription factors in response to inoculation showed approximately the same trend at the four time points in both lines between the inoculated and uninoculated groups. However, type-specific differentially expressed TFs of changes in the clonal line in response to inoculation were less in the resistant clonal line than in the susceptible clonal line, but there was greater variation in the resistant clonal line.
The phenylpropanoid pathway generates secondary metabolites such as lignin, flavonoids, and phytoalexins, which are involved in plant defense against pathogens [52]. Flavonoid biosynthesis begins with the amino acid phenylalanine, and the terminal products include anthocyanins, flavones, isoflavones, and condensed tannins. Phytoalexins such as pisatin and lignans are well-known defense metabolites because of their potent antifungal activities and their ability to inhibit secreted fungal enzymes [53]. In our study, secondary metabolism was the category with the largest number of enriched genes in addition to metabolic processes ( Figure S4). Biosynthesis of phenylpropanoids (MS.gene71 342, MS.gene29407, MS.gene08187, MS.gene44869, etc.), flavonoids (MS.gene75693, MS.gene 89654, and MS.gene02599), lignin (MS.gene04038 and MS.gene66786), and anthocyanins (MS.gene00229, MS.gene05556, and MS.gene02621) were detected in both clonal lines after inoculation ( Figure S4).
The Fusarium genus contains filamentous ascomycete fungi that can infect a diverse range of plants, and a large number of similar phytopathogenic genes undergo diverse selection during host-pathogen interactions [54]. Evolution has equipped Fusarium pathogens with a wide variety of infection strategies [55]. These include the production and secretion of proteins and other effectors to successfully facilitate the infection process by reprogramming the host metabolism and parasitic colonization by manipulating the host cell's immune response [56]. In the inoculation group, we detected the genes related to lipids (gene-FPRO_00503 and gene-FPRO_01999), glycosides, polysaccharides, peptides (gene-FPRO_03131), and proteins (gene-FPRO_00137, gene-FPRO_13465, and gene-FPRO_ 03476) of the fungus that recognized the different secreted effectors of Fusarium pathogens (Table S6).

Conclusions
To the best of our knowledge, this is the first study to systematically study the defense transcriptome and mining-resistant-related genes response to Fusarium spp. infection in alfalfa. In summary, DEGs were identified in control samples at 24 h, 48 h, 72 h, and 7 d after inoculation with F. proliferatum L1 in resistant and susceptible clonal lines of alfalfa. Importantly, we studied plant-pathogen interaction genes, plant hormone signal transduction, the MAPK signaling pathway, secondary metabolism, multiple disease resistance proteins, TFs, and genes involved in cell wall expansion and antioxidant processes that were modulated by inoculation with F. proliferatumL1. Overall, this study extended our understanding of the d molecular defense of two clonal lines with different genetic backgrounds during F. proliferatum inoculation. We identified several candidate genes that could be useful for future research and expect that these data will provide valuable information for research on alfalfa root rot.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13050788/s1. Figure Figure S5: GO enrichment analyses of overlapping DEGs in response to the inoculation of resistant and susceptible clonal lines. Figure S6: KEGG enrichment analyses performed for overlapping DEGs in response to the inoculation of resistant and susceptible clonal lines. Figure S7: Comparison of KEGG enrichment analyses of clonal line-specific differentially expressed genes in two clonal lines. Table S1: Real-time RT-PCR primer sequences. Table S2: Project sample information, raw data, quality control data, and comparative data. Table S3: Analysis of overlapped differentially expressed genes between two uninoculated groups at four time points. Table S4: DEGs in resistant and susceptible clonal lines in response to F. proliferatum. Table S5: The top 30 DEGs in the inoculated compared to uninoculated groups of the two lines with differential fold change at four time points. Table S6: The DEGs in response to F. proliferatum L1 inoculation. Table S7: Clonal line type-specific DEGs transcription factor family. Table S8: Top 30 fungal genes detected in the two lines at four time points. Table S9: Expression of respiratory burst oxidase homolog protein (RBOH) and HSP90 in both lines. Table S10: LRR receptor-like serine/threonine-protein kinase in both clone lines. Table S11: Expression difference of three CDPK genes in susceptible clone line at four time points between inoculated and uninoculated groups. Table S12: Expression of CERK1 genes in both clone lines.  Data Availability Statement: The original contributions presented for this study can be found in the NCBI database (accession number: PRJNA782344; https://www.ncbi.nlm.nih.gov/bioproject/ ?term=PRJNA782344 (accessed on 1 April 2022)).