The resistant tomato line of Cf-12-tomato (CGN7495) and the susceptible line Moneymaker were obtained from the Chinese Academy of Agricultural Sciences (Beijing, China). They were grown in a greenhouse at the Horticultural station of Northeast Agricultural University (Harbin, China). The growing condition was 16-h light and 8-h darkness at 25°C with an ambient humidity of 95%. C. fulvum physiological race 1.2.3 was acquired from tomato-growing regions in Harbin using the single sporangiophore transfer method, as previously described (Hubbeling, 1971; Wong and Wilcox, 2000); and was further axenically propagated on Moneymaker. At the four-five leaf stage, the abaxial surface of the Cf-12-tomato seedlings was inoculated with a suspension of 1 × 10⁷ sporangia per ml (Wang et al., 2007). Leaf samples were harvested at 0–15 days of post-infection (dpi) for microscopic analysis.

In order to observe the interaction process of C. fulvum on Cf-12 tomato, the lactophenol trypan blue staining method was carried out according Franco's approach (Franco et al., 2008). The fungal tissue and dead host cells would be densely stained, while living host cells would impart a translucent and slightly brown color. The leaf samples were harvested at 0–15 dpi, immediately stained, clarified overnight in chloral hydrate solution (2.5 mg/ml) (Keogh et al., 1980), and examined using an EVOS® microscope (ThermoFisher, USA) and an OLYMPUS SZX10 dissecting microscope (Olympus, Japan).

Total RNA from Cf-12-tomato leaves was extracted and analyzed, as previously described (Fang et al., 2015). The integrity of the isolated RNA was calculated, as previously described (Schroeder et al., 2006); and samples with an RNA integrity number >9.0 were used for libraries construction. The libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA). The clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumia). After cluster generation, the libraries were sequenced on an Illumina Hiseq 4000 platform conducted by the Novogene Bioinformatics Institute (Beijing, China); and 150-bp paired-end reads were generated.

Raw sequence data were processed using Perl scripts (http://www.perl.org), developed by Novogene Bioinformatics Institute (Beijing, China), in order to remove reads that contained adapter fragments and ploy-N stretches (the number of ploy-Ns is >10%). Phred quality scores (Q20: ratio of an error rate 1%, Q30: ratio of an error rate 0.1%) and GC-content were calculated, and only the data with a quality score (Q_phred) ≥ 30 (Q30) were used for further analyses.

The reference genome and gene model annotation files were downloaded from the Ensembl Genomes Databases (ftp://ftp.ensemblgenomes.org/pub/release-23/plants/fasta/solanum_lycopersicum/dna/). An index of the reference genome was built using Bowtie v2.2.3 (Langmead and Salzberg, 2012), and paired-end clean reads were aligned to the reference genome using TopHat2 v2.0.12 (Kim et al., 2013). Cufflinks v2.1.1 (Trapnell et al., 2010, 2012) was used to construct and identify both known and novel transcripts from TopHat2 alignment results.

HTSeq v0.6.1 (EMBL, Heidelberg, Germany) was used to count the read numbers mapped to each gene. Gene expression levels were calculated based on the length of the gene, sequencing depth and read count mapped to this gene using the Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM; Trapnell et al., 2010) method. Genes with FPKMs in intervals of 1–3, 3–15, 15–60, and beyond 60 were considered to be expressed at low level, medium level, high level and very high level, respectively.

After calculating the gene expression levels, a differential expression analysis of two conditions or groups was performed using the DESeq R package (1.18.0) (Anders and Huber, 2010). DESeq provides statistical strategies to determine differential gene expression using the negative binomial distribution model (K_ij ~NB[μ_ij, σij2]). The resulting P-values and fold-changes were adjusted using the Benjamini-Hochberg's approach to control the false discovery rate (FDR ≤ DRe5; Benjamini and Hochberg, 1995). Genes with an adjusted P-value (padj) < 0.05 were considered to be significant DEGs.

Eighteen DEGs involved in plant disease resistance pathways were validated using quantitative real-time PCR (qRT-PCR). The primer pairs of the selected genes were designed using Primer Premier 6.0 (Premier Biosoft, Canada; Table S1). The tomato actin gene (U60478.1) was used as a reference control. The qRT-PCR was performed using AceQ® qPCR SYBR® Green Master Mix (Vazyme, USA) on an iQ™ 5 Multicolor Real-time PCR Detection System (Bio-Rad, USA). The reaction parameters were as follows: 95°C for 7 min, and 40 cycles of 95°C for 10 s, 58°C for 30 s, and 72°C for 20 s. Each sample was repeated three times, and relative expression levels were evaluated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001).

