Colletotrichum fructicola strains, PAFQ31 and PAFQ32, were obtained from pear leaves (P. pyrifolia cv. Cuiguan, susceptible to pear anthracnose) that caused TS and BnL symptoms on pear, respectively (Fu et al., 2019). Plant tissue-cultured materials from P. pyrifolia cv. Cuiguan (3 years old) were grown in the greenhouse of the College of Plant Science and Technology of Huazhong Agricultural University. Pathogenicity analysis was performed on the attached leaves of P. pyrifolia cv. Cuiguan. The inoculation protocols have been previously described (Zhang et al., 2015). Briefly, the attached pear leaves (approximately 4-week-old) were selected for uniform size, color, and an absence of visual defects. Each leaf was surface-sterilized with 75% ethanol, washed with sterile water, air-dried, and inoculated with 40 μL of conidial suspensions (1 × 10⁶ conidia/mL) of C. fructicola by spraying the conidial suspensions with a handheld sprayer over the leaves. Control leaves were inoculated in parallel with sterile water. Five leaves were inoculated per biological replicate with three biological replicates per treatment. The inoculated leaves were subsequently incubated at 25°C, about 85% relative humidity (for 1 week) with a 12/12 h light/dark photoperiod for disease development.

For the RNA-sequencing (RNA-Seq) experiments, the samples of pear leaves infected with PAFQ31 and PAFQ32 were collected at the quiescent stage (QS) and necrotrophic stage (NS), respectively, and frozen in liquid nitrogen. Mock inoculation with sterile water at QS was chosen as the pear leaf control. Three independent biological replicates were sequenced for each treatment.

