ARTICLE   Open Access    

Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock

More Information
  • Apple (Malus domestica) is an important multi-functional horticultural crop. However, thus far, there have been few reports concerning the molecular mechanism by which interstock regulates rootstock and scion growth. MicroRNAs (miRNAs) are a type of post-transcriptional regulator that regulates apple growth and development, fruit quality, hormone signal transduction, and stress response. Interstock grafting resulted in differences in miRNA distribution between apple tissues. The regulatory roles of miRNAs in the grafting of apples are not yet clear. In this study, 15 libraries were constructed using the phloem of an apple grafted on two types of interstock and rootstock by high-throughput sequencing. A total of 281 miRNAs were catalogued into 80 families, and 159 novel ones were identified. Compared with the control, the grafting combination with the interstock resulted in 79 differentially expressed miRNAs (DEMs) in the scion, with 36 up-regulated and 43 down-regulated miRNAs, and 57 DEMs in the rootstock, including 22 up-regulated and 35 down-regulated miRNAs. Grafting with the dwarfing interstock led to DEMs in the whole plant, including a decrease in the expression of mdm-miR156 in the scion and mdm-miR172 in the rootstock. Predictive analysis of the DEMs and their target genes suggested that miRNAs mediate scion growth through several aspects, such as plant hormone synthesis and signal transduction, and nutrient absorption and balance. Through combined miRNA and mRNA analysis, the dwarfing effect of the interstock may affect the expression of genes related to the starch and sucrose metabolism pathways in the rootstock and the plant hormone signal transduction pathway in the scion. The available evidence facilitates a better understanding of the role of miRNAs in the response of apples after grafting interstock. One possible reason for the stunting of apple trees and the promotion of early flowering caused by the dwarfing interstock is the decrease in expression levels of mdm-miR156 and mdm-miR172.
  • Rice (Oryza sativa L.) is a world staple crop that feeds over half of the world's population[1,2]. In recent years, high-temperature events are becoming more frequent and intensive as a result of global warming, which can severely affect rice grain yield and quality[3,4]. During flowering and grain filling stages, high-temperature stress can result in a significant reduction in seed setting rate and influence amylose content, starch fine structure, functional properties and chalkiness degree of rice[57]. Transcriptome and proteome analysis in rice endosperm have also been used to demonstrate the differences in high-temperature environments at gene and protein expression levels[810]. In addition, as an important post-translational modification, protein phosphorylation has proven to be involved in the regulation of starch metabolism in response to high-temperature stress[11]. However, little is known about whether protein ubiquitination regulates seed development under high-temperature stress.

    Ubiquitination is another form of post-translational modification that plays key roles in diverse cellular processes[12]. Several reports have described the functions of ubiquitination in rice defense responses based on ubiquitome analysis. Liu et al.[13] investigated relationships between ubiquitination and salt-stress responses in rice seedlings using a gel-based shotgun proteomic analysis and revealed the potential important role of protein ubiquitination in salt tolerance in rice. Xie et al.[14] identified 861 peptides with ubiquitinated lysines in 464 proteins in rice leaf cells by combining highly sensitive immune affinity purification and high resolution liquid chromatography-tandem mass spectrometry (LC-MS/MS). These ubiquitinated proteins regulated a wide range of processes, including response to stress stimuli. A later study revealed the relationships between ubiquitination status and activation of rice defense responses, and generated an in-depth quantitative proteomic catalog of ubiquitination in rice seedlings after chitin and flg22 treatments, providing useful information for understanding how the ubiquitination system regulates the defense responses upon pathogen attack[15]. Although many studies have shown that ubiquitination plays improtant roles in the heat response of plant[16,17], there has been little systematic discussion on the ubiquitome of rice endosperm in the context of global climate change.

    In this study, we examine the high-temperature induced ubiquitination change in two rice varieties with different starch qualities, through a label-free quantitative ubiquitome analysis. This study provides a comprehensive view of the function of ubiquitination in high-temperature response of rice developing seed, which will shed new light on the improvement of rice grain quality under heat stress.

    Two indica rice varieties with different starch quality, 9311 and Guangluai4 (GLA4), were used as materials. 9311 is a heat-sensitive variety, which displays low amylose content with good starch quality; while GLA4 is known to be the parental variety of HT54, an indica breeding line with heat tolerance, and thus GLA4 is possibly heat tolerant, which shows high amylose content with poor starch quality[18,19]. Rice growth conditions, sample treatment and collection were conducted as previously described[11].

    Husk, pericarp and embryo were detached from immature rice grains on ice[20]. Rice endosperm was then ground with liquid nitrogen, and the cell powder was sonicated 10 rounds of 10 s sonication and 15 s off-sonication on ice in lysis buffer (6 M Guanidine hydrochloride, pH 7.8−8.0, 0.1% Protease Inhibitor Cocktail) using a high intensity ultrasonic processor. Following lysis, the suspension was centrifuged at 14,000 g for 40 min at 4 °C to remove the debris. The supernatant was collected, and the protein concentration was estimated using BCA assay (Pierce BCA Protein assay kit, Thermo Fisher Scientific, Waltham, MA, USA) before further analysis.

    The protein mixture was reduced by DTT with the final concentration of 10 mM at 37 °C for 1 h, alkylated by iodoacetamide with a final concentration of 50 mM at room temperature in the dark for 0.5 h, and digested by trypsin (1:50) at 37 °C for 16 h. Then the sample was diluted by adding trifluoroacetic acid (TFA) to the final concentration of 0.1%. The enzymatic peptides were desalted on a Sep-Pak C18 cartridge (Waters, Milford, MA, USA), concentrated by lyophilization and reconstituted in precooled IAP buffer (50 mM MOPS-NaOH PH 7.2, 10 mM Na2HPO4, 50 mM NaCl) for further analysis.

    The peptides solution was incubated with prewashed K-ε-GG antibody beads (PTMScan Ubiquitin Remnant Motif (K-ε-GG) Kit), and gently shaken at 4 °C for 1.5 h. The suspension was centrifuged at 2,000 g for 30 s, and the supernatant was removed. The Anti-K-ε-GG antibody beads were washed with IAP Buffer three times and with ddH2O three times. The peptides were eluted from the beads with 0.15% trifluoroacetic acid (TFA). Finally, the eluted fractions were combined and desalted with C18 Stage Tips.

    LC-MS/MS analysis were performed using the methods of Pang et al.[11]. Raw mass spectrometric data were analyzed with MaxQuant software (version 1.3.0.5) and were compared with the indica rice protein sequence database (Oryza sativa subsp. indica-ASM465v1). Parameters were set according to Pang et al.[11]. All measurements were obtained from three separate biological replicates.

    Quantification of the modified peptides was performed using the label-free quantification (LFQ) algorithm[11]. Differentially ubiquitinated sites (proteins) in response to high-temperature were identified by Student's t-test (p < 0.05, log2(fold-change) > 1) with at least two valid values in any condition or the ubiquitination sites that exhibited valid values in one condition (at least two of three replicates) and none in the other.

    Subcellular localization was performed using CELLO database (http://cello.life.nctu.edu.tw). Gene Ontology (GO) annotation proteome was derived from the AgriGO (http://bioinfo.cau.edu.cn/agriGO/). The differential metabolic profiles were visualized with MapMan software (version 3.6.0RC1). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation was performed by using KEGG Automatic Annotation Server (KAAS) software. A p-value of < 0.05 was used as the threshold of significant enrichment. SWISS-MODEL was used to generate the tertiary structure of GBSSI (SWISS-MODEL, http://swissmodel.expasy.org/). The figures were annotated with Adobe Illustrator (Adobe Systems, San Jose, CA, USA).

    To elucidate how high-temperature stress influences rice developing endosperm at the ubiquitination level, a label-free analysis was performed to quantify ubiquitome from two indica rice varieties under normal (9311-C and GLA4-C) and high-temperature conditions (9311-H and GLA4-H). The distribution of mass error of all the identified peptides was near zero and most of them (71.5%) were between -1 and 1 ppm, suggesting that the mass accuracy of the MS data fits the requirement (Fig. 1a). Meanwhile, the length of most peptides distributed between 8 and 42, ensuring that sample preparation reached standard conditions (Fig. 1b).

    Figure 1.  Characteristics of the ubiquitinated proteome of rice endosperm and QC validation of MS data. (a) Mass error distribution of all identified ubiquitinated peptides. (b) Peptide length distribution. (c) Frequency distribution of ubiquitinated proteins according to the number of ubiquitination sites identified.

    In all endosperm samples, a total of 437 ubiquitinated peptides were identified from 246 ubiquitinated proteins, covering 488 quantifiable ubiquitinated sites (Supplemental Table S1). Among the ubiquitinated proteins, 60.6% had only one ubiquitinated lysine site, and 18.7%, 8.1%, 5.3%, or 7.3% had two, three, four, or five and more ubiquitinated sites, respectively. In addition, four proteins (1.6%, BGIOSGA004052, BGIOSGA006533, BGIOSGA006780, BGIOSGA022241) were ubiquitinated at 10 or more lysine sites (Fig. 1c, Supplemental Table S1). The proteins BGIOSGA008317 had the most ubiquitination sites with the number of 16. It was noted that besides ubiquitin, two related ubiquitin-like proteins NEDD8 and ISG15 also contain C-terminal di-Gly motifs generated by trypsin cleavage, and the modifications of these three proteins cannot be distinguished by MS[21]. Here, the di-Gly-modified proteome therefore represents a composite of proteins modified by these three proteins. However, the sites from NEDD8 or ISG15 modifications were limited because they mediate only a few reactions in cells[21].

    To better understand the lysine ubiquitome changes in rice endosperm induced by high-temperature, we performed a Gene Ontology (GO) functional annotation analysis on all identified ubiquitinated proteins (Fig. 2a). In the biological process GO category, 'metabolic process' and 'cellular process' were mainly enriched, accounting for 75.1% and 74.1% of ubiquitinated proteins, respectively. In addition, 34.6% proteins were associated with 'response to stimulu', emphasizing the regulatory role of ubiquitination modification in response to high-temperature stress. From the cellular component perspective, ubiquitinated proteins were mainly associated with 'cellular anatomical entity' (99.4%), 'intracellular' (84.4%) and 'protein-containing complex' (29.9%). The molecular function category result suggested that these proteins were largely involved in 'binding' (62.7%), 'catalytic activity' (43.4%) and 'structural molecule activity' (16.5%). Furthermore, subcellular location annotation information indicated that 34.7%−39.4% proteins were located in the cytoplasm, and other were mostly located in the nucleus (23.5%−27.7%), plasma membrane (9.4%−11.4%), and chloroplast (9.6%−12.8%) (Fig. 2b). It is noteworthy that the ubiquitinated proteins located in the cytoplasm were decreased in high-temperature environments in both varieties.

    Figure 2.  Analysis of ubiquitinated proteins and motifs. (a) Gene ontology (GO) functional characterization of ubiquitinated proteins. (b) Subcellular localization of ubiquitinated proteins. From the inside out, the ring represents 9311-C, 9311-H, GLA4-C and GLA4-H, respectively. (c) Motif enrichment analysis of ubiquitinated proteins.

    The following two significantly enriched motifs from all of the identified ubiquitinated sites were identified using MoDL analysis: [A/S]xKub and Kubxx[E/Q/R/V]x[E/G/L/P/Q/R/Y], which covered 84 and 100 sequences, respectively (Fig. 2c). Further analysis showed that the conserved alanine (A) and glutamic acid (E) were included in upstream and downstream of the ubiquitinated lysine sites in rice endosperm. A similar phenomenon also occurred in rice leaf[14], wheat leaf[22], and petunia[23], indicating that alanine (A) and glutamic acid (E) were likely to be the specific amino acids in conserved ubiquitination motifs in plants. Additionally, serine (S) was enriched at the position -2 (upstream) of the ubiquitinated lysine, while various amino acids such as arginine (R), glutamic acid (E), glutamine (Q), valine (V) were found at positions +3 and +5 (downstream).

    To detect possible changes in rice endosperm ubiquitome attributable to high-temperature stress, we then performed LFQ analysis on all quantifiable ubiquitination sites within our dataset. As shown in Fig. 3a, more ubiquitinated proteins, peptides and sites were detected in the treatment groups (9311-H and GLA4-H), suggesting that exposure to high-temperature stress may increase the ubiquitination events in rice endosperm. Only 282 common ubiquitinated sites in 158 proteins were quantifiable for all sample groups due to reversible ubiquitination induced by high-temperature (Fig. 3b). Principal component analysis (PCA) showed that three repeats of each sample clustered together, and four groups were clearly separated (Fig. 3c). Furthermore, the differentially expression profiles of ubiquitination sites (proteins) in 9311 and GLA4 under high-temperature stress were depicted to further understand the possible changes (Fig. 3d). Where LFQ values were missing, the data were filtered to identify those ubiquitination sites with a consistent presence/absence expression pattern. These analyses yielded 89 ubiquitination sites that were only present in 9311-H and six that were only present in 9311-C (Fig. 3d, Supplemental Table S2). Similarly, 51 differentially expressed ubiquitination sites were present in GLA4-H and 13 ubiquitination sites only occurred in GLA4-C (Fig. 3d & Supplemental Table S3). Beyond that, a total of 113 and 50 significantly changed ubiquitination sites (p < 0.05, log2(fold-change) >1) were screened out in 9311 and GLA4, respectively (Fig. 3d, Supplemental Tables S4 & S5). For subsequent comparative analysis, the ubiquitination expression profiles with consistent presence/absence and ubiquitination sites with significant differences in statistical testing were combined and named as 9311-Up, 9311-Down, GLA4-Up, and GLA4-Down, respectively (Fig. 3d). The number of significantly up-regulated ubiquitination sites was far greater than down-regulated ubiquitination sites in both 9311 and GLA4 varieties. These findings indicated that high temperature not only induced the occurrence of ubiquitination sites, but also significantly upregulated the intensity of ubiquitination. Beyond that, the magnitude of the up-regulation in 9311 was higher than that in GLA4 (Fig. 3b & d), indicating that the ubiquitination modification of heat-sensitive varieties was more active than heat-resistant varieties in response to high-temperature stress.

    Figure 3.  A temperature regulated rice endosperm ubiquitome. (a) The number of ubiquitinated proteins, peptides and sites detected in four group samples. (b) Venn diagram of ubiquitination sites (proteins) detected in four group samples. (c) PCA based on ubiquitination intensity across all four sample groups with three biological repetitions. (d) Differentially expression profiles of ubiquitination sites (proteins) in 9311 and GLA4 under high-temperature stress. The expression profiles of selected ubiquitination sites (p < 0.05, log2(fold-change) >1) were normalized using the Z-score and presented in a heatmap.

    To further investigate the ubiquitination regulatory pattern under high temperature stress in two varieties, four groups with significantly regulated sites were analyzed. There were 37 ubiquitination sites showed the same regulatory trend in 9311 and GLA4, accounting for 17.8% and 32.5% of the total differentially expressed sites in 9311 and GLA4, respectively. Among them, 36 ubiquitination sites were upregulated and one site was downregulated (Fig. 4a). In addition, 159 upregulated ubiquitination sites and three downregulated sites were only present in 9311, while 53 upregulated sites and 15 downregulated sites were only present in GLA4. Moreover, nine ubiquitination sites showed opposite regulatory trends in 9311 and GLA4. A similar regulatory trend of ubiquitination proteins is shown in Fig. 4b. It is noted that some proteins had both upregulated and downregulated ubiquitination sites (Supplemental Tables S6 & S7), indicating that significant differences in ubiquitination were, to some extent, independent of protein abundance.

    Figure 4.  Comparison of differentially ubiquitinated sites and proteins in 9311 and GLA4 under high-temperature stress.

    To understand the function of ubiquitination in response to the high-temperature stress of rice endosperm, we conducted GO enrichment-based clustering analysis of the differentially ubiquitinated proteins in 9311 and GLA4 at high temperature, respectively (Fig. 5). In the biological process category of 9311, proteins were relatively enriched in the carbohydrate metabolic process, polysaccharide metabolic process, starch biosynthetic process, cellular macromolecule localization, protein localization, intracellular transport, and phosphorylation (Fig. 5). For the molecular function analysis, we found that the proteins related to kinase activity, nucleotidyltransferase activity, phosphotransferase activity, and nutrient reservoir activity were enriched (Fig. 5). The two principal cellular components were intrinsic component of membrane and integral component of membrane (Fig. 5). There was no significantly enriched GO term in the GLA4 group due to the dataset containing relatively few proteins, and thus, further enrichment analysis was conducted on the proteins that were common to both varieties. The results showed that proteins were over-represented in carbon metabolism, including starch biosynthesis and metabolism, glucan biosynthesis and metabolism, and polysaccharide biosynthesis and metabolism (Fig. 5), indicating the importance of carbohydrate synthesis and metabolism in the ubiquitination regulatory network.

    Figure 5.  Enrichment analysis of differentially expressed ubiquitinated proteins based on Gene Ontology (GO) terms.

    To identify pathways which were differentially ubiquitinated under high-temperature stress, the KEGG pathway-based clustering analysis was conducted. The results showed that the differentially ubiquitinated proteins in both 9311 and GLA4 were mostly abundant in the pathways of carbohydrate metabolism, starch and sucrose metabolism, folding, sorting and degradation, translation, ribosome, and protein processing in endoplasmic reticulum (Fig. 6a). In the 9311 group, the pathways of carbohydrate metabolism, starch and sucrose metabolism, glycosyltransferases, glycolysis, and energy metabolism were enriched in the differentially ubiquitinated proteins (Fig. 6b); while there was no significantly enriched KEGG pathway in the GLA4 group. We further found the proteins that were common to both varieties were only significantly enriched in the starch and sucrose metabolism pathways (p = 0.04). The ubiquitination proteins involved in the starch and sucrose metabolism mainly include: sucrose hydrolysis (SUS, FK, UGPase), and starch synthesis (AGPase, GBSSI, BEI, BEIIb, PUL, PHO1), which are discussed below.

    Figure 6.  KEGG classification and enrichment analysis of differentially ubiquitinated proteins. (a) Number of differentially ubiquitinated proteins based on KEGG classification in 9311 and GLA4. (b) KEGG enrichment analysis of differentially ubiquitinated proteins in 9311.

    Although many reports have described specific examples of ubiquitination in rice defense responses[13,15,16], our knowledge on global changes in the developing endosperm ubiquitome under high-temperature stress is still lacking. In this study, a label-free quantitative proteomic analysis of ubiquitination was applied to examine the high-temperature induced ubiquitination change of two indica rice varieties (9311 and GLA4) with distinct starch quality. We identified many new lysine modification sites on proteins involved in various pathways, highlighting the complexity of the ubiquitination-mediated regulatory system in high-temperature stress responses in rice.

    Heat shock proteins accumulate under various stresses and play important roles in plant defenses against abiotic stresses[24,25]. Research has shown that a number of heat shock proteins were prominent in the rice ubiquitome network, of which OsHSP71.1, and OsHSP82A showed increased ubiquitination levels under chitin and flg22 treatment[15]. Here, seven lysine residues on five heat shock proteins possessed ubiquitination modification in rice endosperm. Three sites (BGIOSGA011420-K78, BGIOSGA026764-K99, BGIOSGA029594-K106) showed significant up-regulation in 9311 under high-temperature stress, while in GLA4, the ubiquitination level of BGIOSGA011420-K78 was down-regulated. This differential ubiquitination of heat-tolerant and heat-sensitive varieties provided a basis for studying the regulation of post-translational modification of heat shock proteins under high-temperature stress, despite the regulatory role of those heat shock proteins being still unclear.

    Transcription factors (TFs) play an essential role in the regulation of gene expression. A total of three transcription factors were identified in the ubiquitination dataset, of which two were NAC family members. As one of the largest plant-specific TF families, NAC is involved in the responses to abiotic and biotic stresses[26]. The ubiquitination modification of K173 in BGIOSGA018048 was specifically expressed in the 9311-H group, which may affect the stress resistance level in high-temperature environments. In addition, two sites K148 and K149 of ERF TF family member BGIOSGA024035 was downregulated in GLA4 under high-temperature stress. This differential ubiquitination was likely to affect the expression of related genes regulated by this transcription factor.

    The results of GO and KEGG enrichment analysis of differential ubiquitination proteins indicated that the sucrose and starch metabolic pathway was largely affected by ubiquitination regulation under high-temperature stress (Figs 5 & 6). The ubiquitination sites involved in sucrose and starch metabolism are listed in Table 1. To assess how high-temperature stress affects the crucial pathway, the significantly differential ubiquitination sites in 9311 and GLA4 were displayed in the heatmap of specific proteins (Fig. 7).

    Table 1.  Ubiquitination sites related to sucrose and starch metabolism in rice endosperm.
    Gene nameAnnotationProtein entryModification site(s)
    SUS1Sucrose synthase 1BGIOSGA010570K172, K177
    SUS2Sucrose synthase 2BGIOSGA021739K160, K165, K176, K804
    SUS3Sucrose synthase 3BGIOSGA026140K172, K177, K541, K544, K588
    FKFructokinaseBGIOSGA027875K143
    UGPaseUDP-glucose pyrophosphorylaseBGIOSGA031231K27, K150, K303, K306
    AGPS1ADP-glucose pyrophosphorylase small subunit 1BGIOSGA030039K94, K464, K484
    AGPS2ADP-glucose pyrophosphorylase small subunit 2BGIOSGA027135K106, K132, K385, K403, K406, K476, K496
    AGPL2ADP-glucose pyrophosphorylase large subunit 2BGIOSGA004052K41, K78, K134, K191, K227, K254, K316, K338, K394, K396, K463, K508, K513
    AGPL3ADP-glucose pyrophosphorylase large subunit 3BGIOSGA017490K509
    GBSSIGranule bound starch synthase IBGIOSGA022241K130, K173, K177, K181, K192, K258, K371, K381, K385, K399, K462, K517, K530, K549, K571, K575
    BEIStarch branching enzyme IBGIOSGA020506K103, K108, K122
    BEIIbStarch branching enzyme IIbBGIOSGA006344K134
    PULStarch debranching enzyme:PullulanaseBGIOSGA015875K230, K330, K431, K736, K884
    PHO1Plastidial phosphorylaseBGIOSGA009780K277, K445, K941
     | Show Table
    DownLoad: CSV
    Figure 7.  Sucrose and starch pathway at the ubiquitination levels in rice endosperm under high-temperature stress.

    In cereal endosperm, sucrose is the substrate for the biosynthesis of starch. The formation of glucose 1-phosphate (G1P, used in starch synthesis, see below) from sucrose requires a series of enzymes[27]. Here, we found that sucrose synthase 1 (SUS1), SUS2, SUS3, fructokinase (FK), and UDP-glucose pyrophosphorylase (UGPase) were ubiquitinated among all sample groups (Table 1, Fig. 7). In the ubiquitome of seedling and leaf in japonica rice, ubiquitination sites have been found in SUS1, SUS2, UGPase, and FK, which were related to sucrose hydrolysis[14,15]. SUS catalyzed the process of cleaving sucrose into UDP-glucose (UDPG) and fructose. Two ubiquitination sites, K172 and K177, were identified in SUS1 in rice endosperm, which were also found in rice leaves[14]. A total of four ubiquitination sites were identified in SUS2, two of which were also reported in rice seedling and leaf, indicating the conservation of the lysine residues in different rice tissues. It was noted that all four ubiquitination sites in SUS2 were upregulated in high-temperature environments, although the regulated sites of 9311 and GLA4 were different. In 9311, the ubiquitination levels of K160, K174, and K804 were increased, while GLA4 was only upregulated in K176. The ubiquitination sites K541, K544, and K588 in SUS3 were screened from developing rice seeds for the first time. In addition, SUS3 had two completely overlapping sites K172 and K177 with SUS1, and it was difficult to determine which enzymes the two sites belonged to. The ubiquitination levels of SUS3-K541 and SUS3-K544 in 9311 significantly increased in high-temperature environments, while there was no significant difference in the ubiquitination level of SUS3 in GLA4. Overall, the ubiquitination sites of SUS in rice endosperm were located in the functional domain except for SUS2-K804, reflecting the importance of ubiquitination regulation in SUS.

    UGPase catalyses the conversion of glucose 1-phosphate and UTP into UDPG[28]. Research has shown that the mutation of UGPase gene lead to chalky endosperm[29]. As shown in Table 1, four ubiquitination sites K27, K150, K303, and K306 were identified in rice endosperm, which were completely inconsistent with the seven ubiquitination sites in rice seedlings and two in leaves[14,15], reflecting the tissue specificity. We speculated that the UGPase with different modification sites may play different regulatory roles in metabolic pathways in different tissues. Under high-temperature stress, the ubiquitination level of UGPase-K27 was 8.1-fold up-regulated. Liao et al.[10] demonstrated that the expression of UDPase was down-regulated in both heat-tolerant and heat-sensitive rice lines under high temperature conditions, which could reasonably explain the significant up-regulation of UGPase-K27 ubiquitination level. The ubiquitination site K143 of FK was also reported in seedling tissues[15].

    The AGPase reaction represents the first committed step of starch biosynthesis[27]. A total of 22 lysine ubiquitination sites were identified in four AGPase subunits (AGPL2, AGPL3, AGPS1, AGPS2). AGPL2 had 13 ubiquitination sites, of which six were located in NTP_transferase domain, including K254, K338, K191, K134, K227, and K316. High-temperature stress resulted in an increase in the ubiquitination level of K254 in both 9311 and GLA4, and significant upregulation of K508 and K513 in 9311, as well as K191, K227, and K316 in GLA4. In contrast, AGPL2-K394 were significantly downregulated in GLA4. AGPL3 contained one ubiquitination site K509, and the modification level of AGPL3-K509 was up-regulated in high-temperature environments in 9311. AGPS1 had one specific ubiquitination site K464 and another two sites K94 and K484 that completely overlapped with AGPS2-K106 and AGPS2-K496, respectively. The modification levels of K464 and K484 significantly increased in high-temperature environments in 9311, and K94 was significantly up-regulated in both varieties. There were seven ubiquitination sites in AGPS2 in rice endosperm, which were different with the sites found in rice leaves[14]. In addition to the two sites that overlapped with AGPS1, AGPS2 had another two ubiquitination sites (K406 and K132) that upregulated in high-temperature environments.

    Amylose content is one of the key determinants that strongly influence rice grain quality[30]. The biosynthesis of amylose requires the catalytic effect of granule-bound starch synthase I (GBSSI)[30,31]. Here, a total of 16 ubiquitination sites were identified in GBSSI (Table 1, Fig. 8a). Among these ubiquitination sites, six lysine residues (K130, K173, K177, K181, K192, K258) were located in glycosyltransferase 5 (GT5) domain, and three sites (K399, K462, K517) were located in GT1 domain (Fig. 6), indicating the important role of ubiquitination regulation of GBSSI. Under high-temperature stress, the ubiquitination levels of six sites (K130, K177, K399, K381, K385, K549) increased in two indica rice varieties, while one sites (K258) showed downregulation in 9311 (Fig. 8a). Numerous studies had described that the amylose content was reduced under high-temperature stress in rice[5,7], which might be due to the degradation of GBSSI proteins caused by the increased significantly up-regulated ubiquitination sites. These ubiquitination sites identified in rice GBSSI with significant differences under high-temperature stress were expected to become a new breakthrough point for the improvement of starch quality.

    Figure 8.  Structure of GBSSI. (a) Domain structure of GBSSI and ubiquitination sites with significant differences in response to high-temperature stress. (b) 3D model of GBSSI and the relationship between ubiquitination sites K462 and ADP, SO4 (salt bridge or hydrogen bond).

    To further determine the regulatory role of the ubiquitination sites in GBSSI, SWISS-MODEL was used to predict 3D structural model. As shown in Fig. 8b, GBSSI had three SO4 (sulfate ions) and one ADP ligand. These ligands interact with GBSSI through hydrogen bonds and salt bridges. Three sites, K447, R458, and K462, were associated with SO4 through salt bridges, while G100, N265, Q412, K462, and Q493 interact with the hydrogen bonds of ADP in GBSSI[32,33]. Based on this finding, it can be reasonably inferred that the K462 site with ubiquitination modification located in the GT1 domain played an important role in the interaction between GBSSI, SO4, and ADP. An in-depth investigation was necessary to gain a more comprehensive understanding of the regulatory function of ubiquitination modification at GBSSI-K462, although there was no significant difference in the ubiquitination level under high-temperature stress.

    Amylopectin, the major component of starch, is synthesized by the coordinated action of multiple enzymes including soluble starch synthase (SSs), starch branching enzyme (BEs), starch debranching enzyme (DBEs), and phosphorylases (PHOs or Phos) with ADPG as a substrate. In this study, ubiquitination sites were detected in BEs, DBEs, and Phos.

    BEs, covering two isoforms, BEI and BEII, are responsible for catalyzing the formation of α-1,6-glucosidic linkages of amylopectin[34]. There were three ubiquitination sites (K103, K108, and K122) identified in BEI (Fig. 9a). K122 was the first amino acid in the carbohydrate-binding module 48 (CBM48) domain. Sequence alignment analyses of BEs from eight plants revealed K122 was conserved among all plants' BEI (Fig. 9a), suggesting a high probability of the functional effects of ubiquitination modification at this site. In high-temperature environments, ubiquitination levels of K108 and K122 were significantly up-regulated in 9311, while no significantly regulated ubiquitination sites of BEI were observed in GLA4. Only one ubiquitination site, K134, was found in BEIIb (Fig. 9a). The ubiquitination levels showed a slightly upward trend with no significant differences in high-temperature environments in both varieties. These changes could be one of the reasons for increased gelatinization temperature and relative crystallinity of rice starch in response to high-temperature[5].

    Figure 9.  Domain structure of (a) BEs, (b) PUL and (c) Pho1 as well as their ubiquitination sites with significant differences in response to heat stress. Residues in red indicate the ubiquitination site. Non-ubiquitinated residues are shown in dark grey.

    DBEs consists of isoamylase (ISA) and pullulanase (PUL) with catalytic function for hydrolyzing α-1,6-glucosic linkages[35]. In the present study, we found that only PUL was ubiquitinated in rice endosperm (Fig. 9b). Among five ubiquitination sites (K230, K330, K432, K736, and K884) identified in PUL, K230 was located in the PULN2 domain, while K330 was in the CBM48 domain. Under high-temperature stress, K330 showed completely opposite regulatory trends in two cultivars. In addition, the ubiquitination level of K884, located in the DUF3372 domain, was significantly up-regulated in 9311. Previous study has reported that the expression of PUL was significantly up-regulated in 9311 under high-temperature stress, while GLA4 showed down-regulation in PUL abundance[11]. Consequently, there might be two possible functions of these ubiquitination sites. One possibility is that ubiquitination sites were unrelated to protein degradation; instead, they regulated the biosynthesis of amylopectin by affecting other functions of the protein. Secondly, ubiquitination sites were associated with protein degradation, and the levels of ubiquitination modification were based on protein abundance, resulting in a completely consistent regulation of ubiquitination modification and protein abundance under high-temperature stress.

    PHOs, including two types, Pho1/PHO1 and Pho2/PHO2, are responsible for the transfer of glucosyl units from Glc-1-P to the non-reducing end of a-1,4-linked glucan chains[36]. Pho1 is a temperature-dependent enzyme and considered crucial not only during the maturation of amylopectin but also in the initiation process of starch synthesis[37,38]. The three ubiquitination sites (K277, K445, K941) identified in Pho1 were located in two phosphorylase domains. We found that two sites, Pho1-K277 and Pho1-K445, were only ubiquitinated in high-temperature environments in 9311 and GLA4, respectively. Pang et al.[11] has demonstrated that the protein abundance of Pho1 decreased under high-temperature stress, especially in GLA4. Satoh et al.[38] reported that the functional activity of Pho1 was weakened under conditions of high temperature and its function might be complement by one or more other factors. Hence, these ubiquitination modifications that specifically occurred in high-temperature environments might be related to the degradation of Pho1 proteins.

    As a factory for protein synthesis in cells, the ribosome is an extremely crucial structure in the cell[39]. It has been proven that multiple ribosomal subunits were abundantly ubiquitinated in Arabidopsis and wheat[22]. In the present study, 57 ubiquitination sites involving 33 ribosome subunits were identified in 40S and 60S ribosome complexes in rice. Under high-temperature stress, the ubiquitination levels of some sites were significantly upregulated or downregulated, implying that ubiquitination of ribosomal proteins is likely to be an important regulatory mechanism in high-temperature response in rice endosperm. The results of GO and KEGG enrichment analysis indicated that the ribosome system was one of the most active systems for ubiquitination regulation under high-temperature stress. We speculated that the ubiquitin-proteasome system might be involved in the removal of subunits or entire ribosomes that were improperly folded in high-temperature environments. As shown in Fig. 10, the S10e, L18Ae, S27Ae, L9e, S3e, S28e, S20e, and S2e subunits were significantly up-regulated in 9311, while L13e subunits showed a completely opposite regulatory trend at the ubiquitination sites K81 and K88. In GLA4, the ubiquitination levels of S10e, S27Ae, L10Ae, L9e, S3e, S2e, and L4e showed a significant increase, while the ubiquitination level of L17e was significantly down-regulated under high-temperature stress. A total of seven ubiquitination sites involving S10e, S27Ae, L9e, S3e, and S2e subunits were jointly up-regulated in both two varieties. These sites might be related to the degradation of improperly folded ribosome subunits under high-temperature stress, while other ubiquitination sites with variety specificity might be associated with ribosomal function.

    Figure 10.  Ribosome system at the ubiquitination levels in rice endosperm under high-temperature stress. Grey shadings represent ubiquitinated proteins with no significant differences under heat stress. Red and blue shadings indicate up-regulated and down-regulated ubiquitinated proteins, respectively. Orange shading displays a combination of up- and down-regulated ubiquitinated sites in the same ubiquitinated protein.

    In conclusion, this study provides the first comprehensive view of the ubiquitome in rice developing endosperm, and demonstrated that ubiquitination has diverse functions in the high-temperature response of rice endosperm by modulating various cellular processes, especially the sucrose and starch metabolism. Comparative analysis of the temperature-induced ubiquitination status revealed some similarities and more interesting differences between 9311 and GLA4. These differences might be the reason for the different qualities formation of the two indica rice varieties, which could provide potential genetic resources for the improvement of the heat resistance in rice. Considering the diversity of ubiquitination modification, it is worthwhile to further validate and explore the function and regulatory mechanism of the key targets and key pathways. The findings provide valuable insights into the role of ubiquitination in response to high-temperature stress and lay a foundation for further functional analysis of lysine ubiquitination in rice.

    The authors confirm contribution to the paper as follows: study conception and design: Bao J, Pang Y; data collection: Pang Y; analysis and interpretation of results: Pang Y; draft manuscript preparation: Ying Y; Revised manuscript preparation: Ying Y, Pang Y, Bao J. All authors reviewed the results and approved the final version of the manuscript.

    The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

    This work was financially supported by the AgroST Project (NK2022050102) and Zhejiang Provincial Natural Science Foundation (LZ21C130003).

  • The authors declare that they have no conflict of interest.

  • Supplemental Table S1 Primer sequences for real-time fluorescence quantitative PCR.
    Supplemental Table S2 Statistics on clean reads mapping against reference genome.
    Supplemental Table S3 Expression of miRNAs in each sample (TPM).
    Supplemental Table S4 The miRNA family analysis.
    Supplemental Table S5 All DE-miRNA.
    Supplemental Table S6 Analysis of the DEMs co-expression trends.
    Supplemental Table S7 Targeting relationship between miRNA and target genes.
    Supplemental Fig. S1 BR vs MR's miRNAs and target genes expression levels and targeting relationships.
    Supplemental Fig. S2 BS vs MS′s miRNAs and target genes expression levels and targeting relationships.
  • [1]

    Rogers K, Chen X. 2013. Biogenesis, turnover, and mode of action of plant microRNAs. The Plant Cell 25:2383−99

    doi: 10.1105/tpc.113.113159

    CrossRef   Google Scholar

    [2]

    Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. 2002. MicroRNAs in plants. Genes & Development 16:1616−26

    doi: 10.1101/gad.1004402

    CrossRef   Google Scholar

    [3]

    Han J, Fang J, Wang C, Yin Y, Sun X, et al. 2014. Grapevine microRNAs responsive to exogenous gibberellin. BMC Genomics 15:111

    doi: 10.1186/1471-2164-15-111

    CrossRef   Google Scholar

    [4]

    Li Y. 2022. Identification of cucumber miR160 gene family and functional analysis of Csa-miR160d. Thesis. Henan University of Science and Technology.

    [5]

    Zhang F. 2020. Molecular mechanisms of DNA methylation and microRNA regulation of chrysanthemum inflorescence development. Thesis. Beijing Forestry University.

    [6]

    Wang M, Wang Q, Zhang B. 2019. Identification and characterization of miRNAs in apple rootstock 'M26' and their potential roles in regulating grafting process. Gene 684:108−18

    Google Scholar

    [7]

    Zhang L. 2022. Comparison of physiological responses of two wheat cultivars to saline-alkali stress and screening of differential microRNA. Huaibei: Huaibei Normal University. https://doi.org/10.27699/d.cnki.ghbmt.2022.000149

    [8]

    Li N, Wang J, Wang B, Da Q, Huang S, et al. 2021. Research progress on microRNA regulation of growth, development, and stress responses in tomato. Xinjiang Agricultural Sciences 58:474−82

    Google Scholar

    [9]

    Zhu H, Xia R, Zhao B, An Y, Dardick CD, et al. 2012. Unique expression, processing regulation, and regulatory network of peach (Prunus persica) miRNAs. BMC Plant Biology 12:149

    doi: 10.1186/1471-2229-12-149

    CrossRef   Google Scholar

    [10]

    Zhang J, Ai X, Guo W, Peng S, Deng X, et al. 2012. Identification of miRNAs and their target genes using Deep sequencing and degradome analysis in trifoliate orange [Poncirus trifoliate (L.) Raf]. Molecular Biotechnology 51:44−57

    doi: 10.1007/s12033-011-9439-x

    CrossRef   Google Scholar

    [11]

    Niu Q, Qian M, Liu G, Yang F, Teng Y. 2013. A genome-wide identification and characterization of mircoRNAs and their targets in 'Suli' pear (Pyrus pyrifolia white pear group). Planta 238:1095−112

    doi: 10.1007/s00425-013-1954-5

    CrossRef   Google Scholar

    [12]

    Ye K, Chen Y, Hu X, Guo J. 2013. Computational identification of microRNAs and their targets in apple. Genes & Genomics 35:377−85

    doi: 10.1007/s13258-013-0070-z

    CrossRef   Google Scholar

    [13]

    Wang M, Li T, Wu Y, Song S, Bai T, et al. 2021. Genome-wide identification of microRNAs involved in the regulation of fruit ripening in apple (Malus domestica). Scientia Horticulturae 289:110416

    doi: 10.1016/j.scienta.2021.110416

    CrossRef   Google Scholar

    [14]

    Varkonyi-Gasic E, Gould N, Sandanayaka M, Sutherland P, MacDiarmid RM. 2010. Characterisation of microRNAs from apple (Malus domestica 'Royal Gala') vascular tissue and phloem sap. BMC Plant Biology 10:159

    doi: 10.1186/1471-2229-10-159

    CrossRef   Google Scholar

    [15]

    Yu X, Hou Y, Chen W, Wang S, Wang P, et al. 2017. Malus hupehensis miR168 targets to ARGONAUTE1 and contributes to the resistance against Botryosphaeria dothidea infection by altering defense responses. Plant and Cell Physiology 58:1541−57

    doi: 10.1093/pcp/pcx080

    CrossRef   Google Scholar

    [16]

    Xing L, Zhang D, Li Y, Zhao C, Zhang S, et al. 2014. Genome-wide identification of vegetative phase transition-associated microRNAs and target predictions using degradome sequencing in Malus hupehensis. BMC Genomics 15:1125

    doi: 10.1186/1471-2164-15-1125

    CrossRef   Google Scholar

    [17]

    Corbsier L, Vincent C, Jand S, Fornara F, Fan Q, et al. 2007. FT protein movement contributes to long-distance signaling floral induction of Arabidopsis. Science 316:1030−33

    doi: 10.1126/science.1141752

    CrossRef   Google Scholar

    [18]

    Guo X. 2017. The involvement of sRNA in apple (Malus domestica) floral transition. Beijing: China Agricultural University.

    [19]

    Guo X, Ma Z, Zhang Z, Cheng L, Zhang X, et al. 2017. Small RNA-sequencing links physiological changes and RdDM process to vegetative-to-floral transition in apple. Frontiers in Plant Science 8:873

    doi: 10.3389/fpls.2017.00873

    CrossRef   Google Scholar

    [20]

    Song C, Zhang D, Zheng L, Zhang J, Zhang B, et al. 2017. miRNA and degradome sequencing reveal miRNA and their target genes that may mediate shoot growth in spur type mutant "Yanfu 6". Frontiers in Plant Science 8:441

    doi: 10.3389/fpls.2017.00441

    CrossRef   Google Scholar

    [21]

    Yu X, Du B, Gao Z, Zhang S, Tu X, et al. 2014. Apple ring rot-responsive putative microRNAs revealed by high-throughput sequencing in Malus × domestica Borkh. Molecular Biology Reports 41:5273−86

    doi: 10.1007/s11033-014-3399-8

    CrossRef   Google Scholar

    [22]

    Huang YC, Tsay TT. 2017. The roles of rootstock and scion in grafting. Horticulture Plant Journal 3:193−201

    Google Scholar

    [23]

    Davoudi M, Song M, Zhang M, Chen J, Lou Q. 2022. Long-distance control of pumpkin rootstock over cucumber scion under drought stress as revealed by transcriptome sequencing and mobile mRNAs identifications. Horticulture Research 9:uhab033

    doi: https://doi.org/10.1093/hr/uhab033

    CrossRef   Google Scholar

    [24]

    Kviklys D, Samuolienė G. 2020. Relationships among the rootstock, crop load, and sugar hormone signaling of apple tree, and their effects on biennial bearing. Frontiers in Plant Science 11:1213

    doi: 10.3389/fpls.2020.01213

    CrossRef   Google Scholar

    [25]

    Pant BD, Buhtz A, Kehr J, Scheible WR. 2008. MicroRNA399 is a long-distance signal for the regulation of plant phosphate homeostasis. The Plant Journal 53:731−38

    doi: 10.1111/j.1365-313X.2007.03363.x

    CrossRef   Google Scholar

    [26]

    Zhang F, Zhong H, Zhou X, Pan M, Xu J, et al. 2022. Grafting with rootstocks promotes phenolic compound accumulation in grape berry skin during development based on integrative multi-omics analysis. Horticulture Research 9:uhac055

    doi: 10.1093/hr/uhac055

    CrossRef   Google Scholar

    [27]

    Thieme CJ, Rojas-triana M, Stecyk E, Schudoma C, Zhang W, et al. 2015. Endogenous Arabidopsis messenger RNAs transported to distant tissues. Nature Plants 1:15025

    doi: 10.1038/nplants.2015.25

    CrossRef   Google Scholar

    [28]

    Wang W. 2017. Genome-wide analysis of mRNAs, miRNAs and ALDH gene family in soybean response to drought stress. Nanjing: Nanjing Agricultural University.

    [29]

    Buhtz A, Springer F, Chappell L, Baulcombe DC, Kehr J. 2008. Identification and characterization of small RNAs from the phloem of Brassica napus. The Plant Journal 53:739−49

    doi: 10.1111/j.1365-313x.2007.03368.x

    CrossRef   Google Scholar

    [30]

    Tzarfati R. 2013. Graft-induced changes in microRNA expression patterns in citrus leaf petioles. The Open Plant Science Journal 7:17−23

    Google Scholar

    [31]

    Xu Y. 2018. Effect of grafting and stress on miRNAs expression in Citrus. Chongqing: Southwest University.

    [32]

    Liu N, Yang J, Guo S, Xu Y, Zhang M. 2013. Genome-wide identification and comparative analysis of conserved and novel microRNAs in grafted watermelon by high-throughput sequencing. PLoS ONE 8:e57359

    doi: 10.1371/journal.pone.0057359

    CrossRef   Google Scholar

    [33]

    An N, Fan S, Yang Y, Chen X, Dong F, et al. 2018. Identification and characterization of miRNAs in self-rooted and grafted Malus reveals critical networks associated with flowering. International Journal of Molecular Sciences 19:2384

    doi: 10.3390/ijms19082384

    CrossRef   Google Scholar

    [34]

    Zhang Q, Yang X, Li F, Deng Y. 2022. Advances in miRNA-mediated growth and development regulation in horticultural crops. Acta Horticulturae Sinica 49:1145−61

    Google Scholar

    [35]

    Li Q, Gao Y, Wang K, Feng J, Sun S, et al. 2023. Transcriptome analysis of the effects of grafting interstocks on apple rootstocks and scions. International Journal of Molecular Sciences 24:807

    doi: 10.3390/ijms24010807

    CrossRef   Google Scholar

    [36]

    Livak KJ, Schmittgen TD. 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCᴛ method. Methods 25:402−8

    doi: 10.1006/meth.2001.1262

    CrossRef   Google Scholar

    [37]

    Guo H. 2020. Research progress on interactive mechanism between rootstock and scion of grafted plants. Guizhou Agricultural Sciences 48:35−44

    Google Scholar

    [38]

    Pasquinelli AE. 2012. MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship. Nature Reviews Genetics 13:271−82

    doi: 10.1038/nrg3162

    CrossRef   Google Scholar

    [39]

    Wang Y, Wang X, Wu W. 2018. Identification and characterization of microRNAs and their target genes in apple rootstock under waterlogging stress. Gene 651:148−59

    Google Scholar

    [40]

    Tahir MM, Zhang X, Shah K, Hayat F, Li S, et al. 2021. Nitrate application affects root morphology by altering hormonal status and gene expression patterns in B9 apple rootstock nursery plants. Fruit Research 1:14

    doi: 10.48130/frures-2021-0014

    CrossRef   Google Scholar

    [41]

    Wang M, Wang Q, Zhang B. 2013. Response of miRNAs and their targets to salt and drought stresses in cotton (Gossypium hirsutum L.). Gene 530:26−32

    doi: 10.1016/j.gene.2013.08.009

    CrossRef   Google Scholar

    [42]

    Li X, Mao J. 2017. Identification and functional analysis of miRNAs in apple rootstock Malus hupehensis. Plant Molecular Biology Reporter 35:363−74

    Google Scholar

    [43]

    Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, et al. 2007. Target mimicry provides a new mechanism for regulation of microRNA activity. Nature Genetics 39:1033−37

    doi: 10.1038/ng2079

    CrossRef   Google Scholar

    [44]

    Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, et al. 2005. Specific effects of microRNAs on the plant transcriptome. Developmental Cell 8:517−27

    doi: 10.1016/j.devcel.2005.01.018

    CrossRef   Google Scholar

    [45]

    Wang W, Wang S, Li M, Hou L. 2017. Expression analysis and miRNA prediction of wax synthetic genes in cucumber grafted onto different rootstock. Journal of Nuclear Agricultural Sciences 31:1896−903

    Google Scholar

    [46]

    He Q, Zhu S, Zhang B. 2014. MicroRNA-target gene responses to lead-induced stress in cotton (Gossypium hirsutum L.). Functional & Integrative Genomics 14:507−15

    doi: 10.1007/s10142-014-0378-z

    CrossRef   Google Scholar

    [47]

    Sun P, Tahir MM, Lu X, Liu Z, Zhang X, et al. 2022. Comparison of leaf morphological, anatomical, and photosynthetic responses to drought stress among eight apple rootstocks. Fruit Research 2:20

    doi: 10.48130/frures-2022-0020

    CrossRef   Google Scholar

    [48]

    Voinnet O. 2009. Origin, biogenesis, and activity of plant MicroRNAs. Cell 136:669−87

    doi: 10.1016/j.cell.2009.01.046

    CrossRef   Google Scholar

    [49]

    Jones-Rhoades MW, Bartel DP. 2004. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Molecular Cell 14:787−99

    doi: 10.1016/j.molcel.2004.05.027

    CrossRef   Google Scholar

    [50]

    Xing L, Zhang D, Zhao C, Li Y, Ma J, et al. 2016. Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.). Plant Biotechnology Journal 14:749−70

    doi: 10.1111/pbi.12425

    CrossRef   Google Scholar

    [51]

    Chen T, Chen X, Zhang S, Zhu J, Tang B, et al. 2021. The genome sequence archive family: toward explosive data growth and diverse data types. Genomics, Proteomics & Bioinformatics 19:578−83

    doi: 10.1016/j.gpb.2021.08.001

    CrossRef   Google Scholar

    [52]

    CNCB-NGDC Members and Partners. 2022. Database Resources of the National Genomics Data Center, China National Center for Bioinformation in 2022. Nucleic Acids Research 50:D27−D38

    doi: 10.1093/nar/gkab951

    CrossRef   Google Scholar

  • Cite this article

    Li Q, Gao Y, Wang K, Sun S, Lu X, et al. 2023. Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock. Fruit Research 3:34 doi: 10.48130/FruRes-2023-0034
    Li Q, Gao Y, Wang K, Sun S, Lu X, et al. 2023. Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock. Fruit Research 3:34 doi: 10.48130/FruRes-2023-0034

