Transcriptome and miRNA analyses of the response to Corynespora cassiicola in cucumber

Cucumber (Cucumis sativus L.) target leaf spot (TLS), which is caused by the fungus Corynespora cassiicola (C. cassiicola), seriously endangers the production of cucumber. In this assay, we performed comprehensive sequencing of the transcriptome and microRNAs (miRNAs) of a resistant cucumber (Jinyou 38) during C. cassiicola inoculation using the Illumina NextSeq 500 platform. The possible genes related to the response to C. cassiicola were associated with plant hormones, transcription factors, primary metabolism, Ca2+ signaling pathways, secondary metabolism and defense genes. In total, 150 target genes of these differentially expressed miRNAs were predicted by the bioinformatic analysis. By analyzing the function of the target genes, several candidate miRNAs that may be related to the response to C. cassiicola stress were selected. We also predicted 7 novel miRNAs and predicted their target genes. Moreover, the expression patterns of the candidate genes and miRNAs were tested by quantitative real-time RT-PCR. According to the analysis, genes and miRNAs associated with secondary metabolism, particularly the phenylpropanoid biosynthesis pathway, may play a major role in the resistance to C. cassiicola stress in cucumber. These results offer a foundation for future studies exploring the mechanism and key genes of resistance to cucumber TLS.

mRNA, and miRNAs are completely complementary for cleavage 13 . Second, the combination of plant miRNAs and target mRNA is not fully matched; thus, mRNA translation is inhibited but does not affect transcription 14 .
The target genes of miRNAs in plants are mostly transcription factors involved in plant growth and the response to environmental stresses. Therefore, miRNAs play a vital role in plant development and the adaptation to adverse situations 15,16 . Certain plant miRNAs can down-regulate the expression of target genes to cope with the demands of growth and environmental stress". Down-regulation is a very complex process, and a feedback inhibition regulation pathway has been proposed to exist among the miRNAs and target genes 17,18 .
Many miRNAs in plants can be induced following pathogen infections; these miRNAs can interfere with the expression of genes by interacting with target genes 19 . In Arabidopsis thaliana, a significant difference was observed in the expression levels of miR156, miR159, miR166, miR825, miR852 and miR843 after Pseudomonas syringae infection 20 . In Verticillium dahliae-inoculated cotton roots, 65 miRNA families exhibited differences in expression 21 . In Chinese wild grape, researchers have identified miRNAs that may be involved in powdery mildew resistance 22 . Similarly, under powdery mildew stress, 79 miRNAs were found in wheat leaves 23 . In cucumbers, some researchers have identified a few miRNAs associated with the development and response to stress [24][25][26][27][28] . However, there are no reports of miRNAs associated with the response to C. cassiicola.
In this study, high-throughput sequencing was performed to investigate the changes in the transcriptome and miRNAs in cucumber infected with C. cassiicola. This assay screened for differentially expressed genes (DEGs) and found miRNAs involved in the response to disease stress in cucumber, providing a new theoretical basis for studies investigating disease resistance mechanisms and key resistance genes in TLS.

