Mulberry (Morus multicaulis) cutting seedlings were inoculated with phytoplasmas by being grafted with the scions collected from phytoplasma-infected mulberry trees Gai et al. (2018), and the cutting seedlings derived from the same mother tree were grafted with healthy scions and used as controls. The grafted mulberry seedlings were cultivated in a greenhouse under light/dark regime (16 h light/8 h dark) at 24°C with 50–60% humidity. The grafted mulberry seedlings showing Witches' broom disease symptoms were used to detect phytoplasma by PCR amplification of the 16S rRNA gene of phytoplasma using a universal primer set (Forward primer: 5′ TAAAAGACCTAGCAATAGG 3′; Reverse primer: 5′ CAATCCGAACTGAGACTGT 3′) (Ji et al., 2009). The sixth leaves from the top of the grafted healthy and phytoplasma-infected shoots were used for subsequent experiments, and three phytoplasma-infected and healthy plants were used as independent biological replicates for RNA sequencing (RNAseq) and MethylRAD sequencing analyses. A. thaliana (Col-0) and Nicotiana benthamiana seedlings were cultivated in the greenhouse under light/dark regime (16 h light/8 h dark) at 22 and 24°C, respectively, with 50–60% humidity.

Total RNA was extracted using a TRIzol reagent (Invitrogen Corporation, CA, United States), and the poly-A-containing messenger RNAs (mRNAs) were purified using beads with Oligo (dT). The purified mRNAs were used as templates to synthesize complementary DNAs (cDNAs) which were then added to a single “A” base to their ends. These cDNA fragments were subsequently connected with sequencing adapters and then were purified and underwent PCR amplification to create cDNA libraries. Transcriptome sequencing was conducted on the Illumina HiSeq™ 2000S platform, Illumina, San Diego, CA, United States. Raw reads were cleaned and assembled into unigenes and then the assembled unigenes were annotated by BLASTn searches against the M. notabilis genome (https://www.ncbi.nlm.nih.gov/genome/?term=Morus). The expressions of genes were calculated using the reads per kilobase million (RPKM), and if the expression value of a gene changed more than two times (P ≤ 0.05), and the false discovery rate (FDR) was <0.01 between infected and healthy samples, it was designated as a significantly differentially expressed gene (DEG). The identification of DEGs was performed using DESeq software (Zenoni et al., 2010). Three biological replicate experiments for healthy and infected samples and three technical replicates per biological experiment were performed.

Genomic DNA was extracted from the infected and healthy mulberry leaves using the cetyltrimethyl ammonium bromide method (Clarke, 2009) and then was digested with FspEI. The digestion productions were added with adaptor primers and amplified through PCR to create the MethylRAD libraries which were subjected to sequencing on an Illumina HiSeq2000 sequencer, Illumina, San Diego, CA, United States. Input sequencing data were cleaned and then subjected to Pair-End sequencing on a HiSeq X Ten platform. After mapping the paired-end sequencing reads, the MethylRAD sequencing data were cleaned to obtain high quality reads which were mapped to the reference CCGG/CCWGG sites built with the FspEI sites extracted from the M. notabilis genome using the SOAP program. A site with the sequencing depth of not <3 was determined as a reliable methylation site, and the RPKM was used to determine the relative DNA methylation levels of each site. If the methylation level of a methylation site changed more than two times (P ≤ 0.05) between the infected and healthy samples, it was designated as a significantly differentially methylated site. Then the methylation levels of sites that were localized in the gene regions were summed to evaluate the DNA methylation levels of the genes using the R package edge R (Robinson et al., 2010), and the thresholds of statistically significant difference were found to have a p < 0.05 and log2FC > 1.

BlastN searches against the reference M. notabilis database were performed with online tools (https://morus.swu.edu.cn/morusdb/blast) to provide gene ID for the target genes. Gene ontology (GO) analysis was performed with online tools (https://morus.swu.edu.cn/morusdb/searchgo) based on the gene IDs obtained.

Total RNA was extracted using a TRIzol reagent (Invitrogen Corporation, CA, United States) following the instructions of the manufacturer and then it was digested with DNase I. The cDNA was synthesized from 1 μg of total RNA extracted using oligo (dT)₁₈ primers (Invitrogen Corporation, CA, United States) with reverse transcriptase M-MLV (Promega Corporation, Madison, WI, United States) in 20 ml reactions. The qRT-PCR was performed on the CFX96TM Real-time System (Bio-Rad Laboratories, CA, United States) according to the protocol of the kit (SYBR Premix Ex Taq™) with 20 μl of PCR mix. The PCR mix contains 2 μl of synthesized cDNA, 10 μl of 2 × SYBR Premix Ex Taq, 200 nmol l⁻¹ of gene-specific primers, and 0.4 μl of 50 × ROX Reference Dye II. The PCR mix was preheated to 95°C for 30 s followed by 39 cycles of 95°C for 10 s and 58°C for 30 s. Melting curve analysis was conducted from 65°C up to 95°C followed by a 0.5°C incremental ramp to validate product specificity. The two genes, EF1-a and ACTIN, were used as endogenous controls to quantify the gene expression levels by the 2^−ΔΔCT method (Livak and Schmittgen, 2001). All samples were assayed in at least three biological replicates, and three technical replicates were run for each of them. The primers used for quantitative real-time PCR (qRT-PCR) are given in Supplementary Table 1.

Since the restriction endonuclease, FspEI can recognize 5-methylcytosine (5-mC) and 5-hydroxymethylcytosine (5-hmC) in the CmC and mCDS sites in genomic DNA, the genomic DNA double-strand can be cleaved by FspEI at CCGG and CCWGG methylated sites. The extracted genomic DNA was digested with FspEI, and then the digested DNA was used as a template to perform 28 cycles of PCR amplification. The genomic DNA without being digested was also used as a template to perform 28 cycles of PCR amplification using the same primers and PCR conditions. All the PCR experiments were repeated at least three times. The primers used for semi-quantitative PCR are listed in Supplementary Table 1.

The isolated RNA was used to synthesize cDNA using the reverse transcriptase, M-MLV (Promega). The specific primers (Supplementary Table 2) used for PCR amplifications were designed based on the nucleotide sequence of the gene available from the M. notabilis genome (https://www.ncbi.nlm.nih.gov/genome/?term=Morus) database; the PCR products were separated by electrophoresis and the target DNA fragment was recovered and subcloned into the pMD18-T vector (Invitrogen). After transformation into DH5α, the positive clones were identified and selected for further sequencing. The multiple alignments of the deduced amino acid sequences with the sequences from other plants were conducted using the DNAMAN program. SignalP-5.0 Server (http://www.cbs.dtu.dk/services/SignalP/) was used with the default parameters to predict signal peptides, and TargetP-2.0 Server (http://www.cbs.dtu.dk/services/TargetP/) was used to predict the subcellular localization of proteins. Putative conserved domains were detected using the online NCBI program (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). A phylogenetic tree was generated using the MEGA program by the neighbor-joining method, and bootstrapping was run 1,000 times. The three-dimensional (3D) structure of the protein was generated by the SWISS-MODEL pipeline.

In order to elucidate its subcellular localization, the fusion gene of Mu-GsSRK (GenBank: MN364943.1) and the green fluorescent protein gene (GFP) under the control of 35S was generated and cloned into the binary vector, pBI121 to produce the 35S::Mu-GsSRK-GFP expression vector. The N. benthamiana leaf epidermal cells were infiltrated with Agrobacteria, containing the vector with an optical density of 0.4. After filtration for about 48–72 h, small sections of the infiltrated leaves were excised and mounted in water. Confocal imaging was performed using a Bio-Rad MRC1024 confocal laser scanning microscope (Bio-Rad Microscience Ltd., Manchester, United Kingdom).

The promoter of Mu-GsSRK (designed as pMu-GsSRK) was obtained with specific primers (Supplementary Table 2), and designed based on the nucleotide sequence of the gene available from the M. notabilis genome and used to create the promoter expression vector, pMu-GsSRK::GUS by replacing the 35S promoter and fusing it with the GUS gene. The pMu-GsSRK::GUS vector was then introduced into the Agrobacterium tumefaciens strain, GV3101 which was used to infiltrate N. benthamiana leaves as described previously (Arpat et al., 2012), and β-glucuronidase (GUS) expression in the infiltrated tobacco leaves was assessed by histochemical staining (Jefferson et al., 1987).

The coding region of the Mu-GsSRK gene was amplified and ligated into binary plasmid vector, pBI121 under the control of the 35S promoter. Then, the constructed transgenic plant expression vectors were introduced into the A. tumefaciens strain, GV3101 which was used to transform the wild-type (WT) Arabidopsis plants with the floral dip method (Harrison et al., 2006). Transgenic seeds were selected on selection plates [Murashige and Skoog (MS) media supplemented with 50 mg l⁻¹ of kanamycin], and the T3 generation seeds from independent transgenic lines were used for further functional studies.

Jasmonate (JA) and salicylic acid (SA) treatments were conducted by spraying 100 μmol l⁻¹ of JA or 5 mmol l⁻¹ of SA solution onto the adaxial surface of mulberry leaves. Pseudomonas syringae pv. mori inoculation was performed by spraying the bacterial suspension (10⁸ CFU ml⁻¹) onto the adaxial surface of young leaves of mulberry. As for Colletotrichum dematium inoculation, it was conducted by spraying the conidial suspension (2.5 × 10⁶ ml⁻¹ of conidia) onto the adaxial surface of mulberry leaves. The leaves sprayed with sterilized water were used as controls. All of the inoculated and control mulberry seedlings were incubated in a glass chamber for 48 h to maintain sufficient humidity. P. syringae pv. tomato DC3000 (Pst DC3000) inoculation was performed by injecting 50 μl of Pst DC3000 (10⁵ CFU ml⁻¹) bacterial suspensions into the rosette leaves of 4-week-old Arabidopsis plants. The rosette leaves were detached and inoculated with Botrytis cinerea or Phytophthora capsici using 2-mm-diameter mycelium plugs taken from the actively growing strain colonies. The inoculated leaves were placed in covered Petri dishes to maintain high humidity, and disease symptoms were evaluated by determining the diameter of the lesions 1 week after inoculation. Three biological replicates per treatment were assayed and three technical replicates were performed for each sample.

The leaves inoculated with Pst DC3000 were sampled 72 h after inoculation and ground in sterile water, and the suspension was continuously diluted 10 times with sterile water and spread-plated onto King's B medium. Colonies were counted after 48 h of incubation. All the experiments were conducted at least three times.

The Mu-GsSRK gene was cloned into the vector, pROK2 under the control of the 35S promoter to produce 35S::Mu-GsSRK. The pROK2 vector containing 35S::Mu-GsSRK was introduced into A. rhizogenes strains, K599 that are used for the transformation of mulberry seedlings. Gui sang You 62 mulberry seeds were sown in plastic pots containing vermiculite medium and watered with MS medium and grown in an artificial growth chamber at 25°C with a photoperiod of 16 h of light and 8 h of darkness. About 2 weeks after germination, the seedlings with the first two fully expanded true leaves were selected for agro infiltration. The A. rhizogenes strains, K599 harboring vectors, pROK2 containing 35S::Mu-GsSRK were cultured in yeast extract peptone (YEP) liquid medium containing rifampicin (20 mg l⁻¹) and kanamycin (50 mg l⁻¹), and then were centrifuged when its OD₆₀₀ value reached 0.4 followed by re-suspension in 2-N-morpholino ethanesulfonic acid (MES) buffer (10 mmol l⁻¹ MgCl₂, 10 mmol l⁻¹ MES-KOH and 100 μmol l⁻¹ acetosyringone, pH 5.2). About 0.1 ml of A. rhizogenes harboring vectors pROK2 containing 35S::Mu-GsSRK suspension was injected into the junction of true leaves and cotyledons as described previously (Meng et al., 2019). At the same time, the plant injected with A. rhizogenes strains K599 harboring empty vectors, pROK2 were used as the control group. About 1 month later, when the hairy roots developed well, they were detected by PCR and the original roots were cut off.

To obtain mulberry yellow dwarf phytoplasma (16SrI-B)-infected leafhoppers (Hishmonus sellatus), the leafhoppers were transferred to phytoplasma–infected mulberry plants for 2 weeks. A. thaliana and mulberry seedlings were inoculated with mulberry phytoplasma by the sap-feeding method with phytoplasma-infected leafhoppers (Sugio et al., 2011), and the seedlings challenged with uninfected leafhoppers were used as controls. Each of the A. thaliana or mulberry seedling was challenged with six adult leafhoppers for 6 days. All the seedlings were incubated in a growth chamber at 26°C, 60% humidity, and under 12 h of light. Four weeks after inoculation, the DNA was extracted following the method described previously (Prince, 1993), and phytoplasma concentration in the plants was determined using real-time PCR described previously (Christensen et al., 2004). The primers and TaqMan probes used to amplify the 16S rRNA gene of phytoplasma, Arabidopsis, and mulberry are listed in Supplementary Table 2. Three biological replicates per treatment were assayed and three technical replicates were performed for each sample.

Using RNA-seq, a total of 48.9 and 48.7 million clean reads were obtained in the healthy and phytoplasma-infected mulberry (M. multicaulis Perr.) leaf mRNA libraries, respectively. The expression level of genes identified was calculated, and a total of 1,168 genes were found to be differentially expressed between phytoplasma-infected and healthy leaves, among which 769 genes were upregulated and 399 genes were downregulated (Figure 1A). All the detected DEGs are shown in Supplementary Table 1 and the top 40 up and downregulated ones are listed in Table 1. The GO analysis of the DEGs classified them into 11 functional categories (Figure 1B). The first category of gene was involved in stress and environment response, and the genes associated with catabolic process belonged to the second category, and the third category included the genes involved in transcription and post-transcription regulation. The other DEGs belonged to categories, such as secondary metabolite biosynthesis, signaling pathway, growth and development, cell component, and cycle process. The detailed functional GO terms of the DEGs are given in Supplementary Table 3. Interestingly, 6% of the DEGs were found to be associated with signal transduction pathways, indicating that diverse signal transduction pathways were involved in the response of mulberry to phytoplasma infection. Therefore, the regulatory networks involved in response to phytoplasma infection in mulberry plants are found to be complex.

To validate the expression profiles obtained by RNA-seq, a qRT-PCR analysis was performed for 10 genes including the upregulated and downregulated ones in response to phytoplasma infection (Figure 2). The data showed that all the selected genes exhibited similar changes in expression levels between qRT-PCR and RNA-seq results, indicating that the differential expression profiles of these genes obtained by RNA-seq are reliable.

There were 23,056 CCGG and 8,520 CCWGG DNA methylation sites detected in the phytoplasma-infected leaves, and a total of 22,733 CCGG and 7,856 CCWGG DNA methylation sites were found in healthy leaves. When the MethyIRAD reads were aligned onto a unique locus in different regions of the genome, the genome-wide methylation pattern was obtained, and the results showed that the distributions of methylation sites in different components of the genome were differential. However, it was found that the distribution patterns of methylation sites at different elements of genomes were similar in the phytoplasma-infected and healthy mulberry leaves (Figures 3A–F). Most CCGG DNA methylation sites were mainly enriched in the exon regions followed by the intergenic, intron, and upstream regions (Figure 3G). As for the CCWGG DNA methylation sites, the major sites were mainly enriched in the intergenic regions followed by the exon, upstream, and intron regions (Figure 3H).

In this study, a total of 3,676 and 1,366 differentially methylated sites in CCGG and CCWGG, respectively, were found between phytoplasma-infected and healthy mulberry leaves. Most of the differentially methylated sites were found in the intergenic and exon regions, and only about 10% of the differentially methylated sites were found in the upstream and intron regions, respectively (Figure 4). The genes with differentially methylated sites in phytoplasma-infected and healthy mulberry leaves were screened and termed as DMGs. In total, 935 and 347 DMGs were identified in CCGG (Supplementary Table 4) and CCWGG sites (Supplementary Table 5), respectively. The results revealed that 511 DMGs in CCGG sites as well as 220 DMGs in CCWGG sites showed upregulation and 424 DMGs in CCGG sites as well as 127 DMGs in CCWGG sites showed downregulation in the phytoplasma-infected leaves compared to healthy leaves. Interestingly, there were 29 genes which were differentially methylated at both CCGG and CCWGG sites. These results suggested that phytoplasma infection induced the changes of methylation level of some genes.

To validate the results obtained with MethylRAD-seq data, six genes were selected for analysis by PCR. The genomic DNA digested with FspEI was used as the template for semi-quantitative PCR, and there were more amplification products obtained from the samples, in which the genomic DNAs were lower methylated than those from the samples in which the genomic DNAs were highly methylated at CCGG or CCWGG sites (Figure 5). The PCR results showed a high degree of consistency with the MethylRAD-seq data, indicating that the genomic DNA methylation results obtained by MethylRAD are reliable.

Based on the RNA-Seq and MethylRAD sequencing data, the methylation and expression levels of the genes in the infected and healthy leaves were compared. Though a lot of genes were methylated or expressed differently, only 51 genes were differentially methylated and expressed between the infected and healthy leaves. Among these differentially expressed and methylated genes (DEMGs), there were 40 and 11 genes with CCGG and CCWGG differential methylation sites, respectively (Tables 2, 3). The correlation between DNA methylation and gene expression showed that the expression levels of 33 genes were negatively correlated with their methylation levels, while the expression levels of 18 genes were positively correlated with their methylation levels. Moreover, it was found that most of the DEMGs were differentially methylated in the exon regions, and only 5 DEMGs were differentially methylated in the intron regions and one was differentially methylated in the intergenic region (Supplementary Tables 6, 7).

The GO analysis of these DEMGs indicated that these genes were classified into nine functional categories (Figure 6). The first category of genes was involved in metabolism followed by the second category of genes which was associated with stress response, and the third category of genes was involved in signal transduction. The other genes were classified into transportion, RNA processing, regulation of transcription, development, and unknown categories.

Integrated analysis showed that the expression of LOC21383997, which was annotated as Mu-GsSRK was increased, but its methylation level was decreased in the phytoplasma-infected leaves, and these results were confirmed by qRT-PCR (Figure 2) and semi-quantitative PCR analyses (Figure 5). The protein encoded by the Mu-GsSRK gene contains 781 amino acids, and its predicted molecular weight (MW) and isoelectric point (pI) were found to be 87.6 kDa and 6.28, respectively. Multiple sequence alignments revealed the protein-shared homology regions with GsSRKs in various species (Figure 7A). Similar to other typical G-type lectin receptor-like kinases (LecRLKs), Mu-GsSRK protein has some putative conserved domains, such as protein kinase active sites, ATP binding sites, substrate binding sites, and two activation loops. Besides these conserved kinase domains, the protein also contains a B lectin (bulb-type mannose-specific lectin) region and an S-locus-glycoprotein region at the N terminus. In addition, a plant PAN/APPLE-like domain associated with various biological functions by mediating protein–carbohydrate or protein–protein interactions was detected in the Mu-GsSRK protein (Figure 7B).

Phylogenetic analysis of Mu-GsSRK and GsSRKs from other plants showed that there was closest homology between Mu-GsSRK and GsSRK in M. notabilis (Figure 8A). The structural properties of the Mu-GsSRK protein was predicted using SWISS–MODEL, and the result showed that the protein contained 54.93% random coil, 23.30% alpha helix, and 21.77% extended strands and had a β-barrel structure composed of a number of strands in their N-terminal lectin domains (Figure 8B). The prediction of N-terminal extension suggested that Mu-GsSRK contained an obvious signal peptide, but no obvious sublocalization sequence was detected. In order to elucidate its subcellular localization, Mu-GsSRK was fused to a green fluorescent protein gene and introduced into N. benthamiana leaves. The green fluorescent signal was detected within the plasma membrane suggesting the localization of Mu-GsSRK in the plasma membrane (Figure 8C). Therefore, as a protein kinase, Mu-GsSRK may play a role in sensing and transmitting signals arising from different biotic stress conditions.

To explore the function of the Mu-GsSRK gene, its expression pattern in different tissues and organs was analyzed by qRT-PCR. The results showed that Mu-GsSRK was expressed ubiquitously in the investigated parts of mulberry plants, but its expression level was higher in the leaves than that in other parts (Figure 9A). Moreover, the induced expression pattern of Mu-GsSRK was explored by challenging the mulberry seedlings with P. syringae pv. mori and C. dematium and by treating the seedlings with JA and SA, respectively. The results indicated that the expression level of Mu-GsSRK was increased in the leaves challenged by P. syringae pv. mori or C. dematium. Furthermore, it showed that the exogenous application of JA or SA enhanced the expression of Mu-GsSRK in the leaves (Figure 9B). In addition, pMu-GsSRK was fused with the GUS gene and then was transiently expressed in tobacco leaves. The results of GUS staining showed that GUS activity was induced in the leaves infected by B. cinerea or Pst DC3000, as well as in the leaves upon inoculation with JA or SA (Figure 9C). These data indicated that Mu-GsSRK may be associated with stress resistance, and JA and SA may modulate its expression in mulberry.

To explore whether DNA methylation was involved in the expression change of the Mu-GsSRK gene, genomic DNA digested with FspEI was used as the template for PCR analysis. The semi-quantitative PCR results indicated that the Mu-GsSRK gene was methylated at CCGG or CCWGG sites, but the methylation level of Mu-GsSRK did not differ significantly among the different parts of mulberry plants (Figure 10A). Therefore, the expression difference of Mu-GsSRK in the different parts is not caused by DNA methylation. However, the PCR results indicated that the methylation level of the Mu-GsSRK gene was reduced significantly in the leaves inoculated with P. syringae pv. mori or C. dematium or treated with SA, but its DNA methylation level may not be affected by JA treatment. So, the expression change of Mu-GsSRK in the leaves inoculated with Pst DC3000 or C. dematium or treated with SA may be also caused by its methylation level change (Figure 10B).

To explore the defensive role of Mu-GsSRK, wild-type (WT) and transgenic Mu-GsSRK plants were challenged by Pst DC3000. Three days post-inoculation (DPI), the WT plants showed severe disease symptoms. In contrast, there were no evident disease symptoms observed in the leaves of transgenic plants. Meanwhile, the transgenic Mu-GsSRK A. thaliana plants were challenged by B. cinerea or Phytophthora capsici to explore the possible role of Mu-GsSRK in defense against fungal pathogens. Four DPI, B. cinerea or P. capsici were successfully colonized on the inoculation leaf surface of WT plants, and obvious necrotic lesions were observed on the leaves. On the contrary, mild disease symptoms were observed on the inoculated leaf surface of Mu-GsSRK-overexpressing plants (Figure 11A). In addition, the bacterial populations of Pst DC3000 strains in the inoculated leaves were determined, and the result revealed that the population of the strains in the leaves of WT plants was significantly higher than that in the ones of Mu-GsSRK-overexpressing plants at 3 DPI (Figure 11B). Moreover, the lesion areas in the leaves of WT plants inoculated with B. cinerea and P. capsici were larger than those in the leaves of Mu-GsSRK-overexpressing plants (Figure 11C). These results indicate that the overexpression of Mu-GsSRK in Arabidopsis confer enhanced resistance to Pst DC 3000, B. cinerea, and P. capsici, and the Mu-GsSRK gene may play an important role in disease resistance.

To examine whether the expressions of some defense-related genes were elevated in the Mu-GsSRK-overexpressing A. thaliana plants, the expression levels of the pathogenesis-related protein 1 (PR-1) gene, plant defensin gene (PDF1.2), and cytochrome P450 protein CYP82C2 gene in the transgenic Mu-GsSRK plants were evaluated. The data showed that the expression levels of PR-1, PDF1.2, and CYP82C2 genes were all very low both in the WT and transgenic Mu-GsSRK plants in the absence of Pst DC3000 inoculation, indicating that the constitutive overexpression of the Mu-GsSRK gene might not affect the basal expression levels of these defense-related genes. However, after Pst DC3000 treatment, the expression levels of PR-1, PDF1.2, and CYP82C2 genes were higher in the Mu-GsSRK-overexpressing plants than in WT plants at 12 and 24 h post-inoculation (Figure 11D). These results demonstrated that Mu-GsSRK may have a positive role in the regulation of defense gene expressions, but pathogen induction is a necessary step for defense gene expression.

Since the efficient genetic transformation and regeneration system has not yet been established in mulberry trees, in order to examine the role of Mu-GsSRK in the defense response to phytoplasma, transgenic Mu-GsSRK mulberry plants with hairy roots were generated (Supplementary Figure 1A). qRT-PCR showed that the Mu-GsSRK gene was constitutively highly expressed in the transgenic hairy roots (Supplementary Figure 1B). After the confirmation of transgenic Mu-GsSRK hairy roots, the original roots were cut off and the seedlings carrying the transgenic hairy roots were used in the follow-up of phytoplasma challenge experiments using the sap-feeding method. Four weeks post-challenge, the mulberry seedlings carrying the non-transgenic hairy roots showed severe symptoms of Witches' broom disease, such as short internodes, leaf curling, and axillary bud sprouting. Contrarily, the mulberry seedlings carrying the transgenic Mu-GsSRK hairy roots showed mild dwarfism symptoms and enhanced resistance to phytoplasma (Figures 12A,B). Quantitative PCR analysis showed that the concentrations of phytoplasmas in the aboveground parts and roots of the seedlings carrying the non-transgenic hairy roots were about 22.5 × 10⁵ and 1.3 × 10⁵ cells μg⁻¹ plant DNA, respectively. While the concentrations of phytoplasmas in the aboveground parts and roots of the seedlings carrying the transgenic Mu-GsSRK hairy roots were about 6.0 × 10⁵ and 0.4 × 10⁵ cells μg⁻¹ plant DNA, respectively, there was no phytoplasma detected in the plants challenged by healthy leafhoppers (Figure 12C). Therefore, the concentrations of phytoplasmas in the aboveground parts and roots of the seedlings carrying the non-transgenic hairy roots were higher than those in the seedlings carrying the transgenic Mu-GsSRK hairy roots, and the overexpression of the Mu-GsSRK gene in the hairy roots partially inhibited the growth of phytoplasmas and led to enhanced resistance to phytoplasma.

As an important and conserved epigenetic modification, DNA methylation is associated with many important biological processes, and has been extensively studied in recent years. However, few studies have focused on mulberry, and this is the first report about the genome-wide DNA methylation profiles in mulberry. Previous reports showed that DNA methylation is not randomly distributed in genomes. In plants, the gene body regions are often highly methylated, while the transcriptional start sites (TSS) and transcriptional termination sites (TTS) mostly lack DNA methylation. Moreover, many intergenic regions also show hypermethylation; on the contrary, most promoter regions are hypomethylated (Zhang et al., 2006). The results presented here showed that DNA methylation is enriched in intergenic and gene body regions, whereas the promoter regions are also hypomethylated in mulberry genomes. In addition, it was also observed that the methylation level in exons was higher than that in introns and CCGG methylation sites (Figure 4). These results indicated that DNA methylation profiles of mulberry are analogous to those of other plants, and the methylation patterns among different plant species may be conservative. It is widely accepted that DNA methylation is a major transcription silencing pathway, which is negatively correlated with gene expression in plants (Zhang et al., 2006). However, it has also been reported that DNA methylation is positively correlated with gene expression (Lou et al., 2014). Analysis of gene expression differences between the healthy and infected leaves indicated that DNA hypermethylation is associated with activation of several genes during phytoplasma infection (Tables 2, 3). Therefore, the mechanism of gene regulation mediated by DNA methylation is very complicated and further study is necessary to elucidate it in detail.

Previous studies have shown that DNA methylation patterns can be changed when plants are infected with pathogens, and it was reported that the DNA methylation levels decreased in Paulownia plantlets infected with phytoplasma (Cao et al., 2014). However, no significant difference in the average DNA methylation level between healthy and phytoplasma-infected mulberry samples was detected (Figure 3). The transcriptome data also indicated that phytoplasma infection did not lead to significant changes in the expression of methyltransferase and demethylase in the infected samples. Although the average DNA methylation level was not changed, the methylation levels of more than 1,253 genes were changed significantly in the infected samples (Supplementary Tables 4, 5). So the methylation levels of certain genes may be dynamically regulated to control plant resistance against phytoplasma infection, but the genome was closely monitored to maintain stability. However, whether there is a link between loci-specific methylation and phytoplasma infection remains to be established.

It has been reported that DNA methylation is associated with gene expression changes in response to phytoplasma infection (Jagoueix-Eveillard et al., 2001; Pracros et al., 2006; Ahmad et al., 2012; Cao et al., 2014; Pagliarani et al., 2020). In this study, 51 DEMGs were identified in the infected leaves. Among these, there were some receptor kinase genes, such as LRR receptor-like serine/threonine-protein kinase, Mu-GsSRK, and serine/threonine-protein kinase SAPK2e (Tables 2, 3), which play a central role in signaling during the pathogen recognition in plants (Bouwmeester and Govers, 2009; Wang and Bouwmeester, 2017). By adjusting the methylation levels of these genes, plants can regulate their expression levels to enhance the perception of phytoplasma infection. On the other hand, phytoplasma may change the methylation level of these genes to affect their expression, thus avoiding plant perception (Zhai et al., 2010; Cao et al., 2014). In addition, there were 22 DEMGs involved in metabolic, growth and developmental processes were found in the infected leaves; the changes of these genes may disturb the normal metabolic growth and developmental processes leading to various symptoms. Although these symptoms may hinder the normal growth and development of mulberry plants, they may promote parasite multiplication of phytoplasma in the infected plants. Moreover, among these DEMGs, there were 13 genes associated with stress response (Tables 2, 3). This suggests that mulberry plants can change the expression of some disease resistance genes by altering their methylation status, so as to improve their resistance to phytoplasma. Therefore, our evidence supported that the DNA dynamic methylation levels were associated with expression changes of the genes involved in response to phytoplasma infection. Although there were many genes which were simultaneously methylated and expressed differently in the infected and healthy leaves, there is a wide variety of mechanisms that control gene expression, and DNA methylation may be the reason for the change of gene expression, but the change may be affected by other factors rather than DNA methylation.

The lectin receptor-like kinases (LecRLKs) are categorized into three sub-classes: G-, L-, and C-type depending on the features of their N-terminal lectin domains (Vaid et al., 2012), and the Mu-GsSRK protein has some domains conserved in G-type LecRLKs (Figure 7B). So it might have similar biological functions with other G-type LecRLKs. Due to the resemblance of the extracellular domain with the lectin protein, known to bind to fungal and bacterial cell wall components, LecRLKs are hypothesized to predominantly participate in biotic stress tolerance (Wang and Bouwmeester, 2017). There were many G-type LecRLKs that have been found to be associated with the responses to pathogen infections, and some LecRLKs have been reported to confer resistance to a variety of pathogens (Chen et al., 2006; Sanabria et al., 2012; Vaid et al., 2013; Ranf et al., 2015). However, to the best of our knowledge, the current study is the first to report that G-type LecRLK genes are associated with the response to phytoplasma infection.

It was reported that some members of the LecRLK family are located on the plasma membrane and can sense signals arising from different biotic stress conditions and transmit the stress signal to the nucleus (Vaid et al., 2013). So some LecRLK family genes were proposed to play a role in mediating and strengthening the cell wall–plasma membrane (CW–PM) links and continuum which is essential for defense against pathogens (André et al., 2005; Bouwmeester et al., 2011; Singh et al., 2012; Ranf et al., 2015). Our results showed that the Mu-GsSRK protein is localized on the plasma membrane (Figure 8C), and it can significantly increase the expressions of disease resistance-related genes upon pathogenic infection (Figure 11C) and enhanced resistance to the pathogens (Figure 11A). Moreover, our results showed that the overexpression of Mu-GsSRK in the roots not only enhanced the roots but also the resistance of the aboveground parts to phytoplasma (Figure 12). Therefore, the Mu-GsSRK protein may be able to sense the signals arising from phytoplasma infection and transmit the signals to the nucleus and regulate the expressions of disease resistance-related genes, and the signal might be transmitted a long distance and stimulate the immune response of the whole plant.

It was reported that there was an increase in endogenous SA in the phytoplasma-infected plants (Dermastia, 2019). In addition, our results showed that the methylation level of the Mu-GsSRK gene was reduced significantly in the leaves treated with SA (Figure 10B). Therefore, the endogenous SA may be increased in the phytoplasma-infected mulberry leaves, which may reduce the methylation level and enhance the expression of Mu-GsSRK. Moreover, our results showed that the activity of the Mu-GsSRK promoter was induced by SA (Figure 9C). So, the activities of the Mu-GsSRK promoter may be induced by the increase of endogenous SA caused by phytoplasma infection and enhance the expression of the Mu-GsSRK gene. In response to phytoplasma infection, SA may be an important factor in the regulation of methylation and the expression level of Mu-GsSRK in mulberry; further study is required to elucidate the precise mechanism behind this process.