Lignin Synthesis Related Genes with Potential Significance in the Response of Upland Cotton to Fusarium Wilt Identified by Transcriptome Profiling

Fusarium wilt, caused by the fungus Fusarium oxysporum Schlecht. f. sp. vasinfectum (Fov) is a destructive soil-borne cotton disease. To profile the genes and pathways responding to Fov infection, we compared transcriptomic responses before and after F. oxysporum inoculation in a highly resistant cotton cultivar, Yumian21, and a highly susceptible cultivar, Jimian11. Although the overall gene expression pattern was downregulated in both cultivars, the global gene expression in the resistant cultivar was stronger than that in the susceptible cultivar. In addition, the expressed genes of two cultivars mostly differed in “cellular process,” “single-organism process,” “metabolic process,” and “response to stimulus” functional groups in the biological process Gene Ontology category: the upregulated differentially expressed genes (DEG) were largely enriched in the resistant cultivar, while the downregulated DEGs were largely enriched in the susceptible cultivar. Phenylpropanoid biosynthesis and phenylalanine metabolism are the key metabolic pathways in cotton in response to Fov. We found that lignin plays a potential role in cotton resistance to Fov. Two coding genes, caffeic acid 3-O-methyltransferase and peroxidase2, as well as the two transcription factors MYB46 and MYB86, are possibly involved in the accumulation and synthesis of lignin. Furthermore, the result showed that the quantification of lignin could be potentially used as a selection tool to identify Fusarium wilt resistant cotton.


