Four-year-old Chinese chestnut (Castanea mollissima Blume) plants of two varieties, namely, “HL” and “SW,” were planted in Chinese chestnut germplasm resources unit, Nanjing, Jiangsu Province, China, under same growth conditions (i.e., temperature, humidity, and light). The variety “HL” was susceptible to GWDK infestation, while “SW” was partially resistant to GWDK infestation as manifested by the appearance of less and smaller galls on the resistant plants following the methodology by Acquadro et al. (2020). Buds with galls were harvested from 12 plants at 3 time points starting from budburst from April 7 until April 26 [Initiation stage (A) at April 7, growth (B) at April 15, and maturation (C) at April 26] in 2019 (Figure 1). Each sample contained three biological replicates. Bud samples of chestnut were harvested from 12 plants (four plants per replicate) and used for the enzyme activity test, RNA extraction, and RNA-sequencing. The fresh buds were sampled and frozen in liquid nitrogen and stored at –80°C until further use.

The H₂O₂ content (n mol g^–1 fresh weight) from chestnut buds was assessed using the xylenol orange method, as described by Loreto and Velikova (2001). The absorbance was measured using spectrophotometer at 560 nm, and the experiment was repeat at least 3 times.

All the antioxidant enzymes described below were extracted by homogenizing 0.3 g of frozen buds in phosphate buffer solution (pH 7.4 and 4% polyvinylpolypyrrolidone) and then centrifugated at 16,000 × g for 30 min at 4°C; the supernatants were collected and used for antioxidant enzyme activities assays. The activities of SOD, CAT, PPO, and POD were assayed according to the methodologies proposed by Beers and Sizer (1952), Giannopolitis and Ries (1977), Wissemann and Lee (1980) and Amako et al. (1994), respectively, with some modifications. The activities of CAT, POD, PPO, and SOD were determined from the measurement of absorbance at 240, 420, 410, and 560 nm, respectively. The results of aforementioned enzymes were expressed as U mg^–1 FW; three biological replicates were utilized to make each of the measurements.

A total of 100 mg frozen bud samples for 2 species of 3 development stages were ground in a mortar under liquid nitrogen, and the resultant powder was used for RNA extraction with the TRIzol^® extraction reagent (Invitrogen, United States.). Nano Drop and Agilent 2100 bioanalyzer (Agilent Technologies, United States) were used to assess the quality of the RNA sample, and high-quality RNA samples with an A260/A280 ratio of 1.8–2.1 were used as template material for RNA-Seq and quantitative reverse transcription PCR (qRT-PCR) analyses. In total, 18 libraries were prepared (i.e., 3 biological replicates per sample) and sequenced using the SE 50 using Illumina Hiseq2500 platform (BGI, Shenzhen, China); the raw reads were filtered by removal of adapter and low-quality sequences using SOAPnuke v1.5.2. Data are available under NCBI SRA accession PRJNA791965. The clean reads were mapped to Chinese chestnut genome¹ using HISAT2 (v2.0.4); fragments per kilobase per million reads (FPKM) values were applied to calculate the gene expression levels. The DESeq (V1.4.5) package was used to identify the differentially expressed genes (DEGs); false discovery rate (FDR) value ≤ 0.05 and |log2Fold Change| ≥ 1 were set for significant gene expression differences between two samples. To gain insights into the changes in phenotype, the enrichment analysis of the DEGs was performed with gene ontology (GO)² and Kyoto Encyclopedia of Genes and Genomes (KEGG³) annotation. Annotations of some DEGs also referred to the genome of hardwood⁴. The significant levels of the terms and pathways were corrected using a rigorous threshold (q-Value ≤ 0.05) via Bonferroni correction.

Expression pattern verification of 8 randomly selected DEGs from RNA-Seq was performed using qRT-PCR, as described previously (Zhu et al., 2019). Complementary DNA (cDNA) was transcribed using the Takara Prime Script TM-RT PCR reagent Kit (Takara, Japan) according to the manufacturer’s instructions. Specific primers (Supplementary Table 1) for qRT-PCR were designed using the Primer Premier 5 software. Actin gene was selected as internal reference in this study. qRT-PCR was conducted using the ABI7500 RT-PCR system according to the manufacturer’s instructions (Applied Biosystems, Foster City, CA, United States). All reactions and non-template controls were performed in triplicate. Relative transcription levels were calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Data are expressed as mean ± standard error (SE) (n = 3).