Results
Pathogen invasion and plant response. After inoculation with C. cassiicola, resistant (Jinyou 38) and susceptible (Ludixianfeng) cucumber varieties were investigated to detect changes in phenotypes and lesions (Fig. 1). The symptoms of the resistant variety were mild (Fig. 1A). Compared with the resistant variety, the susceptible variety showed more severe symptoms (Fig. 1B). Because the ITS fungal sequence is specific, we used PCR to quickly detect pathogen invasion; the amplified fragment was 291 bp. Approximately 6 h post-inoculation (hpi), the ITS sequence of the pathogen was amplified in the both varieties by PCR (Fig. 2). As shown in Fig. 3, after inoculation with C. cassiicola, the red complex began to appear in both varieties of cucumber leaves at 12 hpi, which indicated that lignin was deposited, and the lignin accumulation increased over time. At 12-24 hpi, staining of the resistant variety was stronger than that of the susceptible variety. This finding indicated that the rate of lignin production in the resistant variety was faster under C. cassiicola stress than that of the susceptible variety. As shown in Fig. 4, brown spots began to appear in the resistant and susceptible cucumber varieties at 6 hpi, which indicated that H 2 O 2 began to accumulate and gradually increase in the leaves. At 24 hpi, a large area in the leaves was stained, showing that the leaves accumulated a higher level of H 2 O 2 . However, at 48 hpi, the accumulation of H 2 O 2 showed a decreasing trend. During the early stage of the pathogen treatment (6-12 hpi), the staining in the leaves of the resistant variety was deeper, indicating that the outbreak of H 2 O 2 in the resistant variety was stronger than that in the susceptible variety. However, at 24 hpi, compared with the resistant cultivar, the susceptible cultivar showed leaves of a deeper color and greater H 2 O 2 accumulation. In the NBT staining experiment, blue spots began to appear in both cucumber varieties at 6 hpi, and large areas of blue began to appear in both varieties at 24 hpi; however, at 48 hpi, the O 2 .− accumulation began to gradually decrease. Similarly, during the early stages of infection, the accumulation rate of O 2 .− was faster and the accumulation was larger in the resistant variety (Fig. 5) Table S1). All quality reads were assembled and annotated using BLAST against the cucumber genome database. eggNOG analysis of functional categories. The eggNOG database can cluster eukaryotic orthologous proteins. The annotated genes were searched in the eggNOG database to determine their functional classification. In total, 9527 genes were classified into 26 categories according to eggNOG. The largest group was "signal transduction mechanisms" (1269, 13.32%), followed by "posttranslational modification, protein turnover, chaperones" (988, 10.37%), "general function prediction only" (839, 8.81%), "carbohydrate transport and metabolism" (605, 6.35%) and "secondary metabolite biosynthesis, transport, and catabolism (569, 5.97%)" (Fig. 6).
Analysis of DEGs. In total, 21 503 genes were reviewed, and pairwise comparisons were performed between each time point (0 hpi vs 6 hpi and 6 hpi vs 24 hpi). Genes with at least a twofold change were considered DEGs (P-value < 0.05).
To further explore the functions of the DEGs, we performed a Gene Ontology (GO) enrichment analysis (Fig. 7). The GO terms of the DEGs were significantly enriched in "metabolic process", "response to stress", "transport", "extracellular region" and "molecular function" (Table 1).
We also performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the DEGs (Fig. 8). "Energy metabolism", "amino acid metabolism", "metabolism of terpenoids and polyketides" and "biosynthesis of other secondary metabolites" were the enriched categories for the DEGs (Table 2).     A Venn diagram analysis was performed on the DEG sets (0 hpi vs 6 hpi and 6 hpi vs 24 hpi) (Fig. 9). We analyzed the intersection of the DEGs between different time points to investigate the DEGs in the resistant variety over time; this analysis further helped our understanding of the response to C. cassiicola infection in cucumber. In total, 146 DEGs were included in the intersection, and these genes are currently the focus of our future studies (Supplementary Table S2). According to these experiments, we selected and classified certain candidate genes that may be associated with the response to C. cassiicola in cucumber ( Table 3).

Verification of candidate genes by qRT-PCR.
To validate the reliability of the transcriptome data and analyze the expression of stress-responsive genes, 17 candidate genes were selected for qRT-PCR assays. We conducted experiments at different time points (i.e., 0, 3, 6, 12, 24, 48 and 72 hpi) using the resistant and susceptible varieties (Fig. 10). The results of the qRT-PCR were consistent with the transcriptome data. Thus, the expression   Table 3. Candidate genes may be associated with response to C. cassiicola in cucumber.  peroxidase (Cucsa.153390) in the resistant and susceptible varieties were nearly identical, and after 12-24 hpi, their abundance was at a very high level. The qRT-PCR results showed that secondary metabolism-related genes were strongly induced after the C. cassiicola infection. There was no significant difference in the expression of these candidate genes in both varieties at 0 hpi. However, the expression of most of the candidate genes changed significantly at 6-24 hpi. In addition, these genes showed stronger changes in the resistant variety than in the susceptible variety during the early stage of infection. In short, the resistant variety could respond more rapidly to C. cassiicola stress and induce more effective defenses.
Annotations of known miRNAs. To identify the known miRNAs, the clean reads were used in the BLAST search against known mature miRNAs and pre-miRNAs of miRBase (Version 21.0). In total, 64, 61 and 59 miR-NAs with a high sequence similarity to known miRNAs were identified at 0 hpi, 6 hpi and 24 hpi, respectively ( Table 4). Analyses of the known conserved miRNAs were performed by comparisons to major dicotyledonous species in miRbase (Version 21.0). The results showed that the known miRNAs were relatively conserved in cucumber and Arabidopsis thaliana (Supplementary Table S3).