Introduction
Cotton (Gossypium hirsutum L.) is an important industrial crop that provides natural fibers and oilseeds (Ullah et al. 2017). Cotton Fusarium wilt (FW), caused by the soil-borne fungus Fusarium oxysporum Schlecht. F. sp. Vasinfectum  (Assigbetse et al. 1994), is a major disease of cotton capable of causing significant economic loss worldwide (Cianchetta and Davis 2015). Fov infects the roots and rapidly colonizes the vascular system of susceptible cotton hosts, causing vascular discoloration, wilting, and the eventual death of the plant (Dowd et al. 2004). FW is difficult to prevent because of long-term persistence of chlamydospores in soil and Fov's resting structures in host xylem vessels (Mei et al. 2014). In China, the spread of FW is of great concern due to the limited knowledge about key aspects of disease epidemiology and the lack of effective management models, including resistant cultivars and soil management approaches.
Plants and their pathogens often evolve through competitive interactions and form complex relationships (Chisholm et al. 2006;Dangl and Jones 2001). The first line of plant defense system is a physical barrier composed of the waxy cuticle and cell wall. If the first line of defense is breached, the plant must resort to a different set of defense mechanisms, such as toxins and enzymes. Plants use their innate immune system to fight against pathogens (Malinovsky et al. 2014). With the primary immune system, plants recognize microbe-associated molecular patterns of potential pathogens through pattern recognition receptors that mediate a basal defense response. Plant pathogens suppress this basal defense response by means of effectors that enable them to cause disease. With the secondary immune system, plants have the ability to recognize effector-induced perturbations of host targets through resistance proteins that mediate a strong local defense response that stops pathogen growth (Jones and Dangl 2006). Currently, the possible resistance mechanisms of cotton's defense against FW are poorly understood.
In cotton plants, a series of metabolic changes can be induced during the initial response to pathogen invasion. For example, the affected cells produce secretions which can bury or encase pathogens; they may also block vessels by increasing metabolic activities, thereby limiting the spread of pathogens (Shi et al. 1991;Shi et al. 1992;Shi et al. 1993). Additionally, phenolic compounds such as gossypol and its derivatives may be induced from infected cells, further inhibiting the spread of F. oxysporum (Hall et al. 2011;Mellon et al. 2014;Zhang et al. 1993). Furthermore, many kinds of transcription factors, protein kinases, and phytohormones are involved during the response process of cotton plants to Fov: for example, the cotton MAPK kinase GhMPK20 negatively regulates resistance to F. oxysporum, and the expression of GhMKK6 contributes to the immune response in cotton Wang et al. 2017). Notably, both salicylic acid and jasmonic acid can contribute to the over-expression of pathogenesis-related genes in cotton plants (Akhunov et al. 2008).
Lignin plays important roles in the growth and development of plants, defending against a variety of biotic and abiotic stresses. It acts as a physical and chemical barrier to limit pathogen colonization and restrict pathogen growth in a wide range of plant species (Bonello and Blodgett 2003;Hückelhoven 2007;Zhang et al. 2017). Lignin content has been used as a biochemical marker of an activated immune response in a multitude of species (Adams-Phillips et al. 2010;Kishi-Kaboshi et al. 2010): lignin content in disease-susceptible tomato varieties is significantly lower than that in disease-resistant cultivars (Mandal et al. 2013); the composition and accumulation of lignin strengthens the cell walls and enhances plant resistance to fungi in pepper plants (Novo et al. 2017); and lignin content increases in Chinese cabbage when it is infected by Erwinia carotovora (Zhang et al. 2007). The lignin biosynthesis genes play significant roles in resisting the invasion of pathogens through the salicylic acid defense pathway in Arabidopsis thaliana (Tronchet et al. 2010). When rice suffers pathogen invasion, the activity of cinnamoyl-CoA increases to positively regulate lignin synthesis (Kawasaki et al. 2006) and OsAAE3 negatively regulates lignin synthesis (Liu et al. 2017). In maize, genes HCT1806, HCT4918 and ZmCCoAOMT2 are involved in the resistance of various pathogens by participating in the biosynthesis of lignin Yang et al. 2017). Lignin plays a critical role in cotton resistance to Verticillium dahliae: the resistant cotton varieties show a high level of lignin deposition and ligninlike phenolic polymers, (Smit and Dubery 1997) and the lignin content is positively correlated with resistance to V. dahliae (Xu et al. 2011). The GhDIR1, GbERF1-like, GhUMC1, GhWAT and GhLAC15 genes improve cotton resistance to V. dahliae via activating or strengthening lignin synthesis (Guo et al. 2016;Shi et al. 2012;Tang et al. 2019b;Zhang et al. 2019;Zhu et al. 2018). We hypothesize that lignin synthesis-related genes contributed to the cotton resistance to Fov through positively regulating lignin synthesis and enhancing the lignin content.
In the present study, one highly resistant cotton cultivar Yumian21 and one highly susceptible cultivar Jimian11 were chosen to analyze the molecular mechanism underlying cotton response to Fov. The gene expression patterns and functions of differentially expressed genes (DEGs) were analyzed based on the RNA-seq data. We found that the deposition and accumulation of lignin play a potential role during the response of cotton to Fov invasion. Two coding genes: caffeic acid 3-Omethyltransferase and peroxidase2, and two transcription factors: MYB46 and MYB86 are possibly involved in the accumulation and synthesis of lignin. Lignin quantification can be used as a selection tool to identify FW resistant cotton. The present study can also be served as a resource for future research addressing the resistance mechanism of cotton to Fov.

Resistant Phenotype of Two Cotton Cultivar after Inoculation with Fov
When resistant and susceptible cultivars were inoculated with F. oxysporum, the phenotype of disease began to appear in 16 days. The resistant cultivar survived well and only exhibited a slightly shorter plant size ( Fig. 1A and C). On the contrary, the susceptible cultivar showed several symptoms including dwarfing, leaf wilting, defoliation, vascular discoloration and necrosis, and death of plants ( Fig. 1B and D). The histochemical analysis of lignin in stem cross-sections showed no significant differences of lignin staining between the resistant cultivar inoculated with distilled water (Fig. S1 A), and the susceptible cultivar inoculated with distilled water or Fov ( Fig. S1 B and D). However, the lignin staining of the resistant cultivar inoculated with Fov was significantly darker than the susceptible cultivar inoculated with Fov ( Fig. S1 C and D).

The Quality of Sequencing and Validation by qRT-PCR
The 12 libraries were sequenced using the Illumina HiSeq X Ten sequencing platform. As shown in Table S1 (Table S2).
We compared RNA-Seq data with real-time RT-PCR data for 18 randomly selected DEGs that showed up-regulation or down-regulation in response to Fov infection according to RNA-Seq (Excel S1). The results of the qRT-PCR data for these 18 genes were consistent with those obtained via RNAseq analysis, which indicated that the changes in expression detected by RNA-seq were accurate (Fig. S2).

