Differential gene expression analysis of resistant and susceptible flax cultivars in response to Fusarium oxysporum stress

A plant’s early response to pathogen stress is a vital indicator of its disease resistance. In order to study the response mechanism of resistant and susceptible flax cultivars to Fusarium oxysporum f. sp. lini ( Foln ), we applied RNA-sequencing to analyze transcriptomes of flax with Foln 0.5, 2 and 8 hours post inoculation (hpi). We found a significant difference in the number of differential expression genes (DEGs) between resistant and susceptible flax clutivars. The number of DEGs in the Fusarium -resistant cultivar increased dramatically at 2 hpi, and a large number of DEGs participated in the Fusarium -susceptible cultivar response to Foln infection 0.5 hpi. GO enrichment analysis determined that the up-regulated DEGs of both flax cultivars were enriched such as oxidoreductase activity and oxidation-reduction process. At the same time, the genes involved in diterpenoid synthesis were up-regulated in resistant cultivar, while those involved in extracellular region, cell wall and organophosphate ester transport were down-regulated in susceptible cultivar. KEGG enrichment analysis showed the genes encoded WRKY 22 and WRKY33 which involved in MAPK signaling pathway were up-regulated expressed in S-29 and down-regulated expressed in R-7, negatively regulated the disease resistance of flax; The genes encoded Hsp 90 family which in involved in plant pathogen interaction pathway were up-regulated in R-7 and down-regulated in S-29, which positively regulated the disease resistance of flax; The genes encoded MYC2 transcription factor and 2 TIFY proteins which involved in plant hormone signaling pathway were up-regulated in R-7, and regulated the jasmonic acid metabolism of flax and the signal transduction of plant hormones. Meanwhile seven regulatory genes with the most correlation were screened out, Among Lus10025000.g and Lus10026447.g regulated other genes expressed both in plant hormone signal transduction pathway and MAPK signal pathway. In conclusion, these findings will facilitate further studies on the function of these candidate genes in flax of response to Fusarium stress, and the breeding of disease-resistant flax cultivar. of hormone signaling, activation of pathogenesis-related genes, and changes in secondary metabolism. Studies had shown after stressed in Fusarium oxysporum , the content


Introduction
Flax (Linum usitatissimum) is an important crop for the production of fiber, oil, and nutraceuticals. Among flax diseases caused by fungal pathogens, fusarium wilt, caused by Fusarium oxysporum f. sp. lini (Foln) has been an important factor in limiting yield of this plant [1]. Fusarium wilt, one of the most prevalent diseases in flax production, occurs at all stages of flax growth. The typical annual incidence rate of flax fusarium wilt is 10% ～ 30%, and the highest incidence rate is more than 50% [2]. The pathogen of this disease is soilborne pathogenic fungus Foln, which can be effectively controlled by breeding and growing resistant cultivars. Therefore, studying different resistance mechanism responses of resistant and susceptible flax cultivars to Foln infection can also facilitate the breeding of resistant cultivar.
A host's stress response to pathogens is a complex biological process that involves multiple defense strategies and metabolic pathways [3]. When attacked by pathogens, a host activates a large number of gene expression and metabolic pathways to cope with the biotic stress [4]. Plants respond to pathogen attack by modifying their defense gene expression and inducing the production of myriad proteins and metabolites [5]. During the constant 3 interaction between plants and pathogens, a compatible response leads to pathogen invasion and disease development, while an incompatible response ignites the plant's immune system and prevents pathogen invasion. Due to developments in RNA-seq technology, transcriptome profiling has been determined to be an efficient measure when deciphering the plant-pathogen interaction mechanism based on different gene expressions [6]. Transcriptome analysis showed that 12 genes encoding secretory proteins were significantly expressed in rice after being inoculated with Magnaporthe oryzae for 36 h [7]. A large number of differentially expressed genes (DEGs) in melon were involved in the response to Acidovorax citrulli at 6 hours post inoculation (hpi) and 12 hpi, mainly in the plant-pathogen interaction pathway [8]. Liu et al. [9] screened out powdery mildew resistance-related genes, resistance pathways, and transcriptant factors using a transcriptome analysis of resistant and susceptible grape cultivars in the early stages of Plasmopara viticola biotic stress .
Some scholars have actively studied the interaction between flax and pathogens. Lignin metabolism in flax cell suspensions was increased under the pressed of Fusarium oxysporum [10]. The chitinase gene of flax was up-regulated expressed after inoculation of Fusarium oxysporum [1]. Flax inoculated with Fusarium after 48 h, the gene that incoded formate dehydrogenase regulated the pectin content in the cell wall of flax was decreased, thus reduced the resistante flax to Fusarium [11], and content of lignin and pectinin the cell wall of flax were increased after invasion of Fusarium oxysporum which improved the resistance of flax [12]. Galindo-Gonzalez and Deyholos. [1]had for the first study of the flax- Fusarium oxysporum, and showed the number of the genes expressed more in resistant genotypes than in sensitive genotypes, which genes related to biological stimulation and stress response, defense response, antioxidant activity and cell wall tissue or biogenesis of flax [14] . These studies provided theoretical basis for the resistance mechanism of flax, but the transient stress response and gene expression of flax exposed to Fusarium had not been reported so far. Based on previous studies, this study analyzed the gene expression of flax in the transient response to Fusarium infestation, which could clarify the molecular response mechanism of flax resistant cultivars under pathogen stress, and further reveal the disease resistance mechanism of flax, providing a theoretical basis for breeding of resistance flax cultivar.