Identification of DEMs and their targets.
According to the difference in expression (|fold change| > 2) in the multi-sample group comparison, we screened the differentially expressed mature miRNAs (DEMs). The miRNAs are primarily bound to target sites by complementary pairing, and we used genomic information and transcriptome data to predict the target genes in cucumber. In total, 150 target genes in 34 DEMs were predicted (Supplementary Table S4). Several miRNAs and their target genes that may be involved in the response to C. cassiicola were tested by qRT-PCR at different stages of pathogen invasion in the resistant and susceptible cucumber lines ( Table 5). The qRT-PCR results validated that the expression patterns of the miRNAs and target genes were consistent with the RNA-seq data (Fig. 11). The expression of miR164d and miR395c in both varieties increased noticeably after the pathogen infection. The expression of miR167e was first decreased at 3 hpi but significantly increased at 6 hpi in Jinyou 38. The accumulation of miR171f and miR172c was induced at 3 hpi in Jinyou 38, however, the expression of miR171f and miR172c increased in Ludixianfeng at 3 hpi. The expression of miR390a was significantly increased at 48-72 hpi in both varieties. In the resistant variety, miR396b was markedly up-regulated at 6 hpi, while in the susceptible variety, miR396b obviously increased at 24 hpi. The expression of miR408 obviously increased at 3-6 hpi in Jinyou 38 and was significantly up-regulated at 48-72 hpi in both varieties. Compared with the positive control (0 hpi), the expression of NAC domain containing protein (Cucsa.040380), which is cleaved by miR164d, was significantly up-regulated at 6, 24, 48 and 72 hpi in Jinyou 38 and reached a peak at 48 hpi in both varieties. After the pathogen infection, auxin response factor (Cucsa.047990), which is cleaved by miR167e, was always increased in both varieties, except for the susceptible variety at 24 hpi. The expression of GRAS family transcription factor (Cucsa.320850), which is cleaved by miR171f, was noticeably reduced after the pathogen inoculation in both varieties, while it was markedly up-regulated at 6 hpi in the resistant variety. Ethylene-responsive transcription factor (Cucsa.165940), which is targeted by miR172c, significantly decreased at 12 hpi and noticeably increased at 24 hpi in the resistant variety. The expression of leucine-rich repeat receptor-like protein kinase family protein (LRR-RLK) (Cucsa.164200), which is cleaved by miR390a, was significantly down-regulated after the pathogen infection but recovered to the original level at 12 hpi in the resistant variety; however, LRR-RLK was significantly decreased at  24 and 48 hpi in the susceptible variety. ATP sulfurylase (Cucsa.254710), which is targeted by miR395c, was significantly increased after 24 hpi in Jinyou 38, while in Ludixianfeng, its expression obviously increased at 48 and 72 hpi. The expression of anthranilate phosphoribosyltransferase (Cucsa.098530), which is cleaved by miR396b, noticeably increased at 3 and 6 hpi and then recovered to the normal level in Jinyou 38 but fluctuated slightly after the C. cassiicola infection in Ludixianfeng. Plantacyanin (Cucsa.077170), which is targeted by miR408, was markedly up-regulated after 6 hpi and gradually decreased in the resistant variety, while in the susceptible variety, its expression fluctuated, and the main increments occurred at 12, 48 and 72 hpi. Similarly, most of the miR-NAs and targets changed more significantly in the resistant variety than in the susceptible variety at 6-24 hpi. This finding was also consistent with the previous finding that the resistant variety could respond more rapidly to C. cassiicola stress.

