Plant tissue-cultured material from strawberry cultivars ‘Yanli,’ ‘Toyonoka,’ ‘Sachinoka’ and ‘Benihoppe’ was grown in a greenhouse at the Fruit Molecular Biology Laboratory of Shenyang Agricultural University. The experiment was organized into two groups, with 20 strains per variety in each group. The plants were selected for uniform size and color and an absence of visual defects. The plants were managed under identical conditions. Collected leaves were frozen immediately in liquid nitrogen and stored at -80°C until further use.

The pathogenic fungi used for the experiments were grown in shaking Potato Dextrose (PD) cultures for 7–10 days at a culture temperature of 25°C and then were filtered through four layers of gauze into a diluent. The accession number of the internal transcribed spacer (ITS) sequence of the C. gloeosporioides strain is KT224457. After incubation, the suspension population was counted using a bacterial counting chamber and adjusted to 10⁶ spores ml^-1. Then, the conidial suspension was transferred to a watering can, and the plant leaves and stems were sprayed to cover their surfaces with a uniform layer of the suspension. The inoculated leaves were subsequently incubated at 25°C for disease development. The results were recorded after 0, 24, 72, 96, and 120 h.

External symptoms of disease were indicated by any sign of disease spots on the plants. The disease incidence and disease index (DI) on each leaf were observed and recorded after incubation. When the spot zone beyond the wounded site on a leaf was greater than 1 mm wide, the site was scored as infected. The disease incidence was expressed as the percentage (%) of infected leaves among the total leaves. The disease index of each leaf was rated according to the percentage of observed spots on a scale of 0–5: (0) no disease spots; (1) <5% disease spot area; (2) 5–10% disease spot area; (3) 10–15% disease spot area; (4) 15–20% disease spot area; and (5) >20% disease spot area.

Total RNA was isolated from the plant leaves using the modified CTAB method as described by Chang et al. (2007). The RNA samples were treated with DNase (TaKaRa, Japan) for 4 h. The integrity of the RNA samples was examined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States).