Plant materials
In the present study, Linum usitatissimum 'Jinya 7' (R-7) and 'Shanxi Huang 29' (S-29), two cultivars of flax that are respectively resistant and sensitive to Fusarium oxysporum were used in the present study. Test materials were provided by the High Latitude Crops Institute, Shanxi Agricultural University, Shanxi Academy of Agricultural Sciences, Shanxi Agricultural University. Intact flax seeds were disinfected in 2% sodium hypochlorite by oscillation at 150 rpm for 10 min, and were washed thrice with sterile water. Seeds were germinated on filter paper soaked with distilled water for 5 days in the dark conditions. For continued hydroponic growth, the seedlings were transferred to falcon tubes with filter paper soaked with half strength sucrose-free Murashige and Skoog (MS) medium and grown in a growth chamber at 21/16 ℃ with 16/8 h light/dark photoperiod to four-leaf stage. the Foln in Petri dishes with potato dextrose agar. To induce sporulation, we cultured the isolate in mung bean liquid medium using oscillation at 150 rpm for 7 days at 25 ℃ [15]. The spores were collected by centrifugation at 10,000 rpm for 10 min at 4 ℃, and then washed thrice with sterile water. The density of collected spores was then adjusted to 1×10 6 spores.ml -1 with half-strength MS medium.

Pathogen inoculation and sampling
Flax seedlings at the four-leaf stage were placed in the pathogen suspension(1×10 6 spores.mL -1 ) and treated in sterile water as a blank control. Each treatment was prepared in triplicates, each one had 30 seedlings. Referring to the previous studies on the early infection process of Fusarium oxysporum on host plants, we collected samples of plant rhizome parts at 2 hpi and 8 hpi [16]. The samples were immediately processed to extract RNA or were frozen with liquid nitrogen and stored at -80 ℃ until extracion [17]. To decrease variability, RNA samples were pooled in groups of three, resulting in two pooled biological replicates per treatment and time point, which were used for RNA-Seq [1].

R-7 and S-29 inoculated with
Foln were tested in the field, and the disease of the plants was observed after maturity, the infected samples were selected for fungal isolation [18]. In order to confirm that symptoms were a result of the fungal infection, the inoculated plants were sampled for fungal isolation from the flax roots. Collected root tip sections (3-5cm) were surface sterilized in 10% sodium hypochlorite for 30 s and then rinsed three times in sterile distilled water. Root sections were further cut into 3～5 mm sections and air-dried on Whatman paper. Four to six of these root sections were transferred to sterile Komada medium [19] and grown for 7 days at 25 ℃ under 16 h dark/8 h light cycles. Purified strains 6 were examined for growth, and colonies were subcultured in PDA for 14 additional days at 25℃ under 16 h dark/8 h light cycles.