Identification of novel miRNAs and their targets.
We used the MIREAP platform to predict new small miRNAs and the RNAfold web server to describe their secondary structures. The previously validated DEMs were used as the target genes to identify potential novel miRNAs, and we identified 7 novel potential miRNAs ( Table 6). The secondary structures predicted for these candidate novel miRNAs are shown in Fig. 12. These novel miRNAs were also tested by qRT-PCR, and their expression patterns were analyzed with their targets (Fig. 13).
In the resistant variety, the main reduction in Novel-miR1 occurred at 6-24 hpi, while in the susceptible variety, its expression was noticeably reduced after the pathogen infection. In the resistant variety, Novel-miR2 was significantly increased after the C. cassiicola infection and then recovered to the normal level at 48 hpi; however, in the susceptible variety, Novel-miR2 was significantly decreased at 24 hpi and increased at 72 hpi. The expression of Novel-miR3 was significantly decreased at 24 hpi but markedly increased at 72 hpi in Jinyou 38, while in the susceptible variety, its expression was obviously decreased after 24 hpi and then significantly increased at 72 hpi.

Discussion
In plants, reactive oxygen species (ROS) accumulation and lignin deposition are two early defense responses to pathogen infection 29 . In a previous study, ROS began to increase in the leaves of both cucumber varieties at 6 hpi, and a significant increase was observed at 24 hpi. During the early stage of pathogen infection (6-12 hpi), the ROS content in the resistant variety was relatively higher, indicating that the resistant variety could have a significant hypersensitive response to resist the invasion of pathogens more effectively. At 12 hpi, both cucumber varieties showed lignin deposition, which gradually increased over time. Similarly, the lignin deposition was more obvious in the resistant variety. These results suggested that the resistant variety responds earlier to C. cassiicola than the susceptible variety. According to the ITS sequence PCR results, the target bands were detected in both varieties at 6 hpi. In summary, we conclude that the response of cucumber to C. cassiicola occurs during the early stage of infection. Therefore, we used the resistant variety for RNA-seq to search for genetic resources related to C. cassiicola resistance, and samples at 0, 6 and 24 hpi were used for sequencing. We analyzed the transcriptomes and miRNAs in the cucumber cultivar Jinyou 38 after inoculation with C. cassiicola. In this study, we identified many DEGs in the two groups (0 hpi vs 6 hpi and 6 hpi vs 24 hpi). According to the results of GO term and KEGG pathways enrichment analysis, the secondary metabolism-related genes were strongly induced in Jinyou 38 after C. cassiicola infection (Figs 7 and 8). There were 146 DEGs in the intersection of the two groups (0 hpi vs 6 hpi and 6 hpi vs 24 hpi). To further explore the gene expression pattern differences between the resistant and susceptible cucumber cultivars challenged with C. cassiicola, we selected 17 candidate genes for qRT-PCR assays (Fig. 10).