Patterns of Gene Upland Cotton in Response to Infection by Fov
In our experiment, a false discovery rate (FDR) < 0.05 and fold change (FC) =1.5 were used as the thresholds for assessing significant differences in gene expression. A total of 1055 DEGs were detected in the resistant cultivar, including 477 (45%) upregulated genes and 578 (55%) downregulated genes. A total of 995 DEGs were identified in the susceptible cultivar, including 369 (37%) upregulated genes and  Fig. 2A). Genes were mostly downregulated rather than upregulated when both cultivars were inoculated with Fov. In addition, the number of downregulated genes (578) in the resistant cultivar was less than that (626) in the susceptible cultivar. A total of 146 genes were either up-regulated or down-regulated in resistant and susceptible cultivar: 27 DEGs were upregulated in both cultivars; 103 DEGs were downregulated in both cultivars; ten DEGs were upregulated in the resistant cultivar but downregulated in the susceptible cultivar; and six DEGs were upregulated in the susceptible cultivar but downregulated in the resistant cultivar (Fig. 2B). The concrete DEGs of the resistant cultivar and the susceptible cultivar after treatment with Fov are listed in Excel S2 and S3, respectively.

Functional Classification of DEGs
A total of 1055 DEGs were divided into three main GO categories: biological process, cellular components, and molecular functions for the resistant cultivar. Most of the GO terms were assigned to biological process (56.13%), followed by cellular component (31.30%) and molecular functions (12.57%). Of 995 DEGs in the susceptible cultivar, 54.82%, 32.36% and 12.82% were assigned to biological process, cellular component and molecular function GO terms, respectively (Excel S4 and S5). The greatest proportion of the GO categories consisted of biological process in both resistant and susceptible cultivars.
In the susceptible cultivar, the upregulated DEGs were distributed into 18 groups of the biological process category. Of 1932 upregulated genes, 264 genes (13.66%) were involved in "metabolic process," followed by 262 genes (13.56%) involved in "cellular process." The downregulated DEGs were distributed into 16 groups of the biological process category. Of 3667 downregulated genes, 471 (12.84%) were involved in "cellular process", followed by 469 genes (12.78%) involved in "single-organism process".
In "biological process" GO categories, the number of upregulated genes in the resistant cultivar was significantly higher than that in the susceptible cultivar (Fig. 3). The upregulated DEGs were mostly enriched in four functional groups: the "cellular process" group (349 DEGs in resistant, 262 in susceptible), "single-organism process" group (343 in resistant, 257 in susceptible), "metabolic process" group (329 in resistant, 264 in susceptible), and "response to stimulus" group (307 in resistant, 223 in susceptible). In contrast, the number of downregulated genes in the susceptible cultivar was significantly higher than that in the resistant cultivar (Fig. 3). The downregulated DEGs were mostly enriched in "cellular process" (399 in resistant, 471 in susceptible), "single-organism process" (412 in resistant, 469 in susceptible), "metabolic process" (389 in resistant, 446 in susceptible), and "response to stimulus" (358 in resistant, 408 in susceptible). The resistant cultivar and the susceptible cultivar mostly differed in "cellular process," "single-organism process," "metabolic process," and "response to stimulus" functional groups: the upregulated DEGs were largely enriched in the resistant

KEGG Enrichment Analysis of Cotton Infection of Fov
For the resistant cultivar, 207 DEGs (87 upregulated vs 120 downregulated) were identified in the 74 KEGG pathway by BLAST analysis against the KEGG database. The DEGs were significantly enriched in the first 20 pathways with the lowest p value when the resistant cultivar was inoculated with Fov ( Fig. 4, the specific data was listed in Excel S6). The top three highly enriched pathways with the largest enrichment factors, the largest number of DEGs, and the lowest p value were phenylpropanoid biosynthesis (ko00940), phenylalanine metabolism (ko00360), and pentose and glucuronate interconversions (ko00040) in the resistant cultivar. For the susceptible cultivar, 210 DEGs (88 upregulated vs 122 downregulated) were activated in the 84 KEGG pathway. The 20 with the most reliable pathway functional enrichment of the susceptible cultivar after inoculation with Fov is shown in Fig. S3 (Specific data in Excel S7). The top three highly enriched pathways were phenylpropanoid biosynthesis (ko00940), phenylalanine metabolism (ko00360), and flavonoid biosynthesis (ko00941); the phenylpropanoid biosynthesis (ko00940) and phenylalanine metabolism (ko00360) pathways were the most noteworthy KEGG pathways in both cultivars. These results showed phenylpropanoid biosynthesis and phenylalanine metabolism are the key metabolic pathways in cotton in response to Fov.