Figures(9)  /  Tables(1)

Article Metrics

Article views(4343) PDF downloads(578)

ARTICLE   Open Access    

Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock

Fruit Research  3 Article number: 34  (2023)  |  Cite this article

Abstract: Apple (Malus domestica) is an important multi-functional horticultural crop. However, thus far, there have been few reports concerning the molecular mechanism by which interstock regulates rootstock and scion growth. MicroRNAs (miRNAs) are a type of post-transcriptional regulator that regulates apple growth and development, fruit quality, hormone signal transduction, and stress response. Interstock grafting resulted in differences in miRNA distribution between apple tissues. The regulatory roles of miRNAs in the grafting of apples are not yet clear. In this study, 15 libraries were constructed using the phloem of an apple grafted on two types of interstock and rootstock by high-throughput sequencing. A total of 281 miRNAs were catalogued into 80 families, and 159 novel ones were identified. Compared with the control, the grafting combination with the interstock resulted in 79 differentially expressed miRNAs (DEMs) in the scion, with 36 up-regulated and 43 down-regulated miRNAs, and 57 DEMs in the rootstock, including 22 up-regulated and 35 down-regulated miRNAs. Grafting with the dwarfing interstock led to DEMs in the whole plant, including a decrease in the expression of mdm-miR156 in the scion and mdm-miR172 in the rootstock. Predictive analysis of the DEMs and their target genes suggested that miRNAs mediate scion growth through several aspects, such as plant hormone synthesis and signal transduction, and nutrient absorption and balance. Through combined miRNA and mRNA analysis, the dwarfing effect of the interstock may affect the expression of genes related to the starch and sucrose metabolism pathways in the rootstock and the plant hormone signal transduction pathway in the scion. The available evidence facilitates a better understanding of the role of miRNAs in the response of apples after grafting interstock. One possible reason for the stunting of apple trees and the promotion of early flowering caused by the dwarfing interstock is the decrease in expression levels of mdm-miR156 and mdm-miR172.

    • MicroRNAs (miRNAs) are an extensive class of 22-nucleotide noncoding RNAs thought to regulate gene expression in post-transcriptional regulation[1]. MiRNAs in plants were first found in Arabidopsis thaliana, and they are involved in many basic biological processes, such as plant growth, morphological formation, fruit development, and biological abiotic stress response[2,3] MiRNAs play an important role in the growth and development of plants. A study by Li found that in regulating plant root development, miR160 downregulates the expression of the auxin response factor (ARF) family genes, thereby controlling plant root differentiation and growth[4]. During the development of floral organs, miR172 promotes flower opening and maturation by inhibiting the expression of the Apetala2 (AP2) gene[5]. MiRNAs also play an important role in plant stress responses. In plant drought stress, miR169 enhances plant drought resistance by downregulating the expression of the nuclear factor Y subunit A (NF-YA) family genes[6]. Under salt stress, miR393 enhances plant salt tolerance by regulating the plant hormone response pathway[7]. In addition, miRNAs are involved in regulating plant virus resistance, climate change, and environmental pollution[8].

      MiRNAs in grape, peach, apple, orange, and pear have been studied genome-wide since the development of high-throughput sequencing technology and the completion of fruit tree genome sequencing[3,912]. So far, several hundred plant miRNA family species have been found. Wang et al. identified 154 miRNAs from apples belonging to 26 families, some of which were generated by genome duplication during evolution[13]. Zhu et al. identified 47 peach-specific and 47 known miRNAs that targeted 80 genes, among which 3 miRNAs specifically targeted 49 MYBs, of which 19 regulated phenylpropanoid metabolism, a key pathway related to stone hardening and fruit colour development[9]. MiRNAs in apples show tissue-enriched expression patterns. For example, the conserved miRNAs miR156, miR159, miR160, miR167, and miR172 have different expression patterns in different tissues and developmental stages of apples[14]. miRNA172 is a key factor controlling apple fruit size[15].

      High expression of miR156 and low expression of miR172 may shorten the juvenile phase, and miR160 and miR393 can promote the transformation of fruit trees at the same time by regulating genes related to auxin signal transduction[16]. Corbsier et al. found that the flowering mechanism of some grafted plants was regulated by miRNAs[17]. Guo found significant differences in the expression of 33 known miRNAs and six newly annotated miRNAs between apple leaf buds and flower buds, mostly showing a negative correlation with the process of flower bud formation[18]. After shoot bending treatment in apples, the differentially expressed miRNAs and their target genes were associated with metabolic pathways such as auxin, abscisic acid, gibberellins, and flowering, and they also participated in various biological processes such as cellular biosynthesis and metabolism[19]. Song et al. found that miR159, miR167, miR396, and their target genes could regulate cell differentiation and internode length, while miR164, miR166, miR171, and their target genes could regulate the growth of the stem apical meristem[20]. Yu et al. found that miRNAs also participate in apple disease resistance[21].

      Grafting is a kind of asexual propagation method that can maintain the inherent good properties of scions and is widely used in the production of vegetables, ornamental horticultural crops, fruit trees, and other horticultural crops[22, 23]. Plant grafting has been used for millennia to reduce juvenility and confer biotic and abiotic stress tolerance[24]. Grafting can form a complete vascular bundle between the rootstock and scion, and during the growth and development of grafted plants, large molecules such as proteins, mRNA, and miRNA in the phloem selectively move from the companion cells through plasmodesmata to sieve tubes and are transported over long distances between the rootstock and scion, ultimately reaching the corresponding target cells to regulate the expression of related genes and thereby regulate the development of plant organs. Pant et al. used Arabidopsis thaliana overexpressing miR399 as the scion grafted onto wild-type Arabidopsis and found that the roots of wild-type Arabidopsis accumulated a large amount of mature miR399, while pri-miR399 was barely detected in its roots, indicating that miR399 could move long distances through the phloem[25]. Zhang et al. found in their study that grafting with rootstocks promotes phenolic compound accumulation in grape berry skin during development based on integrative multi-omics analysis[26]. Thieme et al. found that a large number of miRNAs move in specific tissues to regulate plant tissue development[27]. As a signalling molecule, miRNA controls traits such as flower and leaf shape in grafted plants. Wang found differential expression of miR164, miR168, and miR390 in the soybean after drought stress[28]. Buhtz et al. also found differentially expressed miRNAs after grafting in Brassica napus. MiR398, miR399, and miR395 in B. napus may be transmitted from the rootstock phloem to the scion and affect the scion's response to mineral element deficiency[29].

      Compared to vegetables and other horticultural crops, there are fewer studies on miRNAs involved in the grafting of fruit trees due to their long life cycle and complex tree structure. In recent years, there has been some progress in miRNA research on the interaction between rootstock and scion in fruit trees, and it has been found that grafting significantly affects the expression levels of many miRNAs. In grafted Citrus seedlings, the expression of miR156, miR157, and miR894 in the leaf petioles is significantly downregulated compared to that in seedlings grown from seeds[30]. Xu found that miRNAs associated with regulating plant growth and development, stress response, and hormone signal transduction were upregulated in grafted Citrus scions, while they were downregulated in the rootstocks[31]. Liu et al. suggested that changes in miRNA expression in grafted watermelon leaves may be one of the reasons why grafting affects plant growth, development, and stress response[32]. Grafting can induce differential expression of miRNAs in Citrus leaf petioles, and compared with seedlings grown from seeds, the expression of miR156, miR157, and miR894 is significantly downregulated in grafted seedlings[30]. An et al. found that grafting altered the expression patterns of miR156 and miR172, which are related to flowering. Specific miRNAs differentially expressed during grafting may be involved in regulating growth, development, and metabolism in fruit trees after grafting[33]. In the grafting process of apples 'Gala' and 'M9', the expression levels of miR172, miR396, miR399, and miR167 changed, which may be related to plant auxin and carbon-nitrogen metabolism processes[18].

      There are significant differences in the expression levels of miRNAs in apple fruits grafted with different rootstocks. In 'Fuji' apple fruits, the expression levels of miR156, miR172, miR390, miR397, and miR408 in fruits grafted with dwarfing rootstocks 'M9' and 'B9' were higher, while the expression levels of miR164 and miR167 in fruits grafted with 'MM111' rootstock were higher[34]. Many miRNAs have been discovered in horticultural crops, but the target genes and functional characteristics of most graft-responsive miRNAs are still unclear.

    • The following apple grafting combinations were used in the study: the dwarfing interstock combination M ('Gala'/'Mac 9'/Malus baccata (L.) Borkh) and the standard rootstock combination B ('Gala'/Malus baccata (L.) Borkh). The combinations used in this study were provided by the Research Institute of Pomology of the Chinese Academy of Agricultural Sciences (CAAS), Xingcheng, Liaoning, China (120° 44′ E, 40° 37′ N). For each sample, at least three plants were pooled, and three independent biological replicates were used.

    • On June 21, 2016, six healthy trees without pests or diseases were selected for B and M. Each sample obtained 2 cm long and wide bark slices, and two samples (BS and BR) in B, three samples (MS, MI, and MR) were collected in M, with three biological replicates for each tree. Before transferring the bark samples to a −80 °C ultra-low-temperature freezer, they were briefly stored in liquid nitrogen. First, total RNA needs to be extracted from the bark samples. Then, the extracted RNA samples were sent to Biomarker Technologies for small RNA sequencing. After passing quality control, the library was constructed. To obtain RNA with a length of 18–30 nt, PAGE gel electrophoresis was used for gel cutting and separation. Subsequently, the Clean Reads, which were trimmed for adapter sequences, low-quality bases, and contaminants, were aligned to the Rfam13 database for annotation.

      The expression of miRNAs in each sample was calculated using the TPM algorithm and normalised. The equation for TPM normalisation is as follows:

      TPM=Readcount1,000,000MappedReads

      In the equation, readcount represents the number of reads that correspond to the miRNA, and Mapped Reads represent the number of reads that are mapped to all miRNAs. Known and novel miRNAs were identified using the miRDeep2 software.

    • After normalization of expression levels using the TPM algorithm, the influence of sequencing depth differences between different samples can be eliminated, enabling a scientific and accurate comparison of differential expression of miRNAs among samples. To investigate the differential expression of miRNAs between two biological conditions, DESeq was used to perform differential expression analysis between sample groups. The fold change (FC) and false discovery rate (FDR) were used as filtering criteria. In the original hypothesis test, the recognized Benjamini Hochberg correction method was used to correct the significance p-values, and the FDR was ultimately used as the key indicator for filtering differentially expressed miRNAs. We compared the miRNAs of B and M to determine whether there were differences under the intermediate rootstock grafting. In this study, the threshold for defining differentially expressed miRNAs was set at |log2(FC)| ≥ 0.58 and p-value ≤ 0.05.

    • For the sequencing data of this study, we used the Target Finder software to predict the target genes of miRNAs, using known miRNAs and newly predicted miRNAs with corresponding gene sequence information of the species. After completing the prediction of target genes, we used the BLAST software to compare the predicted target genes with the nr, Swiss Prot, GO, COG, KEGG, KOG, and Pfam databases to obtain annotation information for the target genes.

    • DESeq2_edgeR p value = 0.05, FC = 1.5. Analysis of the different patterns of changes in miRNA expression abundance in two different parts of the bark, and grouping mRNAs with the same expression trend into a dataset, creating an expression pattern graph for the dataset. The distance measurement method used is Euclidean distance, and the clustering method is K-means clustering. Finally, the expression pattern will be divided into K categories.

    • Based on the miRNA sequencing data in this study and the transcriptome sequencing in early research using differentially expressed miRNAs as the screening criteria, we searched for mRNAs regulated by these differentially expressed miRNAs, and focused on analyzing the negative regulatory relationship between miRNAs and mRNAs[35]. At the same time, functional enrichment analysis was performed on the resulting mRNAs. The targeting relationship between miRNAs and mRNAs can be visualized using a Sankey diagram.

    • Total RNA was extracted from phloem of rootstock, interstock, and scion of 'Gala'/'Mac 9'/Malus baccata (L.) Borkh and 'Gala'/Malus baccata (L.) Borkh using the Polysaccharide and Polyphenol RNA Extraction Kit (Tiangen, China). Real-time quantitative PCR (qRT-PCR) was carried out to validate the levels of DEGs from phloem of B and M combinations according to the instructions of the Prime ScriptTM RT reagent Kit with gDNA Eraser (TaKaRa, Code No. RR036A), 1,000 ng of the total RNA was reverse-transcribed with random primers. The reaction system included 1,000 ng of RNA, 2 µL of 5 × PrimeScript RT Master Mix (Perfect Real Time), and RNase-free dH2O was added to obtain a 20 µL volume; the reaction conditions were 37 °C for 15 min and 85 °C for 5 s. Four target genes with miRNA were selected from the transcriptome data. Primers were designed with primer premier 5.0 software and relevant primers were synthesized by Sangon Biotech (Shanghai, China). The qRT-PCR was performed using TB Green Premix Ex Taq II (TaKaRa, Code No. RR820B) in a Light Cycler 96 quantitative PCR detection system (Roche, Switzerland). The qRT-PCR aliquot contained 10 µL TB Green Premix Ex Taq II (Tli RNaseH Plus) (2 ×), 6 µL of ddH2O, 2 µL of cDNA, 0.8 µL of each of the forward and reverse primers (10 μM), and 0.4 µL of ROX Reference Day II (50 ×). The reaction conditions included preincubation at 95 °C for 30 s, followed by 40 cycles at 95 °C for 7 s, 60 °C for 30 s, and 72 °C for 30 s, with melting (65~97 °C, +1 °C/cycle; holding time: 1 s). The q-Actin gene of malus was used as the reference to correct expression levels of related genes. All the PCR assays were performed with three biological replicates. The relative expression levels were calculated with the 2−ΔΔCᴛ method[36].

      The reverse transcription of miRNA was usually performed using the stem-loop method, and miRNA primers were designed using Shanghai Biotech's online primer design tool (https://www.sangon.com/class_mirna_stemloop.html) which was a relatively complex process. First, total RNA was extracted, and then miRNA first-strand cDNA was synthesised. The reaction system included 10 l of miRNA L-RT Solution mix, 1.5 l of miRNA L-RT Enzyme mix, 2 l of total RNA, 1 l of stem-loop primer (10 M), and RNase-free water, for a total volume of 20 l. The reaction mixture was incubated at 16 °C for 30 min, then at 37 °C for 30 min, and finally at 85 °C for 5 min to inactivate the enzyme. The resulting cDNA could be used for miRNA qPCR. miRNA qPCR used a miRNA fluorescent quantification PCR kit (dye method), which was similar to the target gene qPCR method. Specific miRNA primers were used, and typically three miRNAs were selected for detection, with 5s rRNA serving as an internal reference gene. The reaction system included miRNA primers, cDNA, 2 miRNA qPCR master mix, ROX Reference dye (L) or (H), and water. The miRNA qPCR reaction was usually pre-treated at 95 °C, followed by 40 cycles, each consisting of annealing at 95 °C for 30 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s. After qPCR detection, the relative expression level of miRNA was calculated using the 2−ΔΔCᴛ method by determining the Ct value and comparing it to the expression level of the internal reference gene, thus inferring the expression level of miRNA. All primers used in reverse-transcription quantitative PCR (qRT-PCR) are listed in Supplemental Table S1. We selected two miRNAs, namely mdm-miR156a, and mdm-miR396a, and their four target genes (At3g21360, CYP707A2, EIX2, and SPL9).

    • The study used 15 phloem samples from five different combinations to create a library. The clean data obtained after quality control processing amounted to a total of 366.87 M, with a minimum of 16.64 M clean reads in each sample. The clean reads from each sample were aligned to the Malus x domestica GDDH13 v1.1 reference genome using Bowtie software, with alignment efficiency ranging from 44.28% to 61.14% (Supplemental Table S2). By comparing mapped reads with mature miRNA in the miRbase (v22) database, the study identified 291 known miRNAs (Supplemental Table S3). The majority of these known miRNAs (72.16%) were 21 nt in length (Fig. 1a). The reads that were not mapped to known apple miRNAs in miRbase (v22) were used to predict novel miRNAs using the miRDeep2 package. As a result, 159 novel miRNAs were identified in the study's materials. Their precursors were folded into the secondary hairpin structure by RNAfold, such as Novel miR6, Novel miR16, and Novel miR158 (refer to Fig. 1b for more information). The majority of these novel miRNAs (88.05%) were 24 nt in length. The study also conducted miRNA family analysis using the search software and identified 281 miRNAs from 80 miRNA families (Supplemental Table S4). Among these families, the mdm-MIR171-1, mdm-MIR159, mdm-MIR172, and mdm-MIR395 families accounted for approximately 8.5%, 7.1%, 5.7%, and 5%, respectively (Fig. 1c).

      Figure 1. 

      The length distribution and family of miRNAs identified in apple. (a) Length distribution of known and novel miRNAs. (b) Secondary structures of some novel miRNAs (Novel miR6, Novel miR16, and Novel miR158). (c) Analysis of the miRNA family.

      Dicer and DCL enzymes were known to have strong sequence cleavage preference for 5'U. Analysis on miRNA base bias was used to compare that of identified miRNA with typical miRNA. First base preference of miRNA and base preference on all sites were shown in Fig. 2. First base preference of miRNA has a strong preference for U and tenth base preference for A, which complied with the base preference rule of miRNA, indicating that the miRNA obtained through sequencing was high quality and confidence.

      Figure 2. 

      Analysis of the miRNA nucleotide bias. (a), (b) First base preference of known and novel miRNA in different length, respectively. (c), (d) Base preference on known and novel miRNA at each position, respectively.

    • Expression of miRNAs in each sample was calculated and normalized by TPM algorithm. All miRNA expression in each sample is listed in Supplemental Table S3. Distribution of miRNA expression described overall miRNA expression pattern in each sample. Distribution of TPM density in each sample is shown in Fig. 3a, b. In each sample, the miRNA’s log10TPM was mainly between 1 and 3. The miRNA counts from each sample with ranging from 353 to 411 (Fig. 3c). The expression levels of the miRNAs in the sample’s three replicates tended to be roughly consistent, indicating that biological reproducibility was reliable. We analyzed co-expressed and specifically expressed miRNAs in different samples with Venn plots, and a total 358 miRNAs were expressed together in all combinations (Fig. 3d).

      Figure 3. 

      Analysis of the miRNA expression status in each sample. (a) Boxplot of the overall distribution of miRNA expression in each sample. (b) TPM density distribution (the curves with different colors in the figure represent different samples, the abscissa of the point on the curve represents the logarithm value of the corresponding sample TPM, and the ordinate of the point represents the probability density). (c) Number of miRNA statistics in each sample. (d) Venn analysis of miRNAs in each combination (the number of miRNAs in each combination were taken as three repeated union sets).

    • MiRNA target genes were predicted based on sequences of known miRNAs, novel miRNAs and gene sequences of corresponding species. The Target Finder software (v1.6) was employed for target gene prediction in this study. Summary of miRNA target gene prediction was shown in Table 1. The results demonstrated that 4,792 genes were predicted to be targeted by 436 miRNAs.

      Table 1.  Predicted known and novel miRNAs with targets.

      TypesAll miRNAmiRNA with targetTarget gene
      Known miRNA2912914,194
      Novel miRNA159145783
      Total4504364,792

      In 4,792 target gene, with 4,648 of them annotated in various databases (Fig. 4), including Non-redundant protein sequence database (nr) (4,641), Kyoto encyclopedia of genes and genomes (KEGG) (3,561), Gene ontology database (GO) (3,964), and Clusters of orthologous groups (COG) (1,675). A total of 4,125 target genes were annotated in Malus domestica, and they accounted for 88.9%. In GO enrichment analyses, the major GO terms for target genes were membrane, membrane part, cell, cell part, and organelle for cellular components (CCs); the major GO terms in the binding and catalytic activity were molecular functions (MFs); and the major GO terms in biological processes (BPs) were cellular process and metabolic process. The top three KEGG pathways were plant-pathogen interaction (177), plant hormone signal transduction (143), and starch and sucrose metabolism (119). The top three COG function classes were general function prediction only (256), signal transduction mechanisms (216), and carbohydrate transport and metabolism (164). In general, these target genes were mostly related to plant signal transduction and material transport.

      Figure 4. 

      The analysis of the miRNA’s target genes. (a) NR homologous species distribution. (b) KEGG pathway enrichment analysis. (Y-axis: name of KEGG pathway annotation; X-axis: number and percentage of genes involved in this pathway over counts of all annotated genes). (c) GO annotation classification statistics. (d) COG function classification of consensus sequence (X-axis: the COG functional classifications; Y-axis: the functional categories corresponding to each classification).

    • We compared the miRNAs of B and M to determine whether there was a difference under grafting interstock. During the detection of DEMs, In this project, threshold for defining DEMs was set as |log2(FC)| ≥ 0.58; p-value ≤ 0.05. Fold change (FC) refers to ratio in expression between two samples (groups). P-value represents significancy of difference in expression. We detected DEMs by comparing five pairs of samples (Fig. 5a). A total of 291 DEMs were identified (Supplemental Table S5). The comparison of BS vs MS (B scion vs M scion) revealed 36 upregulated miRNAs and 43 downregulated miRNAs. BR vs MR (B rootstock vs M rootstock), in contrast, had 57 DEMs, of which 22 miRNAs were upregulated and 35 miRNAs were downregulated. We analyzed the DEMs of the BS vs MS and BR vs MR combinations using Venn diagrams (Fig. 5b). A total of 10 miRNAs were upregulated in the scions but downregulated in the rootstocks. A total of six miRNAs were downregulated in the scions but upregulated in the rootstocks. Additionally, after grafting the interstock, there were one DEM that were commonly upregulated in the BS vs MS and BR vs MR groups and three DEMs that were commonly downregulated in these groups (Fig. 5b). The DEMs were subjected to cluster analysis, and a cluster heat map was created (Fig. 5c, d). After grafting the intermediate rootstock, there will be differential expression of tree miRNAs, with a decrease in the expression levels of mdm-miR156 in the scion and mdm-miR172 in the rootstock. The expression levels of the miRNAs in the group’s three replicates tended to be roughly consistent, indicating that biological reproducibility was reliable.

      Figure 5. 

      The analysis of the DEMs in each sample. (a) Number of DEMs in each group. (b) Venn diagram analysis of DEMs (different colors represent different comparison combinations). (c) Heat map of DEMs with BR vs MR. (d) Heat map of DEMs with BS vs MS.

    • We conducted a co-expression trend analysis of the DEMs from the BR vs MR and BS vs MS combinations to study the expression of these DEMs at different samples, with a total of nine expression patterns (Fig. 6 & Supplemental Table S6). Most of the gene trends a and b indicate that the miRNA expression of scion and rootstock will be affected by the dwarfing interstock. mdm-miR172 was Cluster 9. mdm-miR156 was Cluster 6, as well as mdm-miR167 and mdm-miR395 was Cluster 6.

      Figure 6. 

      A co-expression trend analysis of the miRNAs in each sample. (a)−(j) Results plots for each class of the co-expression trend clustering in all DEMs with Cluster 1−9, respectively.

      In this research, we explored the effects of interstock in grafting dwarfing on mRNA and miRNA in the cambium. Generally, individual miRNA not only has no function but is also easily degraded. Only when it forms RISC with Ago protein can it play an effective role, where miRNA plays a guiding role, guiding RISC to specific mRNA regions and exerting its function[37]. Typically, plant miRNA can form complete complementary pairing with the coding region of mRNA and play a role in inhibiting gene expression by inducing mRNA degradation[38]. We used the DEGs from our previously published article and the DEMs from this experiment for joint analysis[35]. Using DEMs as the screening standard, we searched for the mRNA regulated by these DEMs, and focused on analyzing the negative regulatory relationship between miRNA and mRNA (Fig. 7). Meanwhile, KEGG functional enrichment analysis was performed on the resulting mRNA. The relative expression levels of miRNA and target genes are shown in Supplemental Figs S1, S2.

      Figure 7. 

      Sankey diagram of the top 20 miRNA and their target genes. (a) BR vs MR’s top 20 miRNA and their target genes. (b) BS vs MS’s top 20 miRNA and their target genes.

      As shown in Fig. 8, after the joint analysis of differential miRNA and transcriptome sequencing in the cambium of two grafted combinations, 401 target genes were annotated in KEGG, involving a total of 100 pathways. Among them, the most significant enrichment was the starch and sucrose metabolism pathway (ko00500), which involved 42 related genes. In the scion, a total of 371 target genes were annotated in KEGG, involving 97 pathways, mainly enriched in plant hormone signal transduction (ko04075, 30 target genes). The KEGG enrichment results indicate that the interstock may affect the expression of related genes in the root starch and sucrose metabolism pathway and the scion plant hormone signal transduction pathway, thereby affecting the growth and development of the plant, resulting in a dwarf phenotype.

      Figure 8. 

      The KEGG enrichment analysis of DEM’s target genes. (a) BR vs. MR’s DEMs enrichment bubble plots for top 20 KEGG pathways (b) BS vs. MS’s DEMs enrichment bubble plots for top 20 KEGG pathways.

    • To verify miRNA-seq results, we selected three miRNAs, namely mdm-miR156a and mdm-miR396a, and their four target genes (EIX2, CYP707A2, At3g21360, and SPL9) at B and M stages were selected for qRT-PCR. The targeting relationships between these two miRNAs and four genes are shown in Supplemental Table S7. Expression patterns of the three miRNAs and four target genes in B and M were consistent with results of the RNA-seq analysis (Fig. 9). Thus, the transcriptome data were correct and reliable and truly reflected the transcription levels of genes involved in phloem development after grafting dwarfing interstock.

      Figure 9. 

      The linear relationship between RNA-seq and qRT-PCR data. Represented by the linear relationship between RNA-seq and qRT-PCR Data. (a) miR156a, (b) miR396a, (c)–(f) were EIX2, CYP707A2, At3g21360, and SPL9, respectively.

    • MiRNA is a type of small single-stranded endogenous non-coding RNA, usually about 20 to 24 nucleotides long, that regulates gene expression levels by inducing inactivation or degradation of target mRNAs[33]. Through high-throughput sequencing technology, researchers have successfully identified many miRNAs in grafted apples. miRNAs affect the growth and development of grafted apple scions and rootstocks by regulating the expression of target genes[39]. Wang et al. identified 96 miRNAs in the apple rootstock 'M26' through small RNA library construction and Illumina high-throughput sequencing, of which the expression of 23 miRNAs was influenced by grafting[6]. Wang et al. compared the miRNA libraries of apple rootstocks 'M26' and 'JM2', and identified a total of 277 miRNAs, of which 23 miRNAs showed significantly altered expression during grafting[39]. In addition, researchers have also identified some miRNA target genes in grafted apples through target prediction and validation[3941]. For example, the miR156 a/b and the Squamosa promoter binding protein like (SPL) gene family exhibited opposite expression patterns in the rootstock 'M26'. Meanwhile, the study also found that miR395 was expressed at a low level in the rootstock 'M26', possibly because its target gene ATP Sulfurylase 3 (ATPS3) was inhibited during grafting. Li & Mao further identified miRNA target genes in the apple rootstock Malus hupehensis, including transcription factors, hormone synthesis, and signal transduction-related genes, among others[42].

      Grafting can shorten the juvenile phase of fruit trees and improve their stress resistance. It was found that miR156, miR172, and miR395 can transport long distances in grafted plants, and the expression of many miRNAs changed significantly after grafting[37, 4347]. miR156 is an important regulatory factor that participates in the transition from vegetative growth to reproductive growth in plants. It is the main regulatory factor for the transition from the juvenile phase to the adult phase, and has been shown to be the most conserved miRNA in Arabidopsis, maize, rice, and other plants[43]. Overexpression of miR156 under long-day conditions can cause delayed flowering and a decrease in the expression of SPL family members[44, 48]. In this experiment, the expression level of mdm-miR156 decreased in the interstock combination scion, which may promote the transition of the plant to the adult phase and obtain early flowering characteristics. Low expression of miR172 may shorten the juvenile phase[16]. In this study, mdm-miR172 was downregulated in the interstock of the grafting. miR395 is a key factor that regulates plant sulfur metabolism. It targets three ATP sulfurylase genes, ASP1, ASP3, and ASP4. When miR395 is induced by low sulfur stress, the transcription level of ASP1 decreases. When sulfur supply is normal, the transcription level of ASP1 increases, and the expression of mdm-miR395 is completely suppressed[49].

      In this study, it was found that after grafting the interstock, the expression level of mdm-miR395 was downregulated, and its effect on plant growth and development was not clear. According to the regulatory mechanism between miRNA and target genes, miRNA and their corresponding target genes show a negative correlation in expression regulation. This study found that the expression of most miRNAs and their target genes showed a negative correlation after grafting, while a small number of miRNAs and target genes showed a positive correlation or no correlation. Generally, gene expression is influenced by more than one factor, and the expression of miRNA target genes is also regulated by factors other than miRNA. To date, the complex regulation of miRNA on target genes has not been fully elucidated[50]. This result preliminarily reveals the expression patterns of miRNA and their target genes between grafted apple rootstocks. Through KEGG enrichment analysis, the effects of these differentially expressed genes on plant physiology and growth and development mainly manifested in the pathways of root starch and sucrose metabolism and the plant hormone signaling pathways in the scion. The problems of how to regulate gene expression and transcription need further research. In this study, to clarify the genetic interaction mechanism between the grafted rootstocks and to explore miRNAs that play important roles in plant growth, stress resistance, and anti-aging, are of great practical significance for improving grafting resistance, disease resistance, yield, and quality through genetic improvement of rootstocks and the interaction between rootstocks and scions.

      Compared with the control, there were 79 differentially expressed miRNAs in the scion and 57 differentially expressed miRNAs in the interstock, including 36 upregulated miRNAs and 43 downregulated miRNAs in the scion and 22 upregulated miRNAs and 35 downregulated miRNAs in the interstock. Through KEGG enrichment analysis, the interstock may affect the expression of related genes in the root starch and sucrose metabolism pathways and the scion plant hormone signaling pathways, thereby affecting plant growth and development and causing plant dwarfing. The downregulation of mdm-miR156 in the scion and mdm-miR172 in the interstock may promote the plant's transition to the adult phase and obtain early flowering characteristics. This study not only provides responsive miRNAs and target genes to apple grafting, but also enriches and supplements the mechanism of apple grafting. It also provides candidate miRNAs and target genes that respond to apple grafting-induced dwarfing and early flowering, enriches the research on the mechanism of apple grafting-induced dwarfing and early flowering, and provides a new theoretical basis for the creation of new apple interstock varieties.

    • Grafting is a common technique in horticulture for improving fruit tree performance. A total of 281 miRNAs were categorized into 80 families, including 159 newly discovered ones. Grafting with an interstock led to 79 differentially expressed miRNAs (DEMs) in the scion, with 36 up-regulated and 43 down-regulated miRNAs, and 57 DEMs in the rootstock, with 22 up-regulated and 35 down-regulated miRNAs, when compared to the control. Grafting using a dwarfing interstock influenced DEMs in the entire plant, causing a decrease in mdm-miR156 expression in the scion and mdm-miR172 in the rootstock. Analysis of DEMs and their target genes indicated that miRNAs influence scion growth by affecting aspects such as plant hormone synthesis, signal transduction, and nutrient balance. Combined miRNA and mRNA analysis revealed that the dwarfing effect of the interstock might impact genes related to starch and sucrose metabolism in the rootstock and plant hormone signaling in the scion. KEGG enrichment analysis identified pathways affecting plant physiology and growth, including root starch and sucrose metabolism, and hormone signaling. Notably, differentially expressed miRNAs in the scion and interstock influence related genes and contribute to plant dwarfing. The downregulation of mdm-miR156 and mdm-miR172 could promote the transition to the adult phase and early flowering. This research sheds light on miRNA regulation in grafted apples, enhancing our grasp of grafting-induced changes and providing a foundation for creating new apple interstock varieties.

    • The authors confirm contribution to the paper as follows: study conception and design: Gao Y, Feng J, Wang D; data collection: Gao Y, Wang K, Li L; analysis and interpretation of results: Li Q, Liu Z, Li Z, Wang G, Wang L, Lu X; draft manuscript preparation: Li Q, Feng J, Sun S, Lu X, Guo H. All authors reviewed the results and approved the final version of the manuscript.

    • The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) in National Genomics Data Center (Nucleic Acids Res 2022), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA008505) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa[51, 52].

      • We would like to thank OmicShare tools, a free online platform for plotting (https://www.omicshare.com/tools, accessed on 13 February 2023), the ImageGP, an easy-to-use data visualization web server for scientific researchers (http://www.ehbio.com/ImageGP/, accessed on 26 February 2023), and the Hiplot team (https://hiplot.com.cn, accessed on 2 February 2023) for providing technical assistance and valuable tools for data analysis and visualization. This work was supported by the Agricultural Science and Technology Innovation Program (CAAS-ASTIP-2021-RIP-02) and Fundamental Research Funds for Central Non-profit Scientific Institution, grant number 161018201614.

      • The authors declare that they have no conflict of interest.

      • Copyright: © 2023 by the author(s). Published by Maximum Academic Press, Fayetteville, GA. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Figure (9)  Table (1) References (52)
  • About this article
    Cite this article
    Li Q, Gao Y, Wang K, Sun S, Lu X, et al. 2023. Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock. Fruit Research 3:34 doi: 10.48130/FruRes-2023-0034
    Li Q, Gao Y, Wang K, Sun S, Lu X, et al. 2023. Identification of microRNAs and target genes in apple (Malus domestica) scion and rootstock with grafted interstock. Fruit Research 3:34 doi: 10.48130/FruRes-2023-0034

Catalog

  • About this article

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return