MicroRNA
Mature  Table 6. Novel miRNAs and their targets. Plant hormones are active substances that can regulate the physiological response to environmental factors 30 . In plants, polyamine plays a role in defense signal transduction after pathogen infection 31 . S-adenosylmethionine decarboxylase (Cucsa.089960) is a key regulatory gene in the polyamine biosynthesis pathway 32,33 . The expression of S-adenosylmethionine decarboxylase was significantly up-regulated at 3 hpi, but no significant difference in expression was observed between the varieties. Gibberellin 3-oxidase (Cucsa.004570) can regulate the bioactivity of gibberellin, which plays an important role in plant responses to biotic and abiotic stress 34 . In the resistant variety, the expression of gibberellin 3-oxidase (Cucsa.004570) and gibberellin regulated protein (Cucsa.343030) was up-regulated earlier and to a higher degree. Thus, the gibberellin signaling pathway may be related to the response to C. cassiicola in cucumber.
MYB transcription factors play a key role in the plant defense response and can activate certain plant PR genes, leading to an increase in resistance to stress 35 . The expression of MYB transcription factor (Cucsa.121500) significantly increased after 24 hpi in both varieties, but during the early infection stage, its expression was higher in the resistant variety.
C. cassiicola is an obligate fungus that requires host energy for its nutrition; therefore, plants may inhibit the spread of pathogens by reducing energy metabolism. During the early stage of infection, phosphoenolpyruvate carboxylase (Cucsa.147540) and carbonic anhydrase (Cucsa.339710) were significantly decreased. Both genes are associated with photosynthesis 36,37 , supporting the hypothesis that hosts may limit energy metabolism to resist the pathogen. The calcium signaling pathway plays an important role in plant stress resistance because Ca 2+ can activate and regulate the expression of downstream genes as a secondary messenger 38 . Calmodulin-like is a key factor in plant-pathogen interactions. When a plant is subjected to pathogenic stress, pathogen-associated molecular patterns (PAMPs) produce specific Ca 2+ signals, and calmodulin-like mediates NO signaling to induce the plant hypersensitive response (Fig. 14). The significant increase in calmodulin-like (Cucsa.254730) expression first occurred at 6 hpi in the resistant variety, while in the susceptible variety, the increase occurred at 24 hpi. Therefore, calmodulin-like in the resistant cultivar could respond to stress and regulate the downstream pathway earlier.
We also identified several defense-related genes, most of which were associated with the cell wall and reactive oxygen scavenging. Xyloglucan endotransglucosylase/hydrolase family protein (Cucsa.092350) is a cell wall-modified protein with the ability to loosen cell walls 39 . These genes have been recently shown to be involved in multiple growth regulation processes, including cell extension and the stress response 40,41 . In this study, these genes responded earlier to stress in the resistant varieties. We speculate that these genes may provide resistance to pathogen stress by altering the toughness of the cell wall. Chitinase can hydrolyze chitin in the cell wall of the pathogen, and the digested oligomeric product can also induce further defense as a signal molecule 42 . Unfortunately, the chitinase (Cucsa.152090) expression patterns were similar in both varieties. The PR protein is the terminal gene of the plant defense pathway  and can induce resistance in the response to pathogens 43 . Thaumatin is a protein with multi-biological activity and important function in the plant defense class of sweet protein, belonging to the PR proteins 44 . The eukaryotic aspartyl protease family protein is an important hydrolase that is used to synthesize specific disease resistance-related proteins and hydrolyze the proteins secreted by invading pathogens 45 . Under pathogen stress, many ROS in plant outbreaks can improve plant resistance but can also destroy the organizational structure of the plant. Plants often produce peroxidase to remove excess ROS, thereby relieving the damage caused by oxidative stress to the plant 46 . After 12-24 hpi, the expression of PR protein (Cucsa.133220), thaumatin (Cucsa.302870), aspartyl protease family protein (Cucsa.043900) and peroxidase (Cucsa.153390) in both varieties was nearly identical, and after 12-24 hpi, their abundances were at a very high level. We also speculate that these genes play an important role in resistance to C. cassiicola, although no difference was observed between the varieties.
Secondary metabolites are key factors inducing systemic resistance 47 . The pleiotropic resistance protein is an ATP-binding cassette (ABC) type transporter that plays an important role in the regulation of secondary metabolism and environmental adaptation 48 . The significant accumulation of pleiotropic resistance protein (Cucsa.017550) in the resistant variety appeared at 24 hpi, but in the susceptible variety, this accumulation occurred at 48 hpi. Cytochrome P450 can catalyze certain secondary metabolites, such as indole, phenylpropane, phytohormones and alkaloids, to improve plant immunity 49,50 . In our study, the expression of cytochrome P450 CYP2 (Cucsa.342000) was higher in the resistant variety, particularly at 24 hpi. Phenylalanine ammonia-lyase is closely related to plant stress resistance as follows: phenylalanine ammonia-lyase acts as the initial enzyme in flavonoid biosynthesis and promotes the synthesis of flavonoids, thereby allowing plants to resist external stress 51 . In addition, 4-coumarate: CoA ligase is an important enzyme in the phenylpropanoid pathway in plants that affects the synthesis of lignin 52 . These two genes belong to the pathway of phenylpropanoid synthesis. As shown in Fig. 15, phenylalanine ammonia-lyase (PAL) first catalyzes phenylalanine to cinnamic acid. Cinnamic acid is involved in the synthesis of ubiquinone and can be catalyzed by 4-coumarate: CoA (4CL) ligase to form cinnamoyl-CoA. Cinnamoyl-CoA is then involved in the synthesis of lignin and flavonoids. The expression patterns of phenylalanine ammonia-lyase (Cucsa.124480) and 4-coumarate: CoA ligase (Cucsa.261210), which are the two key genes in this pathway, were very similar. After 12 hpi, the expression levels of these genes were significantly higher in the resistant variety. Therefore, under the stress of C. cassiicola infection, secondary metabolism in cucumber plays an important role in disease resistance.
The interactions among the miRNAs and their targets can offer a complex mechanism for the response to stress 23,53 . In this study, miRNAs were also sequenced to explore the role of the miRNAs in the regulation of cucumber mRNAs under C. cassiicola infection stress. In this study, we identified 34 DEMs belonging to 17 families and 7 novel miRNAs with precursors using high-throughput sequencing techniques and bioinformatics analysis. We also predicted many targets for the miRNAs in combination with the transcriptome data. The qRT-PCR results of the candidate miRNAs and targets showed that their expression patterns differed between the cucumber varieties (Fig. 11). In our study, most candidate miRNAs, including miR164d, miR167e, miR171f and miR172c, target various transcription factors, such as NAC, ARF, SCL and AP. In this assay, miR164d targets NAC domain containing protein (Cucsa.040380). NAC is a plant-specific transcription factor family that is not only involved in the process of growth and development but also participates in stress responses 54 . In Jatropha curcas, NAC can improve plant resistance by regulating hormone signaling. The overexpression of the NAC gene can significantly improve the tolerance of Jatropha curcas to Pseudomonas syringae 55 . Several studies have shown that NAC transcription factors are related to lignin synthesis 56,57 . The target gene for miR167e was ARF, which can bind to auxin response elements. ARF can activate the downstream hormone signaling pathway to promote the expression of defense genes in plants. These expression products can be transmitted as signal molecules to allow plants to perceive signals and respond 58 . We predicted that the GRAS family transcription factor (Cucsa.320850) is cleaved by miR171f. GRAS is a plant-specific transcription factor that plays an important role in signal transduction and hormone regulation 59 . Our results show that the target of miR172c is APETALA, which is an ethylene-responsive transcription factor. APETALA plays a very important role in plant growth and development, including the development of flower organs, fruit development, seed formation and stress response 60 . Due to its importance in plants, the regulation of miR172c to the target gene APETALA is essential for plants to resist pathogens and improve their resistance. The target gene for miR390a is LRR-RLK (Cucsa.164200). LRR-RLK proteins constitute a variety of transmembrane receptors involved in biological functions in plants, such as growth and development and response to biological and abiotic stress 61,62 . Unfortunately, in our assays, the LRR-RLK proteins may be inhibited by miR390a; thus, their expression declined in both varieties. We hypothesize that the expression of miR390a can be reduced to increase the expression of LRR-RLK, which may help improve the resistance of the host to C. cassiicola. The target gene of miR395c is ATP sulfurylase (Cucsa.254710), which is the first enzyme in the process of sulfate assimilation. Sulfur is not only involved in the synthesis of sulfur-containing amino acids and proteins but is also an important component of many enzymes and cofactors. Sulfur is also related to metabolic substances produced by plants under stress 63,64 . miR396b is a conserved class of miRNAs in plants. In this assay, the target of miR396b is anthranilate phosphoribosyltransferase (APE), which is related to tryptophan synthesis 65 . In plants, the endogenous jasmonic acid (JA) biosynthesis, plant signal transduction pathway and terpenoid indole alkaloid (TIA) pathway were triggered by tryptophan, which can regulate plant responses to stress 66 . Plantacyanin (Cucsa.077170) is the target gene of miR408. In chloroplast photosynthesis, plantacyanin is a copper-containing protein that serves as an electron transporter. Plantacyanin is indirectly involved in the removal of ROS in plants; therefore, plantacyanin is closely related to plant resistance 67 .
We identified seven novel miRNAs that targeted several previously validated DEMs (Table 6). These targets, such as 4-coumarate: CoA ligase, pleiotropic drug resistance protein, and phenylalanine ammonia-lyase, are mainly associated with secondary metabolism. Since Novel-miR1 and Novel-miR7 can cleave 4-coumarate: CoA ligase and phenylalanine ammonia-lyase separately, we hypothesized that if Novel-miR1 and Novel-miR7 were decreased, the expression of genes in the phenylpropanoid and lignin biosynthesis pathway would be promoted, which can greatly improve the resistance of cucumber to C. cassiicola.
The experimental results showed that the candidate miRNAs are mainly involved in the regulation of secondary metabolism-related genes, so the secondary metabolism-related defense responses in resistant variety occured earlier and faster. These results combined with our previous data indicated that secondary metabolism plays an important role in cucumber resistance to C. cassiicola infection. Based on the above analysis, we summarized these regulatory pathways related to secondary metabolism (Fig. 16). Secondary metabolism was regulated by the interaction of these miRNAs and targets, which could affect the resistance of cucumber to TLS.

