The hidden Markov model profile of PP2C (PF00481) was used to search for PP2C genes in the N11-1 genome of Chinese chestnut (Wang J. et al., 2020). The results were validated by Pfam (http://pfam.xfam.org/search) and NCBI Batch CD-Search with the CDD database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The CmPP2C amino acid sequences were analyzed on ExPASy (https://www.expasy.org) for length, isoelectric point, and relative molecular weight.

We performed multiple PP2C full-length sequence alignments and constructed phylogenetic tree using MEGA X (Kumar et al., 2018). Clustal W was used to create multiple sequence alignments. A neighbor-joining phylogenetic tree was constructed using p-distance substitution model and partial deletion gaps data treatment, and node support was estimated by conducting 1000 bootstrap replicates.

Conserved motifs in CmPP2Cs were analyzed using the MEME web server (http://meme-suite.org/) with maximum motif number set to 20. The exon-intron organization of CmPP2Cs was analyzed using general feature format (GFF3) files and visualized on the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/). The cis-acting elements in gene promoters were identified by PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) using 1,500 bp upstream sequence of the transcription start site of each gene.

We constructed the chromosomal localization map of CmPP2Cs using Mapchart 2.32 software (Voorrips, 2002). The syntenic gene pairs within the chestnut genome were identified by MCScanX, and displayed using Circos (version 0.69-8) software. The substitution rates of nonsynonymous (Ka) and synonymous (Ks) were calculated using the KaKs-calculator (version 2.0) with default genetic code table (Standard Code) and default method for estimating Ka and Ks and theirs references (Model Averaging on a set of candidate models) (Wang et al., 2010). The analysis of synteny between the genomes was performed by the Python version MCscan in JCVI utility libraries v1.0.5 (Tang et al., 2024).

The plant materials used in this study were one-year-old seedlings of ‘Yanbao’ Chinese chestnut (a widely planted cultivar), which were planted in the greenhouse (temperature: 24 ± 2 °C; air humidity: 75 ± 5%; day-night rhythm: 14h light/10h dark) of the Hebei Normal University of Science and Technology. All seedlings were planted in pots 40 cm high and 38 cm in diameter, with a mixed substrate of grass charcoal:perlite:vermiculite ratio of 3:1:1. These seedlings were subjected to low temperature (cold, CD), drought (DT), waterlogging (WL), and exogenous ABA, respectively. The CD treatment involved putting the seedlings in a refrigerator at 0 °C for 1 hour. The DT seedlings were cultivated in soil with a 39% moisture content for 22 days. However, the WL seedlings were planted under the condition of soil moisture content of 100% for one hour. The ABA treatment involved spraying seedlings with exogenous ABA (150mg/L) for 3 days. Additionally, the control group seedlings were cultured at 24 °C, with a 61% soil moisture content. All treatments had three replications. Leaves of the treated seedlings were collected to investigate their expression patterns in response to the respective abiotic stresses. All samples were frozen in liquid nitrogen and stored at -80 °C for total RNA extraction.

Total RNA was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China). High-throughput sequencing was performed using the MGI platform, with PE150 reads length, by Annoroad Gene Technology Co., Ltd (Beijing, China). All clean reads were mapped to the N11-1 Chinese chestnut genome using the TopHat v2.1.1 software (Kim et al., 2013), and the number of reads mapped to each gene was counted using the HTSeq v0.11.3 software (Anders et al., 2015). The values of fragments per kilobase of the exon model per million mapped fragments (FPKM) were obtained through a perl script. Then, the differentially expressed genes (DEGs) were counted using DESeq2 with |log₂(fold change)| ≥ 1 and FDR < 0.05 (Anders and Huber, 2010). The gene ontology (GO) enrichment analysis of DEGs was performed using TBtools software. The weighted correlation network analysis (WGCNA) was performed with all expressed genes (FPKM > 1) using the R package (Langfelder and Horvath, 2008), and the co-expression networks were generated using Cytoscape (Otasek et al., 2019).

The full-length coding sequence (CDS) of the CmPP2C31 with its termination codon was amplified from the C. mollissima cultivar ‘Yanbao’, and cloned into plasmid pBWA(V)HS to produce the 35S::CmPP2C31 vector. Agrobacterium tumefaciens (EHA105) with 35S::CmPP2C31 vector was transformed into ‘Yanbao’ Chinese chestnut seeds following a previously published protocol (Zhao et al., 2017). The magnetic nanoparticles (MNP) and plasmids (DNA) were mixed in a 1:1 mass ratio to obtain the MNP-DNA complex. Then, the MNP-NDA complex was mixed with Chinese chestnut pollen and incubated under a magnetic field. Finally, the pollen was given to the female flowers of ‘Yanbao’ to obtain transgenic Chinese chestnut seeds. The harvested Chinese chestnut seeds were cultivated into seedlings on the WPM medium containing 1.0 mg/L trans-Zeatin and 0.1 mg/L IAA at 25 °C for 30 days. After RT-qPCR verification, positive transgenic seedlings were propagated on a medium containing 1.0 mg/L trans-Zeatin, 0.2 mg/L brassinosteroid, 0.01 mg/L IAA, and 1.0 mg/L 6-benzylamino purine. Finally, PEG (15%) was added to the subculture medium to simulate drought, and the control was grown on a culture medium with no additions. Two weeks later, the contents of H₂O₂ (hydrogen peroxide) and MDA (malondialdehyde) were assayed following a previously described method (Sun et al., 2018).

The CDS of CmPP2C31 was cloned into plasmid 1300-GFP for generating the GFP-fused protein in living cells. A. tumefaciens (GV3101) containing the 1300-GFP-CmPP2C31 vector was used to inject the lower epidermis of tobacco (Nicotiana benthamiana) leaves. Subcellular localization of CmPP2C31 was investigated at 72 h after infiltration. The fluorescence signal was detected using a laser copolymerization cross-fluorescence microscope at 488nm and 405nm excitation intensities for GFP and CFP, respectively.

The cDNA of EVM0007407 was cloned by reverse transcription PCR of ‘Yanbao’ leaves. The cDNA was fused into the pGADT7 vector to construct the pGADT7-EVM0007407 recombinant plasmid. The CmPP2C31 promoter (-2000bp) was cloned from ‘Yanbao’. Then, this promoter fragment was inserted into the pAbAi vector to construct the pCmPP2C31-AbAi recombinant plasmid. Y1H Gold yeast cells were co-transformed with pGADT7-EVM0007407 and pCmPP2C31-AbAi plasmids, and clones containing recombinant plasmids from SD/–Leu/–Ura medium were selected and grown in SD/–Leu/–Ura media containing different Aureobasidin A (AbA) concentrations (0ng/ml, 400ng/ml, 500ng/ml) to detected the interaction between EVM0007407 and CmPP2C31-pro (Yang et al., 2019).

We identified 68 genes from the whole genome of N11-1, a seedling Chinese chestnut cultivar (Wang J. et al., 2020). The 68 genes encode putative PP2C family proteins named CmPP2C01 - CmPP2C68 based on their locus on the chromosomes (Chr) ( Supplementary Figure S1 ; Supplementary Table S1 ). The putative PP2C proteins ranged from 186 to 1079 amino acids (aa) length, with 20.31 kDa to 119.40 kDa predicted molecular weights. The theoretical isoelectric point (pI) of the CmPP2C proteins was 4.66 to 9.30 ( Supplementary Table S1 ). Gene ontology (GO) analysis revealed that these CmPP2Cs enriched ‘protein dephosphorylation’, ‘response to abscisic acid’, ‘response to water deprivation’, and ‘response to oxygen-containing compound’ ( Supplementary Figure S2 ; Supplementary Table S2 ). Therefore, CmPP2Cs may regulate stress resistance in Chinese chestnut by regulating protein phosphorylation modifications.

To investigate the classification and evolutionary relationships of CmPP2C proteins, we constructed an unrooted phylogenetic tree based on the alignments of full-length protein sequences from Chinese chestnut, Arabidopsis, walnut, and rice. The CmPP2Cs were divided into 14 subgroups (A–N) according to their orthologs in Arabidopsis ( Figure 1 ) (Xue et al., 2008). The distribution of PP2C from rice and walnut in different subgroups indicated that different types of PP2C maintained similar functions in species evolution. There were the most (11) CmPP2Cs in subgroup A ( Figure 1 ; Supplementary Table S4 ). These members should have similar functions and might have evolved with the expansion of gene families caused by genome replication events. However, there was only one member in the N subgroup, CmPP2C36, which may have an independent evolutionary trajectory from other members. Further, the ten conserved motifs of CmPP2C proteins were identified, and showed clade-specificity. Motifs 1, 2, 3, 5, and 8 were widely distributed among CmPP2C proteins and formed two combinations (1-8-3-2 and 5-8-3-2), the possible core domains of the PP2C family ( Supplementary Figures S3 , S4 ; Supplementary Tables S3 , S4 ). Similarly, CmPP2C genes from same subgroup had similar exon-intron structure and differ between subgroups ( Supplementary Figure S3C ; Supplementary Table S4 ). These indicated that the evolution and divergence of CmPP2Cs might have occurred at an early stage.

Whole-genome duplication (WGD) analysis revealed synteny relationships among 15 pairs of CmPP2Cs genes across nine chromosomes ( Figure 2A ). In particular, CmPP2C64 and CmPP2C65 located in 8,039 bp apart of Chr12, indicating a likely tandem duplication event ( Supplementary Figure S1 ; Table 1 ). These observations suggested that both segmental and tandem duplication possibly expanded the CmPP2C gene family. Furthermore, the values of Ka, Ks, and Ka/Ks ratios were calculated to estimate the dates of duplication events of CmPP2C gene pairs. The divergence time varied from 6.18 to 67.86 million years ago (MYa), spanning the Neogene (2.6-23.3 MYa) and Paleogene (23.3-68.5 MYa) periods. Five pairs of duplicated CmPP2Cs might have been positively selected at the most recent divergence (6.18-16.24 MYa (Ka/Ks > 1)). Conversely, ten pairs might have undergone purifying selection in the early stages (40.14-67.86 MYa). The Ka/Ks ratios of CmPP2C59/66 and CmPP2C06/14 were both 0.09, indicating that they underwent the highest selection pressure to adapt to the natural environment in Paleogene ( Figure 2B ; Table 1 ). Therefore, we could speculate that natural environmental pressure (or abiotic stress) promoted the duplication of the CmPP2Cs.

Five comparative syntenic maps of Chinese chestnut association with Quercus robur (oak), Juglans regia (walnut), Malus domestica (apple), Arabidopsis thaliana, and Oryza sativa (rice) were constructed to reveal the evolutionary relationships of the PP2C gene family among different species ( Figure 2C ; Supplementary Figure S5 ). We found that 51 CmPP2Cs with corresponding orthologous genes in the other five genomes. Forty CmPP2Cs (especially CmPP2C31) were syntenic in at least two genomes, indicating their crucial role in the PP2C family evolution ( Supplementary Table S5 ). Furthermore, 22 CmPP2Cs showed syntenic relationships with oak PP2C genes, 30 with walnut, 44 with apple, and 17 with Arabidopsis. Chinese chestnut and rice (the only monocot) shared eight CmPP2C gene pairs, which was much less than between Chinese chestnut and the four dicots. These results may suggest that most orthologous pairs occurred after the divergence of dicotyledons and monocotyledons. The phylogenetic tree of the PP2C family orthologous genes ( Figure 2C ; Supplementary Table S6 ) supported this inference.

The transcriptome of chestnut seedling leaves from CD, DT, WL, ABA, and control treatments (with three replications) revealed the roles of CmPP2Cs in stress resistance and signal transduction ( Figure 3A ; Supplementary Figure S6 ). RNA-seq generated 5.16-7.07 gigabases (Gb) clean data from each of the 15 libraries ( Supplementary Table S7 ). The Chinese chestnut N11-1 genome contained 33,597 annotated genes, and 14,116 genes were almost non-expressed (FPKM < 0.1) across the 15 samples, which included eight CmPP2Cs ( Supplementary Table S8 ). High Pearson correlation coefficients (r > 0.74) indicated high-quality control among the biological replicates ( Supplementary Table S9 ). The hierarchical clustering tree divided CmPP2C genes into five subclades based on their expression profiles. Subclade I had the greatest number (17) of CmPP2Cs, but the CmPP2Cs had the lowest expression (FPKM < 1) across treatments. In contrast, the 14 members of subclade III had the highest FPKM values (10.80 ≤ FPKM ≤ 86.80). The CmPP2Cs in subclade V also showed relatively high expression levels (5.41 ≤ FPKM ≤ 79.63). Interestingly, CmPP2C19 and CmPP2C31 expression were the most sensitive to different treatments (Fold change > 3 or < 0.33) in subclade V ( Figure 3B ; Supplementary Tables S8 , S10 ). Additionally, the relative expression levels measured by RT-qPCR confirmed the FPKM patterns of five CmPP2Cs, with high correlation relationship (r > 0.86) ( Supplementary Figure S7 ).

Furthermore, a total of 6,240 DEGs were identified in the CD, DT, WL, and ABA treatments by comparing with the control ( Supplementary Figure S8 ; Supplementary Table S8 ). Twenty-two CmPP2Cs were included in DEGs. Seven CmPP2Cs were significantly up-regulated, and 14 were significantly down-regulated in one or more treatments. Interestingly, ABA significantly up-regulated CmPP2C31, and CD and DT significantly down-regulated it ( Figure 3B ; Supplementary Tables S8, S10 ).

The WGCNA analysis using the 14,184 genes (FPKM > 1) identified 15 co-expression modules, with gene numbers ranging from 49 (MEmidnightblue) to 3,346 (MEturquoise) ( Supplementary Figure S9 ; Supplementary Table S8 ). Five key modules (MEblue, MEgreen, MEred, MEmidnightblue, and MEyellow) were significantly associated with stress (|r| ≥ 0.94 and p < 0.05). MEmidnightblue was positively related to exogenous ABA treatment (r = 0.97), and MEyellow showed the opposite results (r = -0.96). MEblue, MEgreen, and MEred modules were positively related to CD (r = 0.97), DT (r = 0.96), and WL stresses (r = 0.94), respectively ( Figure 4A ). Furthermore, the module eigengene-based connectivity (KME, |KME| > 0.85) and the weight value (> 0.25) between pairwise genes revealed hub CmPP2Cs in five key modules. MEblue had three hub PP2Cs (CmPP2C38, CmPP2C42, and CmPP2C68) and MEyellow had one (CmPP2C31). MEred, MEgreen, and MEmidnightblue no hub CmPP2C gene. All selected CmPP2Cs were differentially expressed ( Figure 4B ; Supplementary Table S8 ).

The promoters of CmPP2C31, CmPP2C38, CmPP2C42, and CmPP2C68 contained at least one stress-related cis-element, including ACE, ABRE, G-box, MYC, and MYB ( Supplementary Table S11 ). Further, MEyellow and MEblue contained three (EVM0007407, EVM0033064 and EVM0028466) and six (EVM0002559, EVM0005863, EVM0009459, EVM0012518, EVM0030390 and EVM0032542) transcription factors (TF), respectively. These TFs could bind to the cis-acting elements and were significantly co-expressed with CmPP2C31, CmPP2C38, CmPP2C42, and CmPP2C68 ( Figure 4B ; Supplementary Figures S10 , S11 ). Especially, EVM0007407 (NAC072 gene) and EVM0032542 (NAC083 gene) showed the highest expression in MEyellow and MEblue and were co-expressed with CmPP2C31 and CmPP2C42/68, respectively ( Figure 4C ; Supplementary Figure S7 ). These findings suggested that the hub CmPP2Cs in the key modules, regulated by their upstream TFs, may participate in the stress response pathways of Chinese chestnut seedlings.

CmPP2C31 was the only up-regulated PP2C in ABA, the only down-regulated in CD/DT treatments, and the only hub CmPP2C in MEyellow. Thus, transgenic experiment by pollen magnetofection (Zhao et al., 2017) were conducted to explore its function in drought stress resistance. The over-expression vector 35S::CmPP2C31 was generated and then transformed into the seed of ‘Yanbao’ Chinese chestnut. We cultured the embryos of transgenic seeds into seedlings under tissue culture conditions, and 15% PEG was added to the subculture medium to simulate drought stress ( Figure 5A ). The relative expression of CmPP2C31 was significantly increased in CmPP2C31-OE than the WT under normal and drought stress conditions ( Figure 5B ). Drought stress significantly inhibited the growth of Chinese chestnut seedlings and significantly increased the H₂O₂ and MDA content. However, over-expressing CmPP2C31 in Chinese chestnut seedlings significantly reduced these trends ( Figure 5 ). Furthermore, we found that over-expression of CmPP2C31 was able to increase endogenous ABA content through feedback regulation ( Figure 5E ). These results indicated that CmPP2C31 improved the tolerance of Chinese chestnut to drought stress by affecting ABA signaling transduction.

Subcellular PP2C localization and interaction with TF EVM0007407 revealed its regulation mechanism in drought resistance. The full-length CDS of CmPP2C31 was fused to the C-terminus of green fluorescent protein (GFP). The GFP-fusion proteins and GFP alone were transiently expressed in N. benthamiana. The nucleus was identified by fusing the nuclear marker protein GHD7 with the cyan fluorescent protein (CFP). The fluorescent signals of GFP-fusion proteins were exclusively restricted to the nucleus, fully overlapping with the fluorescence of GHD7-CFP ( Figure 6A ), indicating that CmPP2C31 is localized in the nucleus. In contrast, the GFP signal alone was in both the nucleus and cytosol.

The protein of EVM0007407 is highly homologous to AtNAC072 of Arabidopsis and is the TF most likely regulating CmPP2C31 expression under abiotic stress ( Figures 4B, D ). A Y1H identified the potential interaction between EVM0007407 and the promoter of CmPP2C31. Yeast strains containing pGADT7-EVM0007407&pABAi-ProPP2C31 grew normally in media containing 400 and 500 ng/ml Aureobasidin A (AbA) ( Figure 6B ). Therefore, EVM0007407 could directly bind to the promoter of CmPP2C31.