Transcription Factors (TFs) Involved in Cotton Infection of Fov
We investigated the distribution of transcription factors when two cultivars of cotton were infected with Fov. We found that the DEGs have transcription factors belong to AP2/ERF, MYB, NAC, bHLH, bZIP and WRKY families ( Table 1). The numbers of upregulated TFs (30) exceeded downregulated TFs (15) in the resistant cultivar. However, the number of downregulated TFs (28) exceeded upregulated TFs (11) in the susceptible cultivar. We found most TFs belong to AP2/ERF and MYB (18 and 24, respectively). In the AP2/ERF family, twelve genes were upregulated, while one was downregulated in the resistant cultivar. However, only two genes were upregulated, while three were downregulated in the susceptible cultivar. In the MYB family, twelve genes were upregulated and four were downregulated in the resistant cultivar. Two genes were upregulated and twelve were downregulated in the susceptible cultivar. In the resistant cultivar, the total number of upregulated genes (24) of these two families was significantly higher than total downregulated genes (5). In the susceptible cultivar, the total number of upregulated genes (4) of these two families was significantly lower than total downregulated genes (15).

Genes Related to Disease Resistance
We found that two genes: caffeic acid 3-O-methyltransferase ( e v m . T U . G h _ A 0 9 G 1 0 5 3 ) a n d P e r o x i d a s e 2 (evm.TU.Gh_D09G1208) were likely involved in resistance to Fov in cotton. First, we screened 1328 upregulated DEGs in cellular process, metabolic process, response to stimulus and single-organism process functional groups of classification in the resistant cultivar. Second, we screened 7 upregulated DEGs in phenylpropanoid biosynthesis (ko00940) and phenylalanine metabolism (ko00360) based on the analysis of the KEGG enrichment pathway. The result demonstrated that c a f f e i c a c i d 3 -O -m e t h y l t r a n s f e r a s e ( e v m . T U . G h _ A 0 9 G 1 0 5 3 ) a n d P e r o x i d a s e 2 (evm.TU.Gh_D09G1208) were significantly upregulated in the resistant cultivar, but not in the susceptible cultivar (Table 2). We found that two transcription factors MYB46 ( e v m . T U . G h _ D 1 3 G 2 2 6 1 ) a n d M Y B 8 6 (evm.TU.Gh_D08G1266) were involved in resistance. Among 1328 upregulated DEGs in cellular process, metabolic process, response to stimulus, and single-organism process functional groups, 24 TFs belong to MYB family in the resistant cultivar. MYB46 (evm.TU.Gh_D13G2261) and MYB86 (evm.TU.Gh_D08G1266) were uniquely upregulated in the resistant cultivar, but not in the susceptible cultivar (Table 2). When the gene expression levels were normalized to the susceptible cultivar at 0 h, a relative high expression level of caffeic acid 3-O-methyltransferase, Peroxidase 2, MYB46 and MYB86 genes were observed in the resistant cultivar treated with Fov (Table 2). No major shifts in expression levels of these genes were observed in the susceptible cultivar (Table 2) Interestingly, these two genes and two transcription factors are involved in the process of lignin synthesis, the activation of peroxidase activity, and the positive regulation of secondary cell wall biogenesis in GO second level annotation. This indicates that lignin synthesis may play potential roles in defending against Fov.