Bioinformatics analysis
The total RNAs were extracted using Trizol reagent and transcriptome sequencing was carried using the Illumina HiSeq-(TM) 2000 (Meiji Biopharmaceutical Technology Co., Ltd., Shanghai, China). We filtered the raw sequencing reads to eliminate the low-quality reads and adaptor sequences, leaving only the clean reads. Then, we blasted the clean reads against the whole flax genome (https://www.ncbi.nlm.nih.gov/nuccore/ NC_036356.1) to assemble the matched sequence data into transcripts using Cufflinks software. Cufflinks analyzed all the transcripts obtained during this transcriptome sequencing and their corresponding genes in order to confirm the DEGs. We used Transcripts Per Million (TPM) reads as the index when measuring gene expression levels and we set the parameters as Padjust < 0.05 &|log2FC|≧1. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment via GO-seq (Rlease2.1.2) and KOBAS (v2.0) in order to characterize flax's resistance mechanism.

Validating transcript levels using qRT-PCR
In order to validate the expression profile, we selected eight DEGs that were involved in different metabolic pathways for quantitative reverse transcription PCR (qRT-PCR) The flax actin gene (EU830342) was used as the reference gene for normalization. Genespecific primers were designed by Primer Premier 6.0 and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China).
We took the total RNA isolated from the flax root as a template. The first-strand cDNA