Conclusions
The response time is critical for cucumber resistance to C. cassiicola infection, as the results of ROS and lignin staining suggested that the resistant variety Jinyou 38 showed a faster response to pathogen stress. At the early stages of infection, we obtained many DEGs from the resistant variety that were involved in the response to C. cassiicola. The functional analyses showed that these DEGs were mainly concentrated in secondary metabolism. We performed qRT-PCR assays on DEGs related to plant hormone, transcription factor, metabolic correlation, Ca 2+ signaling pathway, secondary metabolism and defense genes. The results showed that the resistance of Jinyou 38 may be mainly attributed to secondary metabolism. In addition to transcriptome sequencing, this study also sequenced known miRNAs and predicted novel miRNAs in cucumber that contribute to the interaction between the host and C. cassiicola. By analyzing the functions of these targets, we found that several known and novel miRNAs may mediate secondary metabolism. Combined with the transcriptome and miRNA data, the results of this study were used to summarize the mechanism of secondary metabolism in the resistance of cucumber to C. cassiicola. Overall, these results provide a new theoretical basis for studies investigating resistance mechanisms and key resistance genes against cucumber TLS at the molecular level.

Methods
Plant growth and C. cassiicola inoculation. The cucumber varieties used in the experiments were Jinyou 38 (preliminary experiments showed that this variety was resistant to TLS) and Ludixianfeng (TLS susceptible), which were planted in a greenhouse at 28 °C under 16:8 light/dark cycles. C. cassiicola was obtained from the Agricultural Culture Collection of China. The true leaves of 6-week-old plants were sprayed with spore suspensions (2 × 10 5 sporangia/ml). Then, the inoculated materials were kept at 100% relative humidity to ensure spore germination. The inoculated leaves were harvested at seven time points (i.e., 0, 3, 6, 12, 24, 48 and 72 h post-inoculation), frozen in liquid nitrogen and stored at −80 °C. All experiments were performed in three independent replicates for each time point. PCR assay for the detection of C. cassiicol. The DNA of the inoculated leaves was extracted using the Plant Genomic DNA Kit (Tiangen, Beijing, China). C. cassiicola internal transcribed space (ITS) primers were used to detect whether the pathogen entered the host. For the qPCR analysis, the following cycling parameters were used: 1 hold at 94 °C for 5 min, 35 cycles at 94 °C (30 s), 60 °C (15 s), 72 °C (1 min), and 72 °C (10 min). The PCR products were detected by agarose gel electrophoresis.
Histochemical staining. Lignified tissues were stained with 1% phloroglucinol solution (Maclin, Shanghai, China). First, the plant leaves were soaked in stationary liquid containing 95% ethyl alcohol and glacial acetic acid (1:1, v/v) for 24 h, rinsed with distilled water, vacuum treated in saturated aqueous solution of chloral hydrate for 10 min, and stored at room temperature for two or three days until they became transparent.
DAB staining was performed to detect H 2 O 2 in the inoculated plants. The cucumber leaves were soaked in a DAB (Sigma, MO, USA) solution (1 mg/ml) for 8 h and placed on the shaker (80~100 rpm) for 4-5 h. The leaves were added to the decoloring solution containing absolute ethyl alcohol, glacial acetic acid and glycerol (3:1:1, v/v/v) and heated until the tissues were transparent for observation.
The O 2 − of the infected leaves was observed by performing NBT staining. The samples were infiltrated with a 0.1% NBT (Ameresco, OH, USA) solution for 4 h and subjected to a vacuum treatment for 15 min. Decolorization is essentially the same as that described above. mRNA-seq and analysis. The three time points of the inoculated leaves were defined as follows: 0, 6, and 24 hpi. Leaves at the three time points were selected for transcriptome sequencing using an Illumina NextSeq 500 platform (Personal Biotechnology, Shanghai, China). The RNA of the three samples at each time point were pooled in equal volumes (the concentrations were adjusted to 10 nM).
The transcripts were annotated via BLAST against the cucumber genome database (https://phytozome.jgi.doe. gov/pz/portal.html). The expression of each gene was calculated according to the reads per kilo bases per million reads (RPKM). The expression level was calculated according to the baseMean value, which is the sequencing depth of each transcript normalized to the library size. The baseMean value and P-value of the genes were calculated by the DESeq platform. Genes with expression changes greater than twofold between two points were considered DEGs (P-value < 0.05). All DEGs were analyzed with eggNOG functional categories (http://eggnog.embl.de/), GO enrichment (http://www.geneontology.org/) and KEGG enrichment (http://www.genome.jp/kegg/) 68-71 . MicroRNA-seq and analysis. The RNA samples used for the miRNA-seq were the same as those used for the mRNA-seq. The sequencing of the miRNAs was performed by Personalbio using the Illumina NextSeq500 system. The expression quantity was calculated according to the transcripts per million (TPM). To identify the known miRNAs, the clean reads were BLASTed against the miRNA database miRbase 21.0 (http://www.mirbase.org/). Because miRBase contains few cucumber miRNA sequences, we used melon data as a reference. We also analyzed the miRNA families and conservation in miRbase. The DEMs were screened according to the difference in expression |fold change| >2). The RNAfold web server (http://nhjy.hzau.edu.cn/kech/swxxx/jakj/dianzi/Bioinf4/ miRNA/miRNA1.htm) was used to map the pre-miRNAs. For sequences that were not annotated to any information, we predicted new miRNAs using mireap (http://sourceforge.net/projects/mireap/). Since the miRNAs are primarily bound to the target site by complementary pairing, based on the transcriptome data, the targets of the mature miRNA sequences were identified using miranda 72 (http://www.microrna.org/microrna/home.do).
Quantitative real-time RT-PCR assay. RNA was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) and synthesized into cDNA using the FastQuant RT Kit (with gDNase) (Tiangen, Beijing, China). Quantitative real-time RT-PCR (qRT-PCR) was conducted using a SuperReal PreMix Plus Kit (SYBR Green) (Tiangen, Beijing, China). The qRT-PCR was performed on a LightCycler 480 (Roche Molecular Systems, CA, USA). Cucumber Actin was used as the reference gene to normalize the data. To validate the presence and relative expression of the miRNAs obtained from the sequencing, the known and novel miRNAs were assayed by qRT-PCR. The reverse transcription reaction was performed using a miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). The qRT-PCR was performed using the miRcute miRNA qRT-PCR Detection Kit (SYBR Green) (Tiangen, Beijing, China), and U6 snRNA was used as the internal control. The relative expression was calculated according to the 2 −ΔΔCt method 73 , and the standard deviation was calculated using three biological replicates. The primers used in the experiment were synthesized by GENEWIZ (Suzhou, China) and are listed in Supplementary Table S5.
Data availability. The mRNA raw data were deposited in the NCBI Sequence Read Archive (SRA) with the accession number SRP117262; the small RNA raw data were deposited in the NCBI Sequence Read Archive (SRA) with the accession number SRP117230.