POD Activity and Lignin Content of Cotton Inoculated by Fov
There was no significant difference of lignin concentration in the resistant cultivar treated with Fov (72.8μg/L) or distilled Fig. 4 Pathway functional enrichment of DEGs in the Yumian21 resistant cultivar inoculated with Fusarium oxysporum. The x-axis represents the enrichment factor (rich factor) which is the ratio of the foreground value (the number of DEGs) and the background value (total gene amount). The y-axis shows the pathway names. A larger value of the rich factor indicates a higher enrichment value. The color indicates the p value, A lower p value refers to a more significant enrichment. Point size indicates DEG number and larger dots refer to higher numbers of DEGs Although the lignin concentration increased with the relative expression level of four genes, no strong correlation (P < 0.05) was found between lignin concentration with caffeic acid 3-O-methyltransferase (r = 0.514), Peroxidase 2 (r = 0.526), MYB46 (r = 0.221) and MYB86 (r = 0.609) genes. The lignin concentration of distilled water treatments was 70.3 μg/L, 75.5 μg/L, 76.9 μg/L and 82.7 μg/L at 12 h, 24 h, 48 h and 96 h, respectively. The concentration of lignin in the Fov treatment was significantly higher than that in the water treatments (Fig. 5A). In the susceptible cultivar, the lignin concentration of distilled water treatments was 68 μg/L, 70.9 μg/L, 77.8 μg/L, 78.9 μg/L and 83.9 μg/L at 0 h, 12 h, 24 h, 48 h and 96 h, respectively.  5B). No significant change was found in the susceptible cultivar treated with distilled water or Fov. As one of the key enzymes in the antioxidant system in plants, peroxidase (POD) plays an important role in the immune defense mechanism, and its expression is induced by fungi. At the same time, it is also an important control enzyme in the lignin synthesis pathway. In this research, there was no significant difference of POD activity in the resistant cultivar treated with Fov (39 mU/L) or distilled water (36.1 mU/L) before Fov inoculation. The POD activity of Fov treatments increased to 59.8 mU/L, 63.1 mU/L, 49.6 mU/L and 43.2 mU/ L at 12 h, 24 h, 48 h and 96 h, respectively. In the resistant cultivar inoculated with Fov, the POD activity showed a statistically significant change over 12 h, 24 h, 48 h and 96 h timepoints (P < 0.01) and it reached the highest level at 24 h. The relative expression level of Peroxidase 2 gene was highly correlated with the POD activity (r = 0.969, P < 0.05) The POD activity of distilled water treatments was 38.5 mU/L, 46.1 mU/L, 42.4 mU/L and 35.8 mU/L at 12 h, 24 h, 48 h and 96 h, respectively. The activity of POD in the Fov treatment was significantly higher than that in the distilled water treatments (Fig. 5C). In the susceptible cultivar, The POD activity of distilled water treatments was 34.4 mU/L, 37.7 mU/L, 43 mU/L, 42.2 mU/L and 38.5 mU/L at 0 h, 12 h, 24 h, 48 h and 96 h, respectively. The POD activity of Fov treatments was 36.3 mU/L, 39.5 mU/L, 47.5 mU/L, 44.6 mU/L and 40 mU/L at 0 h, 12 h, 24 h, 48 h and 96 h, respectively (Fig. 5D). No significant change was observed in the susceptible cultivar treated with distilled water or Fov.