The total RNA was extracted from leaf tissues using PureLink™ Plant RNA Reagent (Cat# 12322012) according to the manufacturer’s instructions (Invitrogen, Carlsbad, CA, United States) and genomic DNA was removed using DNase I (Takara, Shiga, Japan). The RNA quality was determined by a 1% denaturing agarose gel and a NanoDrop 2000 system (Thermo Scientific, Wilmington, DE, United States). RNA integrity was assessed using the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, DE, United States). RNA-Seq transcriptome libraries were prepared following the TruSeq™ RNA sample preparation Kit from Illumina (San Diego, CA, United States) using 1 μg of the total RNA. Library preparation and sequencing were performed by Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China). The libraries were generated using NEB Next^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, MA, United States) and were sequenced using an Illumina HiSeq X Ten platform in 2 × 150 bp pair-end sequencing mode. The RNA-Seq raw data have been deposited in the NCBI Sequencing Read Archive database under the project accession number PRJNA698408. The raw paired-end reads were trimmed, and the quality was controlled by using SeqPrep¹ and Sickle² with default parameters. Then clean reads were aligned to the reference genomes of the pear (P. bretschneideri cv. Suli)³ with orientation mode using HISAT2 (Kim et al., 2019) since the P. pyrifolia genome has not been released when we conducted this research.

To identify DEGs between two different samples, the gene expression levels were calculated according to the fragments per kilobase of transcript sequence per millions of base pairs sequenced (FPKM) method. FeatureCounts⁴ was used to quantify gene abundances. DESeq2⁵ (Love et al., 2014) was utilized for identifying DEGs. A significance analysis (| log₂Fold Change| > 1 and P-value <0.05) of the results was used to identify genes that are strongly upregulated or downregulated by C. fructicola infection using unlogged data. Fold-changes were calculated with average transcript levels compared to control values that were, in turn, log₂-transformed and calculated for Spearman correlation coefficients between treatments. Heatmap analyses were performed using OmicShare tools⁶. Functional-enrichment analyses including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at a corrected P-value ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analyses were carried out by Goatools⁷ and KOBAS⁸. The analyses were based on the online platform of Majorbio Cloud Platform⁹.

Gene co-expression network analysis was performed using the weighted gene co-expression network analysis (WGCNA) package¹⁰ based online platform of Majorbio Cloud Platform. Gene cluster dendrogram was constructed with colors based on the correlations between the expression levels of genes and used to build clustering trees and to divide modules. Besides, the correlation between modules and samples was analyzed using WGCNA. The protein-protein interaction (PPI) networks were constructed for Pyrus based on the data of Arabidopsis thaliana produced by the Search Tool for the Retrieval of Interacting Genes/Proteins database (STRING, https://string-db.org/) and visualized using Cytoscape software (version 3.7.2)¹¹. The nodes of the network represented proteins encoded by DEGs and their functional partners in the predicted pairwise interaction network.

Transcription factor (TF) families were identified using the Plant Transcription Factor Database PlantTFDB 4.0¹² (Jin et al., 2017). The Hmmscan E-value and BLAST E-value are all set to 1e^–05. The TF genes were classified into various plant transcription factor families based on conserved domains predicted in the above analysis.

The complementary DNA (cDNA) synthesis was performed on 1 μg of the total RNA (previous RNA-Seq library construction) with the 5× All-In-One RT MasterMix (with AccuRT Genomic DNA Removal Kit; Applied Biological Materials Inc., United States) according to the instructions of the manufacturer. Samples of cDNA were diluted 1:2 to the final template concentration for quantitative Real-time PCR (qRT-PCR) by using CFX96 Real-Time PCR Detection System (Bio-Rad, United States). PCR amplification was performed with 5 μL of iTaq™ Universal SYBR^® Green Supermix, 3.5 μL of nuclease-free H₂O, 1 μL of cDNA, and 10 pmol of forward and reverse primers in a final volume of 10 μL. The PCR cycles are as follows: 95°C for 20 s, followed by 40 cycles of 95°C for 5 s, 56°C for 15 s, and 72°C for 15 s. At the end of the reaction, melt curves were run for all primer pairs to check for dimerization. The P. pyrifolia beta-adaptin-like gene A (XM_009369647) was used to normalize the RNA samples for each qRT-PCR and relative expression levels were calculated based on the 2^–ΔΔCT method (Schmittgen and Livak, 2008). All qRT-PCR primers are listed in Supplementary Table 1. Each treatment consisted of two biological repeats and three technical replicates.

The contents of JA and SA were detected using previously described methods at Shanghai Applied Protein Technology Co., Ltd (Shanghai, China) (Shao et al., 2019). Briefly, ∼100 mg pear leaf from each treatment was ground to a fine powder in liquid nitrogen using a mortar and pestle. The plant hormones were extracted from the powder at 4°C for 12 h with 1 mL of ethyl acetate supplemented with internal standards and then centrifuged at 14,000 g for 15 min at 4°C. The supernatant was carefully transferred to a new 1.5 mL tube and the pellet was again extracted with 500 μL of ethyl acetate at 4°C for 1 h. The supernatant from the second extraction was collected and pooled with the first extraction. The supernatant was evaporated to dryness under nitrogen gas flow and then dissolved in 500 μL of 50% (v/v) acetonitrile. After being centrifuged (at 14,000 g and 4°C, for 10 min), the supernatant was then analyzed by high-performance liquid chromatography (HPLC)-electrospray ionization-tandem mass spectrometry. The mobile phase consisted of a combination of solvent A [0.05% (v/v) formic acid in water] and solvent B [0.05% (v/v) formic acid in acetonitrile]. The linear gradient was as follows: 2 to 98% B (v/v) for 10 min, 2% B (v/v) for 10.1 min, and hold at 2% B for 13 min. The mass spectrometer (Qtrap 5500 System, AB Sciex) equipped with an electrospray ionization source was operated in positive/negative ionization and multiple reaction monitoring modes. The MS parameters were as follows: source temperature, 500°C; ion source gas 1, 45 psi; ion source gas 2, 45 psi; curtain gas, 30 psi; and ion spray voltage, 5500 V.

Statistical analysis was conducted with IBM SPSS Statistics 21.0 software (IBM Corp. Released 2012. IBM SPSS Statistics for Windows, Version 21.0. Armonk, NY: IBM Corp.) by one-way ANOVA, and means were compared using Duncan’s test at a significance level of P = 0.05. The homogeneity of variance was tested before analysis.

Two pathogenetic strains of C. fructicola identified in our previous study (Fu et al., 2019), PAFQ31 and PAFQ32 causing TS and BnL on pear leaves, respectively, were selected for inoculation on attached pear leaves (P. pyriforia cv. Cuiguan) under the unwounded conditions. As expected, strain PAFQ31 caused TS symptoms at 6 days post-inoculation (dpi), with numbers of the tiny spots rapidly increasing and spreading to the entire leaf area in the following days (Figure 1A2; 9 dpi). At 14–20 dpi, the green parts neighboring the TS started yellowing, subsequently, the veins and petioles turned yellow, and finally, severe defoliation was observed (Figure 1A3). In contrast, strain PAFQ32 did not induce small brown necrotic lesions until 40 dpi, and subsequently, the small lesions expanded into BnL symptoms under appropriate conditions but no defoliation was observed for these diseased leaves (Figure 1A6). In parallel, no lesions were observed on the leaves as inoculated with sterile water. The symptomatic samples were subjected to fungal isolation at the parts neighboring the asymptomatic regions, making colonies matching the inoculated ones according to their morphologies and ITS sequences.

To get insight into the mechanism under the symptoms induced by each strain, the inoculated pear leaves were collected in triplicates for transcriptome sequencing at latent (QS) and evident (NS) stages (Figure 1B), i.e., at 4 and 9 dpi for strain PAFQ31, and 4 and 42 dpi for strain PAFQ32, respectively. Pear leaves without inoculation at QS were involved in parallel as the controls, respectively. In total, 15 RNA samples were prepared for cDNA library construction and sequencing.

In total, 7.0–11.3 Gb of high-quality clean data were obtained from each leaf sample (Supplementary Table 2). The clean reads were aligned to the pear reference genome (P. bretschneideri cv. Suli genome; accession no. PRJNA259338) (Wu et al., 2013), with over 76% reads mapped to the reference genome (Supplementary Table 2). To evaluate the correlations of three independent biological replicates in each treatment and eliminate the possible outliers, principal component analysis (PCA) and hierarchical cluster analysis were carried out for the obtained data and an outlier sample from each treatment was abandoned (Supplementary Figure 1A). Finally, ten samples with high correlations between the biological replicates (R² = 0.969–0.982) were subjected to further analyses (Supplementary Figure 1B).

Based on the DESeq2 analysis, a total of 5,485 (312 in QS and 5,363 in NS) and 7,122 (1,112 in QS and 6,559 in NS) DEGs (no less than twofold change) were revealed in the leaves of P. pyrifolia, with 3,217 (189 in QS and 3,126 in NS) and 3,505 (947 in QS and 2,798 in NS) upregulated in the leaves infected by strain PAFQ31 and PAFQ32 as compared with the controls, respectively (Figure 2A). Venn diagram showed that most of the DEGs were uniquely expressed, and far less overlapped in response to the infection of each strain, with the numbers significantly higher at NS than those at QS (Figure 2B). Of which, only 22 up- and 10 downregulated DEGs were observed in both PAFQ31- and PAFQ32-infected leaves during the infection progress, whereas 3,284 (98 in QS and 3,226 in NS) and 5,182 DEGs (376 in QS and 4,855 in NS) were uniquely expressed in the PAFQ31- and PAFQ32-infected leaves, respectively (Figure 2B). The heatmap analysis further reveals that the DEGs were in a very small amount at QS while dramatically increasing their numbers at NS triggered by each strain, indicating a clear process of P. pyrifolia in response to the fungal infection, from less apparent to dramatically strong along with the infection of C. fructicola (Figures 2C,D).

To characterize the pear response to two C. fructicola strains, we performed KEGG and GO enrichment analyses of the DEGs induced by each strain at QS and NS, respectively. The KEGG enrichment analyses revealed 312 and 5,363 DEGs enriched in 26 and 115 pathways at QS and NS in P. pyrifola leaves as infected by strain PAFQ31, respectively. Of these, only small numbers of DEGs were enriched in plant hormone signal transduction (11 genes; ko04075) and carotenoid biosynthesis (2; ko00906) at QS (Figure 3A), whereas large numbers of genes in plant hormone signal transduction (85; ko04075), plant-pathogen interaction (76; ko04626), phenylpropanoid biosynthesis (50; ko00940), and flavonoid biosynthesis (32; ko00941) metabolic pathways at NS (Figure 3B). Notably, photosynthesis (46; ko00195), carbon fixation in photosynthetic organisms (25; ko00710), and photosynthesis-antenna proteins (15; ko00196) were also significantly enriched at NS (Figure 3B).

As infected by strain PAFQ32, 1,112 and 6,559 DEGs were enriched in 82 and 115 pathways at QS and NS, respectively (Figures 3C,D), with most of the defense-related pathways the same as those of PAFQ31. Besides, stilbenoid, diarylheptanoid, and gingerol biosynthesis (13; ko00945), cutin, suberine, and wax biosynthesis (16; ko00073), porphyrin and chlorophyll metabolism (19; ko00860), and carotenoid biosynthesis (15; ko00906) were also highly enriched at NS (Figure 3D). It is worthy to note that, multiple defense-related pathways were only activated at QS by the infection of PAFQ32 instead of PAFQ31 (Supplementary Table 3). Instead, more of these DEGs at NS are related to the infection by PAFQ31 rather than by PAFQ32. For GO enrichment, the DEGs at NS were mainly enriched in similar metabolic processes for each strain but uniquely enriched in photosynthesis (GO:0015979), photosystem (GO:0009521), photosystem II (GO:0009523), and photosystem II oxygen-evolving complex (GO:0009654) for PAFQ31 infection, suggesting the infection of strain PAFQ31 seriously affect the photosynthesis of plants (Supplementary Figure 2).

The related genes were further analyzed because large DEG numbers were significantly enriched into plant signal transduction pathways. Of which, homologous non-expresser of pathogenesis-related genes (Pp-NPR1 and Pp-NPR2), the key regulators of plant SA-mediated defense response, were mainly upregulated triggered by each strain at QS, whereas the genes at the downstream of SA signaling pathway such as Pp-TGAs (Pp-TGA2.3, Pp-TGA4) and pathogenesis-related gene 1 (Pp-PR1) were upregulated at NS, with higher accumulation levels in the PAFQ31-infected leaves than those in the PAFQ32-infected (Figure 4A). Additionally, the expression of genes involved in the JA, ET, and ABA signaling pathways was specifically upregulated in the PAFQ31-infected leaves (Figures 4B–D). Moreover, the expression levels of genes related to the hormone biosynthesis pathway (the upstream pathway of the hormone signal transduction) were consistent with these related to hormone signal transduction pathways (Figures 4A–D).

The genes related to the auxin signaling pathway, e.g., LAX1 (LIKE-AUXIN1), LAX2, LAX3, IAA2 (auxin-responsive gene), IAA27, and IAA31, were constant at QS while downregulated at NS after the infection by each strain (Figure 4E). In contrast, several genes (LOC103966953, LOC103939002, LOC103934959, LOC103967969, LOC103945994) encoding indole-3-acetic acid (IAA)-induced protein ARG7-like or auxin-induced protein were upregulated at NS after infection by each strain (Figure 4E). Moreover, genes (Gretchen Hagen3.1, GH3.1; LOC103942180, LOC103957346) encoding IAA-amido synthetase GH3.1 were obviously upregulated in the PAFQ31-infected leaves at NS, while downregulated in the PAFQ32-infected (Figure 4E).

To check whether the inoculated leaves behave defensively against the fungal damage, the genes involved in phenylpropanoid metabolism (Figure 5A), an important pathway of secondary metabolism in plants closely related to disease resistance (Wang et al., 2018; Sudheeran et al., 2021), were subjected to an analysis of their expression levels. The related genes, including homologous phenylalanine ammonia lyase (Pp-PAL, LOC103934948) gene, were significantly upregulated at QS or NS after infection by each strain (Figure 5B). Additionally, the genes involved in lignin biosynthesis, including 4-coumarate: coenzyme A ligase (Pp-4CL, LOC103951504), cinnamyl alcohol dehydrogenase (Pp-CAD), p-coumarate 3-hydroxylase (Pp-C3’H), and ferulate 5-hydroxylase (Pp-F5H) gene, were constant at QS but upregulated at NS as compared with the controls (Figure 5B). However, the genes related to anthocyanin synthesis (such as chalcone isomerase/Pp-CHI, LOC103936753; dihydroflavanol 4-reductase/Pp-DFR, LOC103928717; anthocyanidin synthase/Pp-ANS, LOC103952863; anthocyanidin reductase/Pp-ANR, LOC103937289) were significantly upregulated in the PAFQ31-infected leaves (Figure 5B).

To determine the gene regulatory network in pear leaves in response to the C. fructicola infection, a WGCNA was constructed using 7,122 DEGs (after pre-process) in the leaves infected by strains PAFQ31 and PAFQ32. In total, eight gene co-expression modules, containing genes from 48 to 1,691, were identified according to correlations between gene expression levels (Figure 6A). It reveals that each module was correlated with a particular stage (Figure 6B), and some were specifically related to the strains. Of which, the modules ‘MEred’ and ‘MEturquoise’ were highly associated with the infection of strain PAFQ31 at NS, whereas the ‘MEyellow’ and ‘MEbrown’ with strain PAFQ32 at the same stage. KEGG enrichment analyses reveal that the genes in modules ‘MEyellow’ and ‘MEbrown’ were highly enriched in plant-pathogen interaction, plant hormone signal transduction, and phenylpropanoid biosynthesis metabolic pathways (Figure 6C), with the expression patterns as previously indicated (Figures 4, 5 and Supplementary Figures 3A,B). For the modules ‘MEred’ and ‘MEturquoise’, their genes were related to porphyrin and chlorophyll metabolism, carotenoid biosynthesis, circadian rhythm–plant (ko04712), photosynthesis–antenna proteins, photosynthesis, thiamine metabolism (ko00730), sulfur relay system (ko04122), and riboflavin metabolism (ko00740) (Figure 6D). The genes involved in chlorophyll metabolism, carotenoid biosynthesis, and photosynthesis pathways showed different expression patterns between both strains, which were significantly downregulated in the PAFQ31-infected leaves, while upregulated in the PAFQ32-infected (Supplementary Figures 3C–E).

To better understand the phenomenon of early defoliation caused by strain PAFQ31, 225 DEGs of PAFQ31_NS and 243 DEGs of PAFQ32_NS related to plant-pathogen interaction, plant hormone signal transduction, chlorophyll metabolism, carotenoid biosynthesis, photosynthesis, and circadian rhythm pathways were subjected to build PPI networks as referred to the ones in Arabidopsis, respectively (Figure 7). It reveals that most proteins involved in the ET and ABA signaling pathways were upregulated in the PAFQ31-infected leaves, and interacted with the PTI/ETI-related proteins; while the genes involved in the JA and ET signaling pathway were suppressed in the PAFQ32-infected, and the PTI/ETI-related proteins were interactive with the proteins of the SA signaling pathway.

Moreover, we found that several proteins involved in phytohormone signaling, PTI, and ETI are associated with chlorophyll metabolism, carotenoid biosynthesis, or circadian rhythm-related proteins (Figure 7A). Protein phosphatase 2C (PP2C) in ABA signaling interacts with zeaxanthin epoxidase (ZEP) in carotenoid biosynthesis and magnesium-chelatase subunit chlH in chlorophyll metabolism. Moreover, uroporphyrinogen decarboxylase 1 (hemE) in chlorophyll metabolism interacts with different kinds of calmodulin-like proteins (CMLs) and proteins in photosynthesis, likely building a bridge between the plant-pathogen interaction pathway and photosynthesis. By analyzing the expression levels of the corresponding genes of these proteins, it revealed that some of the CMLs interacting with hemE were upregulated and downregulated in the pear leaves after infection by PAFQ31 and PAFQ32, respectively, while corresponding genes related to chlorophyll metabolism pathway and photosynthesis behaved, on the contrary, indicating that these CMLs may serve as key components involved in regulating chlorophyll metabolism and photosynthesis in pear plants and related to the anthracnose-associated early defoliation of pear plants.

To check whether some TFs were involved in the response to the fungal infection, all DEGs were subjected to TF prediction and TF family analysis. A total of 361 and 440 TFs, belonging to 31 and 37 TF families were identified at NS in the leaves infected by strain PAFQ31 and PAFQ32, respectively, and most of the TFs belong to five families, including AP2/ERF, WRKY, bHLH, MYB, and NAC (Supplementary Figure 4A). In total, 34 and 35 TFs belonging to the WRKY family, one of the largest transcription regulator families in plants related to PTI through MAPK signaling cascade regulation and ETI through R protein (Li et al., 2016), were differentially expressed at NS in the leaves infected by strains PAFQ31 and PAFQ32, respectively, with obviously higher expression levels in the PAFQ31-infected leaves (Supplementary Figure 4B). Of these, 21 genes encoding nine WRKY factors (WRKY9, 44, 45, 51, 53, 55, 65, 70, and 75) were significantly upregulated at NS induced by PAFQ31 infection (Supplementary Table 4). It is worth noting that WRKY45, WRKY53, and WRKY75 have been characterized as positive regulators of leaf senescence (Miao et al., 2004; Chen et al., 2017; Guo et al., 2017). Additionally, other TFs, e.g., TF family NAC also known as an immune or leaf senescence regulator, were also differentially expressed at NS in the PAFQ31-infected leaves (Supplementary Figure 4C).

To verify the transcriptome sequencing data, six DEGs were randomly chosen for qRT-PCR analysis. The results showed that all the selected genes exhibited expression patterns similar to those obtained by RNA-Seq, supporting the transcriptome sequencing results (Figure 8A). Correspondingly, gene Pp-PR1 showed significantly higher expression levels at NS in the leaves infected by both strains compared to the healthy leaves, while genes Pp-ERF1 and Pp-PYL were specifically upregulated at NS in the PAFQ31-infected leaves.

To further confirm the phytohormone accumulation triggered by the infection of each strain, endogenous SA and JA were measured in the pear leaves at NS. The results showed that SA accumulation slightly increased (without significance) after infection by both PAFQ31 and PAFQ32 (Figure 8B). JA accumulation was obviously higher in the PAFQ31-infected leaves than in the PAFQ32-infected (Figure 8B).