For the RNA-seq experiments, the strawberry cultivars ‘Yanli’ (resistant) and ‘Benihoppe’ (susceptible) were used. The samples included 0-h-untreated ‘Yanli’ (R1), 72-h-untreated ‘Yanli’ (R2) 72-h-infected (with C. gloeosporioides) ‘Yanli’ (R3), 0-h-untreated ‘Benihoppe’ (S1), 72-h-untreated ‘Benihoppe’ (S2) and 72-h-infected ‘Benihoppe’ (S3). Leaves were randomly collected for each treatment. A total of six pools were submitted to the OEbiotech Company (Shanghai, China) for cDNA library preparation and sequencing. The paired-end library preparation and sequencing were performed following standard Illumina methods using a DNA sample kit (#FC-102-1002, Illumina). The cDNA library was sequenced on the HiSeq^TM 2500 Illumina sequencing platform.

Clean reads were preprocessed by removing the adaptor sequences and discarding the empty and low-quality sequences. The reads were then mapped to the octoploid strawberry genome FAN_r1.1¹ (Hirakawa et al., 2013) using Tophat with the default parameters.

We annotated genes based on a set of sequential BLAST searches (Altschul et al., 1997) to identify the most descriptive annotation for each sequence. The genes were compared with sequences in the National Center for Biotechnology Information (NCBI) non-redundant (Nr) protein databases² and the Swiss-Prot protein database³. The genes were also queried against the Plant Resistance Gene database⁴.

Gene expression levels in the RNA-seq analysis were measured as reads per kilobase of exon model per million mapped reads (RPKM). DESeq software⁵ was used to identify differentially expressed genes (DEGs) via pairwise comparisons, and the results of all statistical tests were adjusted for multiple testing with the Benjamini—Hochberg false discovery rate (FDR < 0.01). Sequences were deemed to be significantly differentially expressed if the adjusted P-value obtained by this method was <0.05 and at least a two-fold change (>1 or <-1 in log 2 ratio value) in the RPKM was observed between the two libraries.

The functions of the DEGs were characterized using Gene Ontology (GO) terms with the agriGO database (Du et al., 2010). The second method used for DEG characterization was based on a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Pathways with Q-values < 0.05 were considered significantly enriched in DEGs. The Q-value was used to determine the P-value threshold (Benjamini and Hochberg, 1995).

First-strand cDNA was synthesized from one microgram of total RNA using the PrimeScript RT reagent kit with gDNA Eraser (TaKaRa, Japan) according to the manufacturer’s instructions.

To determine the expression of DEGs, quantitative reverse transcription PCR (qRT-PCR) was performed on the iQ5 Real-Time PCR Detection System (Bio-Rad, United States) using SYBR Premix Ex Taq (TaKaRa, Japan). The cycling conditions were as follows: 40 cycles of denaturation at 95°C for 30 s, annealing at 58°C for 45 s, and extension at 72°C for 30 s. A dissociation curve was generated at the end of each PCR cycle to verify that a single product had been amplified using software provided with the iQ5 Real-Time PCR Detection System. Each gene was analyzed in triplicate, and then the average threshold cycle (CT) was calculated per sample. All of the data were normalized to the level of the Fv26S internal transcript control. Relative fold changes in gene expression were calculated using the comparative Ct method (2^-ΔΔCT) (Liu C. et al., 2012; Liu T. et al., 2012). Primers used in the PCR analysis are listed in Supplementary Table S3.

The data were statistically analyzed using DPS 7.05 software. After conducting an analysis of variance (ANOVA), Duncan’s test was used to identify significant differences between mean values.

In northern China, ‘Yanli,’ ‘Toyonoka,’ ‘Sachinoka’ and ‘Benihoppe’ are the main strawberry cultivars. The different cultivars responded differently to the pathogen C. gloeosporioides. We observed the leaves of ‘Yanli,’ ‘Toyonoka,’ ‘Sachinoka’ and ‘Benihoppe’ after 72 h of infection with C. gloeosporioides. Disease symptoms appeared in the strawberry cultivars (Figure 1A). The disease incidence and index values were determined at this time point, and statistical analysis was then performed (Figures 1B,C). The disease incidence and index values were the lowest for ‘Yanli’ and the highest for ‘Benihoppe’ among the four cultivars.

To identify the genes that are responsive to C. gloeosporioides infection in strawberry and compare gene expression patterns between resistant and susceptible cultivars, six libraries were constructed from Colletotrichum gloeosporioides-infected leaf tissues and from non-vaccinated leaf tissues as a control. A total of 16.12 G and 16.05 G reads were sequenced from the ‘Yanli’ (R1, R2, and R3) and ‘Benihoppe’ (S1, S2, and S3) cultivars, respectively, and 206433226 clean reads were obtained from the six libraries (Table 1).

The reads were compared with the octoploid strawberry genome FAN_r1.1.⁶ Gene expression levels in the RNA-seq analysis were measured as RPKM (Supplementary Table S1). The analyses revealed pairwise alterations of gene expression. R3 did not cluster with R1 and R2, and S3 did not cluster with S1 and S2 (Figure 2). The heat map, drawn from all gene count data from both genome references based on Euclidean distances, depicts the relationships of all of the transcriptomes.

The genes were annotated based on their similarity to known or putative sequences in web databases. Among the 230838 genes, 163033 (70.63%) had at least one significant match to a gene model according to the BLAST search (Table 2). The expressed genes were categorized according to GO and KEGG pathway analyses. All of the genes were queried against the Plant Resistance Gene database⁷ using BLASTp, and a total of 5227 putative R genes were identified. A list of the identified R genes is provided in Supplementary Table S2.

To identify the C. gloeosporioides-responsive genes in resistant and susceptible strawberry cultivars, R3 was compared with R1 and R2, with a union set of 2707 DEGs as C. gloeosporioides-responsive genes in ‘Yanli’ (Figure 3A), and S3 was compared with S1 and S2, with a union set of 1354 DEGs as C. gloeosporioides-responsive genes in ‘Benihoppe’ (Figure 3B). Then, the C. gloeosporioides-responsive genes were compared between ‘Yanli’ and ‘Benihoppe’ to reveal genes that were highly conserved or divergent in response to infection (Figure 3C). In total, we identified 2481 unigenes in ‘Yanli’ (Supplementary Table S4), 1128 unigenes in ‘Benihoppe’ (Supplementary Table S5), and 226 unigenes that intersected.

The expression patterns of these unigenes were analyzed in ‘Yanli’ and ‘Benihoppe’ during infection with C. gloeosporioides according to their FPKM values. The unigenes responsive to C. gloeosporioides infection in ‘Yanli’ but not responsive in ‘Benihoppe’ are shown in Figure 4A; those responsive in ‘Benihoppe’ but not ‘Yanli’ are shown in Figure 4C. To functionally classify these unigenes, we observed significant enrichment of these genes in multiple GO categories. The top GO terms with unigenes in ‘Yanli’ included cell wall organization process, flavonoid biosynthetic process, cutin biosynthetic process, response to biotic stimulus process and plant-type secondary cell wall biogenesis. The top KEGG pathways with unigenes in ‘Yanli’ included Flavonoid biosynthesis, Phenylpropanoid biosynthesis, Gap junction, Plant–pathogen interaction and Pathogenic Escherichia coli infection (Figure 4B). In ‘Benihoppe,’ 1128 unigenes were enriched in GO terms including response to heat, defense response to other organism, response to stress, and response to cold. The top KEGG pathways with unigenes in ‘Benihoppe’ included Protein processing in endoplasmic reticulum, MAPK signaling pathway, and Estrogen signaling pathway (Figure 4D).

To confirm the RNA-seq results, some genes were validated using qRT-PCR. In the plant–pathogen interaction pathways, nine strawberry genes (CERK1, EDS1, RPM1, CNGCs, CDPK, Rboh, CaMCML, WRKY33 and PR1) were up-regulated. The gene expression patterns (nine DEGs in the plant–pathogen interaction pathways) determined by qRT-PCR showed trends similar to those found in the RNA-seq data, which suggested that our transcriptome analysis was reliable (Figure 5).