Discussion
We found 1055 DEGs in the resistant cultivar and 995 DEGs in the susceptible cultivar. The global gene expression in the resistant cultivar is stronger than that in the susceptible cultivar. Interestingly, more downregulated genes were identified than upregulated genes in both cultivars after inoculation by Fov. Xu et al. found the same pattern with cotton inoculated with Verticillium wilt (Xu et al. 2011). However, the proportion of downregulated genes in the resistant cultivar (55%) was less than that in the susceptible cultivar (63%). Although the developmental processes of two cultivars were all affected by Fov invasion, the resistant cultivar was less affected than the susceptible cultivar.
GO classification showed that cellular process, singleorganism process, metabolic process and response to stimulus were the main enriched GO terms of the biological process group and all of them were involved in response to Fov infection. By further analysis of the KEGG pathway, we found that in both cultivars the most significantly enriched KEGG pathways were "phenylpropanoid biosynthesis," followed by "phenylalanine metabolism." This indicates that these pathways are the key metabolic pathways in cotton in response to FW. Phenylpropanoids can function as inducible antimicrobial compounds as well as signaling molecules in plantpathogen responses (Dixon 2001;Naoumkina et al. 2010). Phenylpropanoid metabolism has been well documented as the most important plant metabolic pathway during plant defense against biotic stress (Cass et al. 2015;La Camera et al. 2004). The phenylpropanoid pathway contains multiple branches, and lignin biosynthesis is a downstream branch of this pathway; specifically, the pathway synthesizes monolignols, which are substrates of lignin polymerization (Boerjan et al. 2003;Mottiar et al. 2016). In resistant cotton, genes involved in the phenylpropanoid pathway and lignin accumulation were significantly induced after V. dahliae infection (Xu et al. 2011). Overexpression of GhLac1 and GhLac15, genes related to phenylpropanoid pathway and lignin biosynthesis, was found to enhance V. dahliae resistance Zhang et al. 2018;Zhang et al. 2019). Consistently, we found that genes in the "phenylpropanoid biosynthesis" pathway can promote lignin biosynthesis, therefore enhancing cotton resistance to Fov.
Plant defense-associated genes are normally regulated by transcription factors. TFs can coordinate and control the activity of multiple stress response genes and establish a complex regulatory network to regulate physiological and metabolic responses to cope with disease stress (Buscaill and Rivas 2014;Nakashima et al. 2009). Our results indicate that six TF families (AP2/ERF, MYB, NAC, bHLH, bZIP and WRKY) are activated in cotton plants following Fov infection. In A. thaliana, the majority of F. oxysporum-responsive TFs also belong to the ERF, MYB, NAC, bHLH, and WRKY families ). Studies showed that MYB genes are involved in a plant's response to stress, such as drought and cold and pathogen (Erpen et al. 2018;Liu et al. 2016). When cotton plants were infected by Fov, most of the MYB genes were upregulated in the resistant cultivar but downregulated in the susceptible cultivar. MYB families may play a potential role in the regulation of defense gene expression when upland cotton is attacked by Fov.
Lignin may play a key role in cotton resistance to Fov. Physiological and biochemical studies have shown that the Fusarium-resistant watermelon varieties have strong cell structures, such as thickened xylem, to prevent entry of the pathogen (Tang et al. 2019a). As the major end product of the phenylpropanoid pathway, lignin is essential in the formation of cell walls, which forms a physical barrier that increase the resistance to pathogen infection. Caffeic acid Omethyltransferase is a methylase in the lignin specific pathway. It catalyzes the multi-step methylation of hydroxylated monomer lignin precursor and plays a major role in lignin biosynthesis (O'Malley et al. 1993). Two genes, caffeic acid 3-O-methyltransferase (evm.TU.Gh_A09G1053) and Peroxidase 2 (evm.TU.Gh_D09G1208), were found upregulated in the resistant cotton cultivar but not in the susceptible cultivar. Strong downregulation of caffeic acid 3-Omethyltransferase decrease lignin content in Alfalfa (Guo et al. 2001), maize (Pichon et al. 2006) and Leucaena (Rastogi and Dwivedi 2006). Peroxidase helps polymerize monolignols into lignin and enhances cell wall during pathogen attack and damage (Marjamaa et al. 2009). Gayoso et al. found that POD activity was upregulated in a resistant tomato cultivar infected by Verticillium wilt (Gayoso et al. 2010). As a key gene in lignin synthesis, the peroxidase2 (evm.TU.Gh_D09G1208) coding gene was activated and upregulated in resistant cotton during Fov infection. Transcription factors can also regulate lignin synthesis, for example: the R2R3-MYB transcription factor regulates lignin biosynthesis through binding to the AC element (Rogers et al. 2005), which exists in the promoter region of gene in most phenylpropanoid pathways (Hatton et al. 1995).We found that transcription factors MYB46 (evm.TU.Gh_D13G2261) and MYB86 (evm.TU.Gh_D08G1266) were upregulated in the resistant cultivar, but not in the susceptible cultivar. Among transcription factors regulating phenylpropanoid and lignin biosynthesis, MYB46 and MYB86 are transcriptional activators of the lignin biosynthetic pathway during secondary cell wall formation in Arabidopsis (Taylor-Teeples et al. 2015;Zhong et al. 2007). In this study, DEGs analysis showed that two genes related to lignin biosynthesis and two transcription factors were upregulated in the resistant cotton cultivar. These observations were also confirmed by qRT-PCR analysis. However, the defense response of plant to biotic stress is a systematic network, which involves many pathways and genes (Dangl and Jones 2001). Lignin synthesis is related to many genes controlling a few pathways, not only two genes and two TF. We believe other pathways and genes may also play important roles in the cotton immune defense mechanism to Fov invasion.
Based on DEG analysis, GO classification, KEGG enrichment, and TFs analysis, we concluded that lignin is possibly involved in upland cotton's response to FW. To further confirm this, we measured POD activity and lignin content before and after F. oxysporum inoculation in the resistant cultivar and susceptible cultivars. Results showed that both POD activity and lignin content were significantly induced in the resistant cotton.
Breeding and utilizing FW-resistant cotton cultivars has proven to be the most cost-effective control method for this disease. Lignin related genes could be good targets for cotton disease resistance breeding. Lignin content of crops could also be modified by transgenic approaches: the manipulation of lignin biosynthesis at the regulatory level, controlling monolignol biosynthetic enzymes, or the modification of lignin polymer structure (Frei 2013). Liu et al. suggest several candidate genes for the genetic modification of lignin towards breeding rice with high lodging resistance . The over-expression of caffeic acid 3-O-methyltransferase ( e v m . T U . G h _ A 0 9 G 1 0 5 3 ) a n d p e r o x i d a s e 2 (evm.TU.Gh_D09G1208) genes would be able to facilitate the breeding of FW resistant cotton. Additionally, our data suggests the potential use of lignin quantification as a selection tool to identify FW resistant cotton. Lignin content has been used as a biochemical indicator to measure immune activation response (Adams-Phillips et al. 2010;Kishi-Kaboshi et al. 2010). High lignin concentration has been found associated with disease resistance in tomato (Mandal et al. 2013), rice ) and soybean (Peltier et al. 2009). The near infrared spectroscopy method is a promising tool for lowcost, and high-throughput lignin analysis of a large number of samples.