Candidate genes associated with the GWDK infestation were identified based on the correlation of transcript levels with antioxidant enzymes activities during berry development among the 2 varieties. Weighted gene coexpression network analysis (WGCNA) was conducted to get the gene coexpression network based on total 29,886 DEGs, Visualization of the gene coexpression network was performed using Cytoscape (version 3.8.2) (Vasaikar et al., 2018). Correlation network was performed using the OmicStudio tools⁵.

Figures were made using GraphPad prism 8.0 (San Diego, CA, United States.). All the analyzed data were expressed as means ± standard error (SE). SPSS Version 17.0 (SPSS Inc. Chicago, IL, U.S.A.) and Excel software were used to perform statistical analyses. A one-way ANOVA with a Duncan’s t-test was used to evaluate the significant differences (p ≤ 0.05). TBtools (v 0.665) software were applied for heatmaps.

The basic physiological indicators of the different development stages of gall formation by D. kuriphilus are shown in Figure 1. For “HL,” the galls are smaller at the initiation stage (A), while during the bud development, the galls become larger and more visible (B). At the maturation stage (C), galls in the buds showed a darker red appearance, while in the buds of “SW,” galls were more difficult to spot during all the development stages, and the average number of galls in the “HL” was significantly higher than “SW” (7.5 times higher).

Hydrogen peroxide content is an important indicator of oxidative damage. Therefore, the H₂O₂ content was investigated during the 3 development stages of “HL” and “SW” buds. At the initiation stage A, the H₂O₂ content in buds of the 2 experiment species was similar, while in the growth and maturation stage (B and C), increasing H₂O₂ was observed compared with that in the initiation stage. The maximum content of H₂O₂ (42.7 n mol g^–1 FW) was observed in the “SW” buds at the growth stage (B), and the content of H₂O₂ was significant higher in “SW” buds at the growth and maturation stages compared with “HL” (Figure 2A).

The antioxidant enzymes system, which includes CAT and POD, PPO, and SOD, play significant roles in preventing or alleviating the herbivory damage of plants resulting from reactive oxygen. These enzyme activities were therefore measured in “HL” and “SW” buds to evaluate the relative antioxidant levels. As shown in Figures 2B–E, at the initiation stage, PPO and CAT activity was higher in “HL” buds, while in contrast, higher POD and SOD activities were observed in “SW” buds. At the growth and maturation stages (B and C), the activity of CAT, POD, and SOD was generally higher in “SW” buds compared with that in “HL.” The only exception was the CAT activity at stage C, which was higher in HL. In contrast to CAT, POD, and SOD, a higher PPO activity was observed in “HL” buds during all three stages, which was statistically significant for stages A and C. Overall, in the majority of cases, the four enzyme activities were lowest at the initially time (i.e., stage A) and increased rapidly toward stage B and remained stable or increase slightly for rest of the examination period. The higher values of H₂O₂ content and antioxidant enzyme activities were generally observed in the resistant “SW” buds, as compared with “HL” buds, apart from the PPO activity. The values of these parameters almost showed an increasing to decreasing trend during the development of both “HL” and “SW” buds.

The samples of the two chestnut varieties at three developmental stages (i.e., initiation, growth, and maturation) were subjected to RNA-Seq transcriptomics analysis. High-throughput RNA-Seq generated between 1.19 and 1.21 Gb of clean base reads from each library (Table 1). After stringent filtering of low-quality sequences, on average of 23.87 M clean reads were obtained with quality scores > Q30, which represented > 91.12% of the reads that were subsequently mapped to the reference genome (see text footnote 1). The aligned ratios of the reads ranged from 71.01% to 74.96% among the eighteen libraries, and 64.95–71.97% of the reads were uniquely mapped (Table 1). These data represented a high sequencing depth and quality sufficient for further transcriptomics analysis.

The expression levels of all the sequenced samples were then calculated via the FPKM method and subjected to principle component analysis (PCA) and Pearson correlation analysis (Figure 3A). All samples showed a degree of correlation (the within-group correlation was at least R² > 0.962). PCA of the expression profiles of the 18 libraries [i.e., 6 samples and 3 biological replicates (Figure 3B)] revealed that the samples from the three different stages and 2 cultivars could be clearly separated, confirming that between-stage variation was relatively high and that the three stages showed distinct global expression patterns, which were suitable for further analysis. PCA data (PCA1 variance explained 66.57% and PCA2 17.79%) represent 84.36% of the variance across the two dimensions. Overall samples HL and SW stages B clustered most closely together and both clustered closely with HL stage C, while the other combinations were more distant.