The GO enrichment analysis of DEGs was performed using the GO-seq based on Wallenius non-central hyper-geometric distribution (Young et al., 2010), in which gene length bias in DEGs were adjusted by Bonferroni correction. GO terms with a corrected P < 0.05 were considered significantly enriched by DEGs. KEGG pathway enrichment (http://www.genome.jp/kegg/) was used to identify significantly enriched signal transduction pathways or metabolic pathways in DEGs. Significantly enriched DEG pathways were identified using the KOBAS 2.0 software (KOBAS, Surrey, UK; Xie et al., 2011), and adjusted by hyper-geometric test and Benjamini-Hochberg FDR correction (FDR ≤ F.05).

Light microscopy was used to observe the interaction process between C. fulvum and Cf-12 tomato, or the Moneymaker leaves. A representative image of the C. fulvum mycelium and spores are shown in Figure 1A. Our results revealed that conidiospores germinated at 2 or 3 dpi (Figure 1B), the hypha started penetrating into the stomata at 4 dpi in both the Moneymaker and Cf-2 tomato (Figures 1C,D). No difference was observed between the two cultivars at this stage. Then, the hypha penetrated into the intercellular space and continued from the substomatal cavity into the intercellular space between the spongy mesophyll cells (Thomma et al., 2005), emerging through the stomata at 8 dpi and gradually plugging the stomata at 10 dpi on the Moneymaker cultivar (Figures 1E,F). However, cells surrounding the stomata appeared to die (necrosis) at 8 dpi on Cf-12 tomato (Figure 1G); and a large number of necrotic spots began to appear at 10 dpi (Figure 1H). The area of necrosis was gradually enlarged at 12–15 dpi (Figures 1I,J). In the advanced stages of disease development, the hyphae formed a thick, gray-white mold layer on the abaxial surface of the leaves at 15 dpi on the Moneymaker cultivar (Figure 1K) and formed yellow necrotic spots on the front side of the leaf on Cf-12 tomato (Figure 1L). This observation is in agreement with the results of a previous study (Thomma et al., 2005). Based on this observation, we selected Cf-12 tomato samples at 0, 4, and 8 dpi for mRNA-Seq and qRT-PCR analysis.

In order to investigate differences in transcriptome between C. fulvum-infected Cf-12 tomatoes (Cf12_B, 4 dpi; Cf12_C, 8 dpi) and non-infected controls (Cf12_A, 0 dpi), RNA from these three samples were sequenced. An average of 59,291,565, 56,131,442 and 70,822,873 raw reads from Cf12_A (0 dpi), Cf12_B (4 dpi), and Cf12_C (8 dpi) were generated, respectively (Table 1, Figure S1). After removing adaptors, low quality reads, duplications and ambiguous reads, an average of 60 million clean reads were obtained for each sample; and 94% of these clean reads were at the Q30 (ratio of error rate ≤ 0.1%) level (Table 1). Furthermore, at least 89% of these clean reads were mapped; of which, ~88% were uniquely mapped and 0.7% were multiple mapped to tomato chromosomes, respectively (Table S2). Pearson correlation coefficient (R²) analysis of the FPKM distribution between different biological replicates revealed a high level of reproducibility of RNA expression patterns (R² = 0.98–0.99, P < 0.001; Table 2, Figure S2). A large majority of reads were mapped to the genome exon regions (Figure S3), and reads density was positively correlated to chromosome length (Figure S4). Cufflinks analysis revealed a total of 15,395 transcripts from the three samples (Cf12_A, Cf12_B, and Cf12_C), including 14,494 known transcripts (Table S3) and 901 new transcripts (Table S4).

DEGs were identified using the DESeq software with a padj <0.05. For Cf12_B vs. Cf12_A, 9100 DEGs were detected, including 4080 upregulated and 5020 downregulated DEGs (Table S5A, Figure 2); for Cf12_C vs. Cf12_A, 8643 DEGs were detected (3999 upregulated and 4644 downregulated DEGs; Table S5B, Figure 2); and for Cf12_C vs. Cf12_B, 2547 DEGs were identified (1729 upregulated and 818 downregulated DEGs; Table S5C, Figure 2). In addition, 736 DEGs were shared among the three groups (Cf12_B vs. Cf12_A, Cf12_C vs. Cf12_A, and Cf12_C vs. Cf12_B). Approximately 1500 DEGs were shared between Cf12_C vs. Cf12_A and Cf12_C vs. Cf12_B, or Cf12_B vs. Cf12_A and Cf12_C vs. Cf12_B. However, 6426 DEGs were shared between Cf12_B vs. Cf12_A and Cf12_C vs. Cf12_A. This further suggests that a common group of genes were activated or deactivated upon C. fulvum infection (Figure 3). In order to observe the overall changes of the quantity of gene expression, a hierarchical clustering of DEGs using FPKM analysis was performed. Results revealed that after the inoculation of C. fulvum on Cf-12 tomato, genes with low expression quantity in Cf12_A increased in expression (in both Cf12_B and Cf12_C). On the contrary, many genes with high expression quantity in Cf12_A decreased in expression in Cf12_B or Cf12_C (Figure 4). Changes in gene expression quantity were also observed between Cf12_B and Cf12_C (Figure 4), suggesting that the response of Cf-12 tomato to C. fulvum infection changes as time progresses.

In order to identify these similar expression patterns, the relative expression level of DEGs were analyzed by K-means clustering algorithm (Hartigan and Wong, 2013). Clustering analysis revealed that six expression patterns (subclusters) of DEGs were identified (Figure 5). The most prominent group was subcluster_5, in which 5343 genes were upregulated after C. fulvum infection. A similar pattern was observed in subcluster_4, where genes revealed a higher expression level in Cf12_B and Cf12_C; however, the number of genes was much lesser than that in subcluster_5. Subcluster_1, subcluster_2 and subcluster_6 revealed a similar expression pattern, in which most of the genes were downregulated in Cf12_B and Cf12_C; but the downregulated level was relatively small. Subcluster_3 contained 41 genes that were downregulated in Cf12_B, but were upregulated in Cf12_C. These dynamic gene expression patterns further suggest that Cf-12 tomato was resistant to C. fulvum via a highly complex process.

In order to validate the RNA-Seq data, qRT-PCR was performed through 18 DEGs using three biological replicates. These 18 genes were selected to reflect some of the functional categories and pathways described below (Section GO and KEGG Enrichment Analysis of DEGs), such as the plant-pathogen interaction pathway and plant hormone signal transduction pathway (Table S1). These qRT-PCR results were compared with the RNA-Seq data. As shown in Figure 6, the trends of these gene expression patterns were consistent and had a strong positive correlation coefficient (R² = 0.9619), indicating that the RNA-Seq data was reliable.

In order to further characterize the functions of DEGs, GO enrichment analysis was performed using GOseq. The top 10 enrichment terms of the biological process, cellular component and molecular function were selected, respectively, as the main nodes of the directed acyclic graph. In the biological process category, significant terms were enriched in the phosphorus metabolic process, protein phosphorylation biological process, cellular protein modification process, salicylic acid metabolic process, response to jasmonic acid, and salicylic acid biosynthetic process. These processes were generally associated with disease resistance (Figure S5A). Components associated with membranes such as thylakoid, plastid thylakoid membranes, and photosynthetic membranes (Figure S5B) were significantly enriched in the cellular component category. In the molecular function category, transferase activity, phosphotransferase activity and protein serine/threonine kinase activity were enriched in catalytic activity; and purine nucleoside binding, anion binding, nucleoside binding, adenyl nucleotide binding and ATP binding were enriched in binding (Figure S5C). These binding terms and catalytic activity terms played an important role in signal recognition and signal transduction.

One interesting finding is that many transcription factor genes were significantly upregulated or downregulated in response to C. fulvum. A total of 406 transcription factor-related DEGs were identified and annotated into 60 families (Pérez-Rodríguez et al., 2010; Jin et al., 2014; Wang et al., 2016b; Table S6). These 60 transcription factor families were grouped into two clusters (Figure S6). The majority of DEGs in cluster I revealed an upregulation trend at the early stage (Cf12_B) and a downregulated trend at the later stages (Cf12_C). In cluster II, DEGs exhibited a higher expression at both stages or only at the later stage. For example, one of the MYB family transcription factors (Solyc05g053150.1) was upregulated to more than 5-fold in both Cf12_B vs. Cf12_A and Cf12_C vs. Cf12_A (Table S6). This protein, like most of the 406 transcription factors, remains uncharacterized; and is predicted to be a MYB23-like transcription factor. The downregulation of transcription factors was also identified. For example, the OFP family Solyc09g065350.1 (Wang et al., 2005), another uncharacterized transcription factor, was downregulated to ~5-fold in both infection stages. This result suggests that transcription factors play a critical role in response, especially in the early response, to C. fulvum infection.

The KEGG pathway enrichment analysis was used to identify the biological pathways of incompatible interaction. It was found that the plant-pathogen interaction pathway was significantly enriched (FDR ≤ 0.05). As shown in Figure 7, the number of genes and the rich factor are significantly higher than the other pathways. Many other disease-resistance pathways, including stilbenoid and gingerol biosynthesis, phenylpropanoid biosynthesis, ubiquinone and quinone biosynthesis, flavonoid biosynthesis, pentose phosphate pathway, and nitrogen metabolism were also enriched. In addition, it was found that the transduction plant hormone signals, the biosynthesis of unsaturated fatty acids, fatty acid metabolism, and carbon metabolism participate in Cf-12 tomato response to C. fulvum infection.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.