Plant Materials and Inoculation
The Jimian11 susceptible cultivar has been set as the susceptible standard control in China's national cotton disease and pest resistance evaluation programs since 2009 (GB/22101.4-2009). The Yumian21 cultivar is currently served as a resistant control to FOV in national cotton regional trials in China (Mei et al. 2014). The cultivars were delinted with sulfuric acid and sterilized with 0.1% HgCl for 10-15 min, and grown in sterilized soil (vermiculite: nutritive soil = 1:1) at 25°C temperature, 70% relative humidity with a photoperiod 16 h light and 8 h dark. In this study, "R" indicates resistant cultivar Yumian21, "S" indicates susceptible cultivar Jimian11, "CK" indicates the control, and "I" indicates inoculation.
The highly aggressive strain of the defoliating fungus F. oxysporum was incubated on potato-dextrose agar (PDA) media for two weeks and then cut into small pieces (1-2 cm) before transference into 250 mL conical flasks containing 110 mL of liquid Czapek medium. Each flask containing three pieces of PDA medium with mycelia was shaken at 120 rpm at 25°C for 3 d on a rotary shaker. The suspension liquid was adjusted to 2 × 10 7 spores per mL with distilled water.
Each seedling was inoculated with a 10 mL F. oxysporum spore suspension of 2 × 10 7 spores per mL by watering roots at the two-true-leaf growth stage (Yao et al. 2019). Control plants were inoculated with sterile distilled water and treated in the same way. The plant roots were collected at 0, 12, 24, 48, and 96 h post inoculation and frozen immediately in liquid nitrogen. Samples were pooled together as a mixed infection RNA library. Three replicates were used for each treatment, and 15 plants were used in each replication.