Resistance characteristics and pathogen isolation of flax R-7 and S-29
After inoculation of Foln at four-leaf-stage, the seeding were placed in the experimental environment and field environment for one month Respectively. As shown in Fig.1, that the incidence of R-7 was significantly lower than that of S-29, R-7 seedlings were grown by hydroponic culture in light incubator, inoculation of Foln at 38 days, that the growth state was basically the same as that of control. (Fig. 1-A, Fig. 1  hpi by high-throughput sequencing platform. The original data and the data after quality control were statistically analyzed ( Table 1). The raw reads were 42834168 -66555862 and clean reads were 42118732 -66033802. About 82.48% -95.06% of the sequences were compared with known genes, and the percentage of Q30 base was more than 94%. The above data showed that the sequencing data was true and reliable, indicating that the accuracy was high, and the quality of sequencing data was enough for further bioinformatics analysis.

Genes function annotation
All the genes and transcripts assembled by transcriptome were compared with six databases (NR, Swiss-Prot, Pfam, EggNOG, GO and KEGG), and the functional information of genes and transcripts was obtained comprehensively. A total of 40641 transcripts (93.46%) and 40628 genes (93.46%) were successfully annotated (Fig. 3). Among them, NR and COG functional annotation accounted for 92.44% and 90.34% respectively, Pfam and Swiss functional annotation accounted for 76.6% and 74.87% respectively, and GO and KEGG functional annotation accounted for 51.36% and 40.23% respectively.

DEGs in resistant and susceptible flax cultivars
When compared with the control groups of the R-7 and S-29 collected at 0 hpi, we found significant differences in the number of up and down-regulated DEGs at 0.5 hpi, 2 hpi, and 8 hpi (Fig 4). Additionally, the number of DEGs gradually increased over time.  We observed a sharp increase in the number of DEGs in different cultivars: between 0.5 and 2 hpi of R-7 and between 2 and 8 hpi for S-29.

GO enrichment analysis of DEGs
In order to compare the GO terms of DEGs in resistant and susceptible cultivars of flax at different times, we performed a GO term enrichment analysis (Fig.5). It showed that some of the same genes were induced or inhibited in both flax cultivars at different times of biotic stress. Genes involved in the oxidoreductase activity and oxidation-reduction process were up-regulated expressed in R-7 and S-29 at the same time. The genes of extracellular region and organic phosphate transport were down-regulated expressed in R-7 and S-29 at the same time. Results showed significant differences in gene expression regulation between resistant and susceptible flax cultivars, except for common gene expression, for example, the upregulated expressed genes in R-7 were enriched in terpenoid metabolicess and aldehydelyase actibity, while the up-regulated expressed genes in S-29 mainly include catalytic activity, regulation of cell death, lipid metalic process. Among the abundant down-regulated genes, R-7 mainly had galacturonidase activity, carbohydrate metabolism process and catalytic Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 15 October 2021 activity; S-29 mainly had stem vascular pattern formation and cotyledon vascular pattern formation, phloem pattern formation, cell wall organization or biogenesis and other related information. Among them, the more of cell wall related biological processes .

KEGG enrichment analysis of DEGs
In order to determine the biological functions of DEGs involved in the metabolic pathways of resistant and susceptible flax cultivars under pathogen stress, we performed a KEGG functional enrichment analysis on DEGs. In total, we found significant enrichment in nine KEGG pathways that were involved in both cultivars of flax, namely photosynthesisantenna proteins, zeatin biosynthesis, plant hormone signal transduction, mitogen-activated protein kinase (MAPK) signaling pathway, plant-pathogen interaction, estrogen signaling pathway, phenylpropanoid biosynthesis, starch and sucrose metabolism, and α-linolenic acid metabolism ( Table 2). The DEGs expression of the three pathways was significantly different between R-7 and S-29, including MAPK signaling pathway, plant hormone signal transduction, plant pathogen interaction and phenylpropanoid biosynthesis.
The genes involved in MAPK signaling pathway were mainly encoding serine/threonine protein kinase (IRAK4) family and encoding heat shock 70 kDa protein (HAPS1s). The differential expression level of Lus10026099.g in S-29 was five times more than that in R-7, which was up-regulated expressed at 0.5 hpi and 2 hpi in S-29 and downregulated expressed at 2 hpi and 8 hpi in R-7.

Screening of regulatory genes
The core regulatory genes were screened out through studying interactions of DEGs in enriched KEGGs, including MAPK signaling pathway, plant hormone signal transduction, and plant pathogen interaction. Three genes involved in plant-pathogen interaction pathway had the most correlation (Fig. 6-A), which were Lus10020644.g, Lus10029896.g and Lus10022776.g. Among them, both Lus10020644.g and Lus 10029896.g encoded respiratory burst oxygen homolog protein B, and GO functions were oxidation-reduction 12 process and peroxidase activity. Lus 10022776.g encoded RPM1-interacting protein 4-like, and GO function was cellular process.
The genes involved in the plant hormone signal transduction pathway showed the highest correlation (Fig. 6-B), which were Lus 10002576.g, Lus 10026447.g and Lus10025000.g, respectively. Among them, Lus10002576.g encoded the protein TIFY 9like isoform X1. Both Lus10026447.g and Lus10025000.g encoded protein phosphatase 2C, and GO functions were hydrolase activity. It shown as Fig6-A, B. that Lus10026447.g and Lus10025000.g regulated the genes in MAPK signaling pathway and plant hormone signal transduction.
The analysis of interaction network of MAPK signaling pathway showed that the Lus10025000.g and Lus10027091.g genes had the highest correlation ( Fig. 6-C). In which, Lus10025000.g encoded protein phosphatase 2C, Lus10027091.g encoded mitogen activated protein kinase homolog MMK1,respectively, and the GO function of them were MAP kinase activity, ATP binding and MAPK cascade.

DEG validation by qRT-PCR
Using qRT-PCR, we found that the relative expression levels of eight DEGs at different post inoculation times were consistent with the findings from the RNA-seq results (Fig 7).
As expected, these data confirmed the reliability of the RNA-seq transcriptome analysis.

RNA-sequencing based on second-generation Illumina high-throughput deep
sequencing technology has became an important means for gene mining. Currently RNAseq technique had been widely applied in the study of the interactions between plants and pathogens, it comprehensively, rapidly, and accurately revealed the transcriptome activities of specific species under certain conditions [22] . Identification of key factors in responsed to flax-Foln interaction transcriptome responses will helps to identify and annotate important Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 15 October 2021 genes in plant defense responses. In this study, transcriptome analysis was performed on the resistant and susceptible flax cultivars inoculated with Foln.

Gene expression changes and GO analysis in response to Foln infection
In this study, the number of DEGs in R-7 and S-29 increased significantly at 2 hpi and 8 hpi, and they formed a kind of saltus at 0.5-2 hpi and 2-8 hpi, respectively. This may have been due to the fact that the early stress response of R-7 to Foln infection was more sensitive and active than that in S-29. However, S-29 had more DEGs than that in R-7 at 0.5 hpi. This may had been because the signal elements of S-29 signalled triggered hypersensitive reaction at 0.5 hpi, while R-7's reactions were not triggered until 2 hpi. The number of DEGs in R-7 was significantly more than that in S-29, which was similar to the previous studied on resistant and susceptible cultivars infected by Plasmodiophora brassicae [23]. Our findings indicated that resistant and susceptible flax cultivars had significant differences in their stress response regulation. The up-regulated expressed genes in R-7 also enriched terpenoid biological process.
Terpenoids are a class of metabolites widely existing in plants. Some terpenoids are necessary for plant growth and development, such as gibberellin, abscisic acid and other plant hormones, carotenoids and so on [25]. Terpenoids can rapidly accumulate in the infected parts of plants to prevent the growth and reproduction of pathogens, thus improving the resistance of plants [26]. These substances were mainly involved in plant hormone signal transduction and photosynthesis. At this time, the up-regulated genes in S-29 were enriched in lipid metabolic process and enzyme activities of some organic substances. These biological processes were involved in the metabolism of plant secondary metabolites, which played an important role in plant disease resistance. Phenylpropanoid, flavonoid and stilbene biosynthesis pathways were involved in the formation of secondary resistant metabolites, such as plant antitoxins, lignin and phenols [8].
The down-regulated expressed genes in R-7 and S-29 were enriched in extracellular region, cell wall and organophosphate ester transport. It found that the main genes involved in the biological process of extracellular region were peroxidase and plant pathogenesisrelated protein. Down-regulated expression of these genes reduced the resistance of plants.
A total of seven genes encoded plant pathogenesis-related protein were found in this study.
Among them, Lus10007102.g and Lus10020480.g were up-regulated in R-7 with the extension of inoculation time, while in S-29 was down-regulated. It was judged that the gene was a positive regulatory resistance gene. The gene with higher expression level was Lus10026729.g encoded pectin methylesterase, which was directly involved in the conversion of pentose and glucuronic acid. The GO function of Lus10026729.g was cell wall modification. Pectin methylesterase affected the formation of pectin. Pectin affected the expression of formate dehydrogenase gene after flax was infected with pathogen. The formate dehydrogenase gene was directly involved in the synthesis of cell wall glycan [11] .
ABC transporter G family member 6-like was the main gene involved in the biological process of organophosphate ester transport. Transporters of G subfamily (ABCG) were considered to be an important part of the plant immune system. They participated in the active transmembrane transport of various secondary metabolites and played an important role in plant microbial interaction [27]. This protein family was a kind of transmembrane transporters widely existed in prokaryotes and eukaryotes. Arabidopsis thaliana has about 130 genes encoded the ABC transporter protein, which played a key role in plant hormone transport, lipid metabolism, detoxification of exogenous toxins, plant disease resistance and so on [28]. It found that the down-regulated genes of S-29 were enriched in cell wall and its related biological processes. Cell wall was the first effective barrier for plants to resisted foreign infection. When plants were infected by pathogens, the cell wall will be partially damaged, which induced cells to self-repair the cell wall, accelerated the deposition of nonprotein substances on the cell wall, and made the cell wall thickened and lignified [29], and down-regulated expression of related genes, indicating that the resistance of S-29 was decreased.  [32]. It has been reported that WRKY8 in tobacco stimulates the production of isoprenoid plant antitoxin, and WRKY33 in Arabidopsis thaliana promotes the synthesis of antitoxin [33][34][35]. Transcription factor WRKY33 was highly expressed in R-7, but the downregulated expression trend was observed at 2 hpi and 8 hpi, indicating that R-7 successfully resisted the invasion of pathogenic bacteria immediately after the contact of Foln. Lei et al.

Role of Transcription Factor WRKY in Plant Stress Resistance
found that the sensitivity of GhWRKY22 silenced plants to Verticillium dahliae was increased, and the expression levels of disease resistance marker genes (PR1, PR3, PR4, PR5, PAL and PDF1.2) were also significantly decreased [36]. Eights Transcription factor WRKY genes were involved in the regulation of Verticillum wilt resistance in cotton, and GhWRKY22 positively regulated the resistance of cotton [36]. In this study, it showed that WRKY22 was up-regulated at 2-8 hpi and positively regulated flax resistance in R-7. Zhou et al.analyzed the expression of four genes in Newhall navel orange and four seasons orange induced by ulcer bacteria, indicating that CSWRKY22 was involved in the susceptible response of the host [37] .In this study, S-29 showed a low expression level of WRKY22, but with the extension of Foln inoculation time, the expression tended to be up-regulated, indicating that Foln was involved in the susceptible response of the host, which was consistent with the results of Zhou's study.  Transcription factor MYC played a core regulatory role in JA hormone stimulation and directly regulated downstream response genes. When plants were stimulated by jasmonic acid, the expression of MYC was up-regulated and abscisic acid was also positively regulated. Currently, transcription factor MYC had been found to be involved in JA and ABA hormone signal transduction in Arabidopsis thaliana [39]. Three genes (Lus10000037.g, Lus10035365.g, Lus10035366.g) encoding transcription factor MYC continuously upregulated expressed in R-7, indicated that MYC was actived in hormone signal transduction.

Role of transcription factors MYC and TIFY in plant disease resistance
JA and SA were also related to pathogen resistance, both of which were response of plant responsed to external damage (machinery, herbivores, insects) and pathogen infection to induced the expression of resistance genes. The signal transduction pathway mediated by JA and SA was closely related to plant resistance [40]. Therefore, the effect of Foln infection on endogenous hormones of flax was maybe an important cause of flax fusarium wilt.
In R-7, three genes that encoded transcription factor MYC and three genes encoded TIFY proteins, were up-regulated expressed after 0.5 hpi to increased flax resistance by regulated the synthesis of jasmonic acid. The TIFY family was a new family of plant-specific genes involved in the regulation of a cultivar of plant-specific biological processes, such as the development and response of plant hormones [41]. TIFY10/JAZ was a key regulatory factor of plant jasmonic acid signal transduction that was continuously up-regulated expressed during the induction of draconis sanguis [42]. After inoculation of Foln, both genes encoded TIFY10A and transcription factor MYC were significantly down-regulated expressed in the process of plant hormone signaling, and the synthesis of jasmonic acid was reduced, which reduced plant resistance in S-29. immune system initiated by resistant proteins [43]. The genes RARI and SGTI also played a role in many R-protein-initiated plant disease resistance pathways. SGTI and RAR had been shown to interact with Hsp 90 in tobacco and Arabidopsis thaliana [44].

Data availability
The datasets generated during and/or analysed during the current study are available in the NCBI repository, https://submit.ncbi.nlm.nih.gov/subs/bioproject/SUB8809588/overview Author contributions, Yingdong Chang, Jianming Wang    Table 1 Transcriptome sequencing quality control data statistics Note: (1) Name of Sample; (2) hpi: Hours post-inoculation;(3)replicate: two biological replicates;(4) Raw reads: the total number of items in the original sequencing data (reads, representing a sequencing read); (5) Clean reads: the total number of items in sequencing data after quality control; (6-7) Q20 (%) and Q30 (%) : the quality of sequencing data after quality control is evaluated. Q20 and Q30 refer to the percentage of the total bases whose sequencing quality is above 99% and 99.9%, respectively. Generally, Q20 is above 85% and Q30 is above 80%. (8) Total mapped: number of Clean reads that can be mapped into the genome; (9) Multiple mapped: number of Clean reads with Multiple mapped locations in the reference sequence; (10) Unique mapped: Number of Clean reads with a Unique alignment position in the reference sequence.