After the estimating of the gene expression levels through the FPKM method, 12,793 protein-coding DEGs were identified (log2 (fold change) ≥ 1 and FDR < 0.05) and further analyzed as the critical genes that associate with GWDK infestation (Figure 4A). Among them, 4,975, 3,017, and 8,293 DEGs were differentially expressed (log₂ (fold change) ≥ 1 and FDR < 0.05) at stages A, B, and C, respectively, when comparing the two cultivars (Figure 4B).

To further investigate the functions of the detected 12,793 DEGs [log₂ (fold change) ≥ 1], GO-based enrichment was performed with a threshold value of p-Value < 0.05 to evaluate significantly enriched GO pathways. A total of 6,674 DEGs (52.16%) were annotated as “biological process,” which contained the majority of GO terms, followed by 4,820 DEGs (37.15%) annotated as “cellular component” and 5,162 DEGs (39.79%) that were annotated as “molecular function.” Through GO assignments, the assigned DEGs were divided into 41 functional groups. The major subcategories along with the analysis of all the transcripts among the three different stages are shown in Figure 5A. Transcripts associated with small molecule binding (GO: 0044283) accounted for the highest number of enriched GO terms among stage A comparisons, whereas catalytic activity (GO: 0003824) include the highest number of DEGs among both stages B and C classification (Figures 5B,C).

Various DEGs are coordinately expressed to perform their different biological functions; therefore, the KEGG pathway-based analysis of DEGs was performed to identify the candidate pathways and genes associated with GWDK infestation. Based on the KEGG analysis, 12,793 DEGs were allocated to 136 pathways with the most enriched pathways as shown in Table 2 (q-Value < 0.05, top 10 numbers of DEGs). In the comparisons of stages A, B, and C, plant hormone signal transduction (ko04075) and MAPK signaling pathway (ko04016) were detected as the top 2 significant KEGG pathways (q-Value < 0.0053) at all three stages and included the highest number of DEGs. Furthermore, in the comparison at stages B and C, the peroxisome pathway (ko04146) was also significantly enriched (q-Value of 0.049 and 1.28E-05, respectively). GO and KEGG annotations for all the DEGs are shown in Supplementary Table 2.

In this study, 157 and 178 transcripts (log₂ (fold change) ≥ 1 and q-Value < 0.05) linked to plant hormone signal transduction (Supplementary Table 3) and MAPK signaling (Supplementary Table 4) were identified by comparing the three development stages of two Chinese chestnut species, respectively. In the hormone signal transduction pathway, the ABA, auxin, brassinosteroid, cytokinin, ethylene, gibberellin, and jasmonic acid (JA) signaling pathways were further analyzed (Supplementary Table 3). Among them, the JA signaling pathways genes were commonly reported and are often associated with plant herbivory defense. In this pathway, 19 DEGs were identified, which were annotated as coronatine-insensitive, TIFY or MYC, and bHLH transcription factors. Across the three stage comparisons, 5 TIFY protein (i.e., EVM0015448, EVM0017479, EVM0018121, EVM0013408, and EVM0017717) were all higher expressed in SW compared with HL (Supplementary Table 3). In the MAPK signaling pathway, the majority of DEGs were annotated as leucine-rich repeat (LRR) receptor-like kinase/receptor like (85 out of 178) and members of the transcription factor family WRKY (17 out of 178) and bHLH [9 out of 178 (Supplementary Table 4)]. Among the 178 DEGs of the MAPK signaling pathway, EVM0018807 and EVM0025267 were detected as the top two highest most highly expressed genes at the maturity stages A and B, respectively, which were annotated as acanthoscurrin-1-like and DNA damage-repair/toleration protein DRT100-like, respectively. EVM0025267 showed consistently higher expression in the HL variety at all three stages, while EVM0018807 was most downregulated at stage B in the HL (33-fold) variety vs. most downregulated at stage C in SW (350-fold).