RNA Sequencing and Data Processing
Total RNA was extracted with an RNA prep Pure Plant Kit (TIANGEN). The RNA quality was assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) and Nano 6000 Assay Kit on Agilent Bioanalyzer 2100 System (Agilent Technologies, CA, USA), respectively.
RNA was fragmented randomly and rRNA was removed from the samples by using a Ribo-Zero rRNA removal kit (Epicentre, Madison, WI, USA). The cDNA was synthesized by using a random hexamer primer, and purified with AMPure XP beads (Beckman Coulter, Beverly, USA). The purified double-strand cDNA was end-repaired, poly-A tailadded, and ligated to the Illumina adapter. Clustering of the index-coded samples was performed on an acBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq X Ten system (Illumina Inc., USA).
The sequence alignment and subsequent analysis of the high-quality data was conducted using the designated genome Gossypium_hirsutum_v1.1 as a reference (http://mascotton. njau.edu.cn/info/1054/1118.htm) via HISAT (Li et al. 2014), while StringTie (Kim et al. 2015) was used to assemble the reads that had been mapped to the designated genome. Based on the selected reference genome sequences, the mapped reads were spliced by StringTie and compared with the original annotation information of the genome to find the previously unannotated transcription regions. In order to obtain annotation information using BLAST software, the mapped reads were also compared to the following databases: the Nonredundant (Nr) (Sayers et al. 2019), Swiss-Prot (Apweiler et al. 2004), Gene Ontology(GO) (Ashburner et al. 2000), Cluster of Orthologous Groups of proteins (COG) (Tatusov et al. 2000), Evolutionary Genealogy of Genes: Nonsupervised Orthologous Groups(EggNOG), Eukaryotic Orthologous Groups (KOG) (Koonin et al. 2004), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2004). The FPKM method was used to indicate the transcriptional or gene expression level (Trapnell et al. 2010). The FPKM value was calculated as follows: FPKM = cDNA fragments/ (mapped fragments (millions) X transcript length (kb), with cDNA fragments being the number of reads that aligned to a specific unigene, mapped fragments being the total number of reads that aligned to all unigenes, and transcript length being the length of the unigene.

Identification of Differentially Expressed Genes (DEG) and Data Validation
The DEGs between the inoculated and control samples were identified by edgeR (Robinson et al. 2009). The FC value indicates the ratio of expression between two samples (groups). The false discovery rate (FDR) control method was applied in multiple hypothesis testing to correct the P value. An FDR < 0.05 and FC =1.5 were set as the thresholds for assessing the significance of the difference in gene expression (Reiner et al. 2003).
The level of gene expression obtained by RNA-Seq was verified using quantitative real-time PCR (qRT-PCR). The primer sequences were designed using Primer3 (http:// bioinfo.ut.ee/primer3/) and synthesized by the Beijing Genomics Institute (BGI, Beijing, China). The qRT-PCR was carried out with three technical replicates using NovoStart SYBR qPCR Super Mix Plus Kit (Novoprotein) on an Eppendorf Mastercycler ep realplex machine (Germany) according to the manufacturer's protocol. PCR cycles were as follows: one cycle of 1 min at 95°C, followed by 40 cycles at 95°C for 20s,60°C for 20s and 95°C for 30s. Following amplification, all products were subjected to melt curve analysis. The cotton UBQ7 gene (GenBank: DQ116441) was used as the reference gene to normalize the total amount of 500 ng RNA in each reaction.

Measurement of Peroxidase and Lignin
The peroxidase activity (POD) and lignin content were measured by using an enzyme-linked immunosorbent assay (ELISA) kit (Shanghai Bio-Tech Company, Ltd) according to the manufacturer's instructions (Kapat and Dey 2000;Lequin 2005). Briefly, blank controls, standard wells, and sample wells were set up. A 50 μl of the ELISA-coated reference standard 40 μl of sample replacement solution and 10 μl of the test sample were used. The solutions were gently mixed and incubated at 37°Cfor 45 min. Each well was washed five times using washing buffer. Samples were incubated with 50 μL Biotinylated anti-IgG at 37°C for 30 min; then incubated with 50 μL streptavidin-HRP at 37°C for 15 min. A total of 50uL chromogen solution A and B solutions were added each well and incubate samples at 37°Cfor 15 min. The reactions were stopped by adding 50 μL stop solution. The optical density (OD) values were measured at 450 nm for the standard curve. The data were collected from three replicated plants at 0 h, 12 h, 24 h, 48 h and 96 h after inoculation. The analysis of variance (ANOVA), and all significant differences were examined according to Tukey test by DPS 6.05 software at p < 0.05 (Zhejiang University, Hangzhou, China) (Tang and Zhang 2013).

Histochemical Test
Hand-cut cross-sections were made from the base of the stem of both the Fov and distilled water treated cotton plants at 16 d after treatment. Lignin histochemistry was examined using safranin O-fast green staining method (Jia et al. 2015). The materials were fixed in 50% FAA for 24 h then dehydrated with ethanol, and embedded in paraffin. The 4 μm thickness cut sections were stained with 1% safranin O in distilled water (pH 6.7) for 10 min then washed and counterstained with a 0.1% solution of fast green in water for 5 min. Tissue sections with lignin were stained and observed under a light microscope (LEICA CTR6000 Germany). The presence of lignin was stained with red.