Transcription factors can play a major role in regulating gene expression of target genes by binding to their promoter regulatory elements and thus affecting transcriptional activity of the targets. A total of 244 DEGs encoding TFs belonging to 20 major transcription factor families were identified across the three stages (Supplementary Table 5). Members of the bHLH (42/244), WRKY (40/244), NAC (30/244), and MYB families (21/244) represented the largest number of DEG transcription factors (Supplementary Table 5). To identify the potential functions of TFs that are associated with antioxidant enzyme-related genes, a correlation analysis of gene expressions between these four TF families and POD-related genes (POD, Supplementary Table 6) was performed (Figure 6 and Supplementary Table 7). CmbHLH130 (EVM0032437) was positively correlated with the expression of POD pathway genes EVM0033037 (r = 0.98), EVM0018305 (r = 0.95), and EVM0010817 (r = 0.92), while CmIRL3 (EVM0030563) was negatively correlated with the expression of EVM0027125 (r = –0.86), EVM0026612 (r = –0.82), EVM0021463 (–0.91), EVM0010016 (–0.93), and EVM0005622 (–0.95). Expression of the POD gene EVM0033037 positive correlated with 20 members of WRKY TFs family (Figure 6D and Supplementary Table 7) (20 out of 40), while the expression of POD genes EVM0018305, EVM0011684, and EVM0010817 correlated with 19, 9, and 18 WRKY TFs, respectively, highlighting that WRKY TF members mostly positively correlated with these 4 POD-related genes. In contrast, only POD gene EVM0007169 has shown a negative correlation with the expression of a large suite of WRKY family genes (12 out of 40). This trend was also apparent in the correlation network of NAC with POD, where EVM0018305 and EVM0010817 (both 13 out of 30) were positive correlated with NAC expression (Figure 6C and Supplementary Table 7). For the MYB family, the majority of members showed a positive correlation with POD-related gene expression (Figure 6B and Supplementary Table 8) with about 70% positive correlations. For example, EVM0018305 only positively correlated with seven MYBs (7/21), while in contrast, EVM0007169 only negatively correlated (5/21).

Based on the transcriptomics data and antioxidant enzyme activities data, we constructed a set of correlation networks between the five main antioxidant enzyme activities (i.e., PPO, CAT, POD, SOD, and H₂O₂) and the related MAPK signaling pathway gene expressions of 6 tested samples (stages A, B, and C of “HL” and “SW”). As shown in Figure 7 and Supplementary Table 8, most of MAPK signaling pathway members showed a positive correlation with PPO and CAT activities, while the majority was negatively associated with POD, SOD, and H₂O₂ activity/content. POD was significantly correlated with the MAPK-related genes basic leucine zipper 61 (EVM0029813, r = 0.84, p-Value < 0.05), receptor-like protein 12 (EVM0032794, r = 0.89, p-Value < 0.05), and the transcription factor MUTE (EVM0020177, r = 0.87, p-Value < 0.05). In contrast, PPO activity was positively correlated to several gene members of the MAPK signaling pathway.

To identify the antioxidant enzyme genes that may correlated to antioxidant enzyme activities, we constructed a set of correlation networks among the five main antioxidant enzyme activities (i.e., PPO, CAT, POD, SOD, and H₂O₂) to evaluate the expression of antioxidant enzyme-related genes. The results revealed that PPO activity was positively correlated with EVM0026323 (r = 0.85, p-Value < 0.05) and EVM0006016 (r = 0.72, p-Value < 0.05); CAT activity was positively correlated with EVM002624729 (r = 0.81, p-Value < 0.05); POD activity was positively correlated with EVM0026612 (r = 0.89, p-Value < 0.05), EVM0005622 (r = 0.84, p-Value < 0.05), and negatively correlated with EVM0011684 (r = –0.71, p-Value < 0.05). SOD activity was positively correlated with EVM0032156 (r = 0.9, p-Value < 0.05), EVM0033049 (r = 0.79, p-Value < 0.05), while EVM0005687 expression was negatively correlated (r = –0.79, p-Value < 0.05) with buds H₂O₂ contents (Supplementary Table 9). All the aforementioned antioxidant enzyme activity-related genes could be important candidates for further Chinese chestnuts “antioxidant enzyme” study.

Terpenes are a class of secondary metabolite volatile compounds derived from the mevalonate and MEP substrate pathways that are often associated with direct and indirect defenses to pests and pathogens. In this study, GWDK feeding exhibits a significant effect on the elicited expression of terpene synthase (TPS) genes (Supplementary Table 10). For example, the expression of TPS gene EVM0000042 was annotated as (E,E)-alpha-farnesene synthase and presented an 8-fold higher expression in HL than SW at the initiation stage A (Supplementary Table 10), while two (3S,6E)-nerolidol synthase genes EVM0024276 and EVM0000019 were induced 3- to 12-fold in both cultivars after infestation (Supplementary Table 10). Interestingly, (α)-farnesene and (E)-nerolidol have both been reported to be associated with herbivory defenses in plants (Jing et al., 2021).

To validate the expression patterns for the DEGs generated by RNA-Seq, the expression levels of 8 randomly selected genes were further analyzed by qRT-PCR across the 6 sequence samples (i.e., 3 biological repeats per samples). The qRT-PCR data showed very similar patterns and trends of expression compared with the RNA-Seq transcriptomics data (Figure 8 lines vs. bar graphs). These results validate that the RNA-Seq data can be considered reliable and reproducible by an independent method (i.e., qRT-PCR) using a random selection of genes (n = 8).