Genome-wide identification and analysis of miR396 family members and their target GRF genes in rubber tree (Hevea brasiliensis)

Abstract

Background

MicroRNAs (miRNAs) are key regulators of gene expression as they play crucial roles at the post-transcriptional level. In particular, the miR319-GRF module is an important gene regulatory network in plants, extensively involved in processes such as plant growth and development. Although miR396 is one of the most conserved miRNA families, its role in rubber trees remains poorly understood. In this study, bioinformatics analysis, including target prediction, was performed to reveal the evolutionary and expression patterns of the Hbr-miR396 family members.

Results

A total of six Hbr-miR396 members were identified, distributed across four chromosomes. Secondary structure analysis revealed that the precursor sequences of the six Hbr-miR396 members could form a typical stem-loop (hairpin) structure. Sequence analysis show that the members of the Hbr-miR396 family form three mature sequences. Furthermore, phylogenetic analysis demonstrated that the Hbr-miR396 family members are closely related to those from cassava. Eight members of the growth regulatory factor (GRF) family were predicted as potential targets of Hbr-miR396. The dual-luciferase assays also confirmed that Hbr-miR396b strongly inhibited the expression of HbrGRF3. Expression analysis of the HbrGRF targets in different tissues revealed that HbrGRFs are mainly expressed in the cambium and flowers. Therefore, Hbr-miR396 may potentially regulate growth and floral organ development in rubber trees by targeting HbrGRFs.

Conclusions

The data presented in this study offer valuable insights into the functional and molecular regulatory mechanisms of the miR396-GRFs module in rubber tree growth and development, laying a foundation for further investigation into its biological roles in enhancing both rubber production and timber quality.

Background

MicroRNAs (miRNAs) are endogenous RNA molecules of approximately 22 nucleotides. They play important regulatory roles in plants and animals by targeting mRNAs for cleavage or translational repression [1]. They were first discovered in the nematode Caenorhabditis elegans in 1993 [2]. However, the first miRNA in plants was reported in 2002, when miRNAs were discovered in Arabidopsis thaliana [3]. Since then, with the development and application of high-throughput sequencing [4], an increasing number of miRNAs have been identified in various plant species, such as in radish [5], tea plant [6], strawberry [7], and so on. To date, miRNAs from 271 species have been studied [8].

The miRNAs are regarded as key regulators of gene expression because they play a crucial role at the post-transcriptional level by blocking translation or promoting the degradation of target mRNAs [1]. miRNAs are involved in various biological processes, including plant growth and development [9,10,11,12], and responses to biotic and abiotic stresses [6, 13, 14]. Among them, the miR396 family is one of the well-conserved families in plants, which was first identified using a devised computational procedure, and the expression of two members, miR396a and miR396b, was validated in Arabidopsis thaliana through Northern blot analysis in 2004 [14]. In the A. thaliana genome, there are two miR396-encoding loci: athMIR396a and athMIR396b [14]. Members of the miR396 family extensively target the growth-regulating factor (GRF) genes to mediate plant growth and development [15,16,17,18]. In Arabidopsis, ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth [19]. Further, miR396 affects the growth and development of lettuce leaves by inhibiting the expression of plant-specific GRFs [20]. In A. thaliana, compared to wild-type plants, the overexpression of miR396 weakens the regulatory effect of GRFs, leading to abnormal root phenotype development and severely impairing root growth [9]. Therefore, the miR396-GRF module plays an important regulatory role in root organ development in plants. The miR396-GRF module is also involved in the development of rice floral organs, when the miR396d overexpression had resulted in the production of more petals as compared to the wild-type plants [12]. Additionally, studies have shown that the heterologous expression of miR396c in Populus results in altered flower development and a significant decrease in the expression of the NtGRF gene in tobacco [10]. In addition, studies have shown that rice Osa-miR396c enhances salt tolerance in the transgenic plants by improving water retention, chlorophyll content, membrane integrity, and Na⁺ exclusion [21]. Recent research has also revealed that the miR396b-GRF6 module enhances rice salt tolerance by reducing ROS accumulation and increasing ROS-scavenging enzyme activity [22]. Recent studies have shown that miR396 also plays a critical role in the resistance of alfalfa against insects, as down-regulation of miR396 expression in MIM396 transgenic plants enhances their resistance to Spodoptera litura, primarily by increasing the lignin content by regulating the related biosynthetic pathways [23]. These studies, thus have revealed the functional mechanisms of the miR396-GRF module.

The rubber tree (Hevea brasiliensis Muell. Arg.), a perennial species of the Euphorbiaceae family, is the primary source of natural rubber [24]. Rubber trees are used to produce a variety of products, such as pulpwood for papermaking, rubber-based composite boards, furniture, and fine woodworking products [25,26,27]. The economic importance of the rubber tree and its growing demand have promoted its extensive domestication. Given the increasing economic and application value of the rubber tree, research on its miRNAs has been reported a few times, most of which is related to high-throughput sequencing [28, 29]. However, studies on the systematic analysis of the miR396 family members in the rubber tree and the prediction of their target genes have remained limited.

This study focused on the rubber tree Hbr-miR396 family, by conducting a systematic analysis of its precursor and mature sequences. Furthermore, target gene prediction was performed by using bioinformatics tools. The interaction between Hbr-miR396b and its target genes was further validated through a dual-luciferase reporter assay. The results of this study provide a theoretical support for further understanding the potential role of Hbr-miR396 family members in the growth and development of rubber trees.

Results

Sequence analysis and chromosomal localization of the Hbr-miR396 family members

The Hbr-miR396 family in the rubber tree genome, consists of six members: Hbr-miR396a, b, c, d, e, and f (Table 1, Fig. 1A). The precursor sequences of the six Hbr-miR396 members range in length from 65 to 120 nucleotides and contain three highly conserved regions. The members of the miR396 family share similar 5’ and 3’ sequences, while Hbr-miR396f shows a reversed situation at both the 5’ and 3’ ends (Fig. 1A).

[IMAGE OMITTED: SEE PDF]

The mature Hbr-miR396 sequences are derived from the highly conserved nucleotides at positions 1–21 of the precursor sequence. Hbr-miR396a, Hbr-miR396c, and Hbr-miR396e share the same mature sequence, while Hbr-miR396b and Hbr-miR396d have identical mature sequences. In contrast, the mature sequence of Hbr-miR396f shows nucleotide variation compared to other members of the Hbr-miR396 family (Fig. 1B; Table 1). The star sequences of all six Hbr-miR396 members consist of 21 nucleotides, which exhibit a high degree of conservation (Fig. 1C; Table 1).

Chromosomal localization results show that Hbr-MIR396a and Hbr-MIR396d are located on chromosome CM021231.1, while Hbr-MIR396e and Hbr-MIR396b are found on chromosome CM021239.1. In contrast, Hbr-MIR396f and Hbr-MIR396c are located on chromosomes CM021229.1 and CM021236.1, respectively (Supplementary Fig. 1, Table 1).

[IMAGE OMITTED: SEE PDF]

Analysis of the precursor sequence characteristics of the Hbr-miR396 family members

The precursor sequences of the six Hbr-miR396 family members vary in length. The sequence of Hbr-MIR396b is the longest (120 bp), while that of Hbr-MIR396f is the shortest (65 bp). Secondary structure predictions showed that all the precursor sequences could form the stem-loop structures, which are commonly called hairpin structures. All the mature Hbr-miR396 family members are located on the 5’ arm of their respective precursors (Fig. 2; Table 1).

[IMAGE OMITTED: SEE PDF]

Phylogenetic analysis of the Hbr-miR396 family members

Phylogenetic analysis of the precursor sequences revealed that the miR396 family members of P. tomentosa, P. euphratica, R. communis, M. esculenta, and A. thaliana could be divided into three groups (Groups I–III). Hbr-MIR396d and Hbr-MIR396b cluster within Group I, while Hbr-MIR396f, Hbr-MIR396c and Hbr-MIR396e cluster within Group II. Compared to P. tomentosa, P. euphratica, and A. thaliana, the rubber tree’s Hbr-MIR396 family has a closer phylogenetic relationship with that of M. esculenta. This is primarily reflected in the clustering patterns, where the rubber tree MIRR396 family members, Hbr-MIR396d and Hbr-MIR396b, group together with M. esculenta miRNAs Mes-MIR396b and Mes-MIR396e. Hbr-MIR396f clusters with Mes-MIR396a and Hbr-MIR396e groups with Mes-MIR396f (Fig. 3).

[IMAGE OMITTED: SEE PDF]

Transcript accumulation analysis of the Hbr-miR396 family members

Analysis of the levels of transcript accumulation of the Hbr-MIR396 family members in the mature and young leaves of the rubber tree showed that Hbr-MIR396b and Hbr-MIR396d are highly expressed in mature leaves; their levels being the highest compared to the other family members. In contrast, Hbr-MIR396f, Hbr-MIR396c, Hbr-MIR396a, and Hbr-MIR396e were highly expressed in young leaves (Supplementary Fig. 2).

Prediction and physicochemical characterization of the GRF genes targeted by Hbr-miR396

Previous studies have shown that miR396 primarily influences the growth and development of plants by mediating the GRF regulatory module. Therefore, this study focused on systematically analyzing the predicted GRF genes that could be targeted by Hbr-miR396. We found that Hbr-miR396 targets 8 members of the HbrGRF family, and that the target sites are all located in the exons. The protein sequences of these 8 HbrGRFs range from 232 to 596 amino acids in length, with theoretical isoelectric points (pI) ranging from 7.75 to 9.43 and molecular weights ranging from 25.51 to 67.01 kDa. Prediction of their subcellular localization indicates that all the 8 HbrGRFs may be localized in the nucleus. HbrGRF6 was also predicted to localize in the chloroplast (Table 2).

[IMAGE OMITTED: SEE PDF]

Gene structure and expression analysis of HbrGRFs

Using the MEME online software prediction platform, the conserved protein sequences of the GRF genes targeted by Hbr-miR396 were analyzed. A total of 10 motifs were identified, and distributed as following: HbrGRF1 and HbrGRF7 all contained Motif 7, Motif 2, Motif 1, Motif 8, Motif 4, Motif 10 and Motif 6; HbrGRF5 and HbrGRF6 both contained Motif 2, Motif 9, Motif 1, Motif 3, Motif 5, Motif 4 and Motif 10; whereas HbrGRF2 and HbrGRF4 only contained Motif 2 and Motif 1.

Using the CDD (Conserved Domain Database) online tool, the conserved domains in the GRF target genes were predicted. A total of four conserved domains were identified. All HbrGRF proteins contained the complete conserved domains at the N-terminus, including the WRC (Motif 1) and QLQ (Motif 2) domains. Additionally, the QLQ domain was located before the WRC domain in these proteins. Moreover, HbrGRF7 also contained the GPHR_N and ABA_GPCR domains (Fig. 4A).

The GRF families in P. trichocarpa, P. euphratica, R. communis, M. esculenta, and A. thaliana consist of 9, 20, 9, 15, and 9 members, respectively. Using the protein sequences of these GRF genes and the 8 GRF genes from H. brasiliensis, a phylogenetic tree was constructed to analyze the evolutionary relationships among the GRF gene families. These GRF genes were clustered into 5 groups (Group I-V). The GRFs from H. brasiliensis, P. trichocarpa, P. euphratica, R. communis, M. esculenta, and A. thaliana were distributed across all the groups, indicating that the GRF genes in these plants may share a similar evolutionary trajectory (Fig. 4B).

miR396 influences plant growth and development by regulating the expression of its target GRF genes [30]. To investigate the expression patterns of HbrGRFs, transcriptome data from seven different rubber tree tissue types (cambia region, female flower, inner bark, leaves, male flower, tapped latex, and virgin latex) were analyzed. The results revealed that the expression patterns of HbrGRFs exhibited significant tissue specificity. These genes were mainly expressed in tissues with highly active cells, such as the cambia region, female flowers, and male flowers (Fig. 4C).

[IMAGE OMITTED: SEE PDF]

Experimental validation of the interaction between Hbr-miR396b and HbrGRF3

Through bioinformatics analysis, the full CDS sequence of HbrGRF3 is comprised of 699 bp, and Hbr-miR396b targets the HbrGRF3 sequence from 534 to 555 bp, where a total of 19 complementary bases is observed (Fig. 5A).

We selected the 534–555 bp sequence of HbrGRF3 and cloned it in the pGreen II 0800-miRNA LUC vector. This created a recombinant luciferase fragment containing the Hbr-miR396b cleavage sequence, which was named HbrGRF3-Luc as the ‘reporter’. The Hbr-miR396b precursor sequence was cloned into the pBI121 vector as the ‘effector’ (Fig. 5B). The pGreen II 0800-LUC empty vector was used as a positive control. Both constructs were transformed into A. tumefaciens (GV1301), which was then infiltrated into tobacco leaves to establish a dual-luciferase reporter system. The Luc/Ren fluorescence ratio was measured to determine the mode of action of Hbr-miR396b on HbrGRF3. Compared to the control group (the co-expression of pGreen II 0800-LUC empty vector and Hbr-miR396b), the experimental group (HbrGRF3-Luc and Hbr-miR396b) exhibited a significant reduction in the Luc/Ren fluorescence expression (Fig. 5C).

In addition, the qRT-PCR analysis of the tissue-specific expression of Hbr-miR396b and HbrGRF3 in tissue-cultured rubber seedlings showed that Hbr-miR396b had the lowest transcript level in the stem, while its target gene HbrGRF3 exhibited the hightest expression level in the strem tissue. These results indicates that Hbr-miR396b and HbrGRF3 exhibit different tissue-specific expression patterns (supplementary Fig. 3).

[IMAGE OMITTED: SEE PDF]

Discussion

In recent years, miRNAs have garnered significant attention as newly discovered gene regulatory factors and have rapidly become a research hotspot in the field of life sciences. They target mRNAs through complementary base-pairing interactions, leading to transcript degradation or inhibition of corresponding genes expression, and thus, they play key roles in regulating various biological processes. With the continuous improvement of plant genome databases in recent years, research in this area has also increased. Among the miRNAs, miR396 is a widely conserved family, and is a crucial component of plants [30, 31]. In this study, bioinformatics methods were employed to analyze the characteristics and functions of Hbr-miR396 family members.

A large number of studies have shown that miR396 regulates the expression levels of its target genes, the GRFs. Thereby, it modulates plant cell division and differentiation. It plays an important role in various biological processes, such as the development of plant floral organs [10, 12], root development [9, 32], leaf and fruit development [33], and responses to both biotic and abiotic stresses [22, 34]. Currently, research on miR396 has mainly focused on model plants like rice, Arabidopsis, and tobacco, while studies in rubber trees are relatively limited. In this study, the results of the evolutionary analysis of the Hbr-miR396 family indicate that miR396 is conserved across species. This study found that the Hbr-miR396 family in rubber trees consists of 3 mature and 6 precursor sequences. Different precursor sequences could produce identical mature sequences. Except for Hbr-miR396b, each mature Hbr-miR396 originates from the highly conserved 1–21 positions at the 5’ end of the corresponding precursor sequence. Apart from Hbr-miR396f, the mature Hbr-miR396 sequences are highly similar to each other, with only 1 or no base differences among them.

The miR396 family exhibits a high degree of evolutionary conservation across species [17]. The results of this study indicate that the precursor sequences of miR396 in rubber trees are more closely related to cassava. This indicates that the evolutionary relationship of the miR396 gene family is consistent with the phylogenetic relationships of the species themselves. Tissue specific expression analysis revealed that Hbr-MIR396b and Hbr-MIR396d are highly expressed in the mature leaves, while Hbr-MIR396f, Hbr-MIR396c, Hbr-MIR396a, and Hbr-MIR396e are highly expressed in young leaves. This indicates that different members of the miR396 family may play distinct roles at different stages of leaf development. This differential expression is likely closely associated with the physiological functions and regulatory mechanisms of the leaves at different developmental stages. This could be mediated by miR396 through modulation of GRFs.

GRFs are plant-specific transcription factors that play a critical role at various stages of plant growth and development. The N-terminus of the GRF proteins typically contains two conserved functional domains, QLQ and WRC [30]. In this study, eight of the Hbr-miR396-targeted GRFs exhibited complete conservation of the domains, WRC (Motif 1) and QLQ (Motif 2), where the QLQ domain was positioned before the WRC domain. It has been reported that miR396 cleaves the conserved WRC domain sequence of GRF genes, which negatively regulates the GRF gene expression [13]. The motif prediction for HbrGRFs show that the 8 HbrGRF proteins possess varying numbers of motifs, which indicates that they may perform different functions in regulating the growth and development of rubber trees.

Hbr-miR396 may primarily influence the growth and development processes by targeting HbrGRFs and by regulating their expression levels. Therefore, analyzing the expression patterns of HbrGRF genes in rubber trees helped us in understanding the function of the Hbr-miR396-HbrGRF module. Expression analysis of HbrGRFs in different tissues revealed a clear tissue-specificity as they were predominantly expressed in the cambia region and the flowers. These findings indicate that HbrGRFs in rubber trees may be mainly involved in regulating the development of floral organs and the growth of the stem. These results are consistent with the previous studies, where members of the GRF gene family are expressed at higher levels in organs or tissues with high developmental activity, and at lower levels in mature tissues or organs [35]. Additionally, in Brassica napus, BnGRF regulate the development of floral organs [36].

As post-transcriptional regulators, miRNAs modulate signal transduction, stress responses, and hormone biosynthesis by either degrading target genes or inhibiting their translation [37]. The qRT-PCR analysis of tissue-specific expression of miR396b and HbrGRF3 in the rubber seedlings showed an inverse correlation between root and stem, which indicated that Hbr-miR396b may negatively regulate the expression of HbrGRF3. Although this association is not observed in leaves, leading to differences in the regulatory relationship between miR396 and GRF, we speculate that this may be a “precise regulatory strategy” formed by plants during long-term evolution. Through “miRNA expression profile differentiation, target gene functional division, and developmental stage adaptation”, the same regulatory module (miR396-GRF) achieves “functional adaptation” in different organs. Moreover, through bioinformatics analysis, we have found that miR396b targets HbrGRF3 at the 534–555 bp position. To further validate this prediction, we conducted a dual-luciferase reporter assay, which showed that miR396b strongly inhibited the luciferase expression of HbrGRF3_Luc. In addition, we performed a multiple sequence alignment on the loci of the 8 HbrGRF genes targeted by members of the miR396 family. The results showed that the target sequences are highly conserved and have a high degree of similarity to the HbrGRF3 gene (Supplementary Fig. 5) which has also been confirmed in soybeans [13]. This may also indicate that other Hbr-miR396 family members regulate the expression of the target HbrGRF genes in the same manner. This further supports the hypothesis that Hbr-miR396 likely exerts its regulatory effect on the target HbrGRFs by cleaving its mRNA sequence, thereby inhibiting its expression, and ultimately regulating plant growth and development. Furthermore, the verification of each Hbr-miR396 member will be a primary focus of our subsequent research, utilizing transient expression systems or VIGS (Virus-Induced Gene Silencing) techniques in the rubber tree to further verify the function of the miR396-GRF module in the native species.

Despite certain limitations in the verification of stable inheritance in the original species, such as gain- or loss-of-function studies in rubber tree, to unequivocally establish the causal relationship between miR396 and HbrGRF. Indeed, genetic transformation and the generation of stable transgenic lines in Hevea brasiliensis remain exceptionally challenging due to its long life cycle, low transformation efficiency, and the recalcitrance of its tissues in culture [38, 39]. We believe that the convergent evidence from heterologous validation, spatiotemporal expression patterns, and evolutionary conservation presents a compelling case for the role of the miR396b-HbrGRF3 module in rubber tree development. Our work serves as a crucial foundation and provides valuable genetic resources (gene sequences, expression profiles) for the community. This information will provide potential strategies and directions for rubber treee breeding programs. The identified genes (Hbr-miR396 and HbrGRFs) serve as direct targets for future genetic engineering or gene editing approaches aimed at modulating growth and development of rubber trees.

Materials and methods

Genome-wide identification and sequence analysis of Hbr-miR396 and its target genes

The stem-loop and the mature sequences of the Hbr-miR396 family members were downloaded from the PmiREN database (https://www.pmiren.com/) [40]. The base conservation of these sequences was analyzed using WebLogo 3 (https://weblogo.threeplusone.com/) [41]. RNAfold WebServe was employed to predict the secondary stem-loop structures of the Hbr-miR396 precursors (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi).

The whole genome sequence and annotation files of the rubber tree were downloaded from the NCBI (National Center for Biotechnology Information) database (http://www.ncbi.nlm.nih.gov/) under BioProject accession number PRJNA587314. The miRNA-precursor sequences were then subjected to a BLAST alignment against the rubber tree genome to determine the gene locus information of the Hbr-MIR396 members. The Gene Location Visualization tool in TBtools was used to construct the chromosome localization map [42].

The psRNATarget tool (http://plantgrn.noble.org/psRNATarget/), with default parameters, was used to predict the target genes of Hbr-miR396 and to analyze the interaction sites between the Hbr-miR396 and the candidate target genes. Subsequently, the online tool, PlantTFDB (https://planttfdb.gao-lab.org/), was used to screen for target genes that belong to the GRF gene family for further analysis [43]. Additionally, the ExPASy tool (https://web.expasy.org/protparam/) was used to predict the molecular weight (MW) and isoelectric point (pI) of the candidate target genes’ proteins. The subcellular localization of HbrGRF genes was further predicted with the help of the Cell-PLoc 2.0 (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/) online tool [44].

Phylogenetic analysis of Hbr-miR396 and its target HbrGRF genes

The precursor stem-loop sequences of miR396 for Populus trichocarpa, Populus euphratica, Ricinus communis, Manihot esculenta, and A. thaliana were downloaded from the PmiREN database. The GRF protein sequences for P. trichocarpa, P. euphratica, R. communis, M. esculenta, and A. thaliana were downloaded from the PlantTFDB database.

The precursor sequences of miR396 and the GRF protein sequences were aligned with the help of ClustalW [45] with default parameters for multiple sequence alignment. A phylogenetic tree was then constructed with the help of the MEGA v7.0 software, using the neighbor-joining (NJ) method with a bootstrap value of 1000 [46].

Sequence and structural analysis of HbrGRF genes family targeted by Hbr-miR396

The conserved motifs of HbrGRFs were analyzed with the help of the Multiple Expectation Maximization Algorithm for motif prediction (MEME, http://meme-suite.org/). Prediction of the conserved domains of the HbrGRF genes, which could be targeted by Hbr-miR396, was performed by using the CDD tool (https://www.ncbi.nlm.nih.gov/Structure/cdd) [47]. Subsequently, the identified motifs and gene structures were further visualized with the help of the TBtools.

Expression analysis of Hbr-miR396 and target HbrGRF genes

The expression data of the Hbr-miR396 family for the young and the mature leaves were downloaded from the PmiREN database for the purpose of differential expression analysis. Transcriptome data of H. brasiliensis ‘Reyan 8–79’ (PRJCA004986) were obtained from the National Genomics Data Center (NGDC, https://ngdc.cncb.ac.cn/) [48]. The raw data in Fastq format were first processed using fastp for read trimming, data filtering, and quality control. Clean reads were obtained by removing reads containing sequencing adapters, reads with poly-N, and low-quality reads. These clean reads were then aligned to the reference genome sequence using HISAT2. The gene expression levels were evaluated using FPKM values, and this step was performed using Stringtie2.1.2. This data was used to analyze the differential expression of the Hbr-miR396 target genes, HbrGRFs, in the seven tissue types: cambium region, female flower, inner bark, male flower, tapped latex, leaves and virgin latex. Transcript levels of Hbr-miR396 and HbrGRFs were estimated with the FPKM values and scale method was normalized, and heatmaps were generated by using the TBtools.

Dual-luciferase assay

To verify the miRNA targets, the binding-site sequences of HbrGRF3 were inserted between XbaI and EcoRI restriction site of the pGreen II 0800-miRNA LUC vector, downstream of the cauliflower mosaic virus (CaMV) 35 S promoter, and were used as reporters. The precursor sequence of Hbr-MIR396b was amplified from the cDNA samples of rubber tree (Renyan 7–33−97) [49]. These PCR products were cloned into between XbaI and SacI restriction site a binary plant expression pBI121 vector to generate 35 S::miR396b, and used as an effector. The plasmids were transferred into A. tumefaciens GV1301 and co-infiltrated into the N. benthamiana leaves. The luciferase activities were measured with the help of the Dual-Luciferase Reporter Assay System (Promega, Madison, WI, USA), according to the manufacturer’s instructions.

RNA extraction and quantitative reverse transcription PCR (qRT-PCR)

The rubber tissue culture seedlings that are 1 year old were selected for sampling of roots, stems, and leaves. Total RNA was extracted by using the TRIzol reagent (Invitrogen), which was followed by a DNase treatment to remove genomic DNA contamination. Complementary DNA (cDNA) was synthesized by using the PrimeScript™ RT reagent kit and gDNA Eraser (TaKaRa), according to the manufacturer’s instructions. The qRT-PCR was performed with the help of SYBR Premix ExTaq™ (TaKaRa) on a thermal cycler. PCR amplification was performed at a thermal cycler heat block and incubated at 95 ◦C for 30 s, followed by 40 cycles of 95 ◦C for 5 s and 60 ◦C for 30 s. The Ubiquitin genes served as internal controls for normalization. Mature Hbr-miR396b was reverse transcribed by using a specific stem-loop reverse transcription primer, and analyzed by stem-loop qRT-PCR detection Kit (Aidlab, Beijing, China). PCR amplification was performed at a thermal cycler heat block and incubated at 94 ◦C for 3 min, followed by 40 cycles of 94 ◦C for 10 s and 60 ◦C for 32 s. The U6 genes served as internal reference gene for normalization.

For melting curve analysis, samples were denatured at 95 ◦C for 30 s, then cooled to 65 ◦C for 30s. Fluorescence signals were collected continuously from 65 ◦C to 95 ◦C at 0.5 ◦C per 5 s. The expression levels of the Hbr-miR396b and HbrGRF3 gene were calculated using the 2–△△Ct method against the internal controls, with the expression levels of Hbr-miR396b and HbrGRF3 in roots are set to 1. Three technical replicates per sample were analyzed to ensure the reliability. The asterisks indicated the significant differences based on Student’s t-test (*, P < 0.05; **, P < 0.01;***, P < 0.001). The gene-specific primers that were used to determine expression levels are listed in Supplementary Tables 1 and Melt Curve and Amplification plots of Hbr-miR396b and HbrGRF3 are presented in Supplementary Fig. 4.

Conclusions

In this study, we have identified 6 Hbr-miR396 family members from the rubber tree genome, which could target 8 HbrGRF genes. Through molecular validation, we also have confirmed the interaction between Hbr-miR396b and HbrGRF3. As post-transcriptional regulators, the Hbr-miR396 family members could primarily regulate plant growth and development by targeting GRF genes. Most HbrGRF genes exhibit higher expression levels in the actively developing tissues and organs, such as the cambia region and flowers. Therefore, we hypothesize that Hbr-miR396 may potentially regulate growth and floral organ development in rubber trees by targeting HbrGRFs. This study is the first comprehensive identification of the Hbr-miR396-HbrGRFs module in the rubber trees. This research provides a foundation for further elucidating the molecular mechanisms of the Hbr-miR396-HbrGRFs module in the growth and development of rubber trees.

Data availability

All data generated or analyzed during this study are included in this paper and its supplementary information files. The whole genome sequence and annotation files of H. brasiliensis were downloaded from the NCBI (National Center for Biotechnology Information) database (http://www.ncbi.nlm.nih.gov/) under BioProject accession number PRJNA587314. Additionally, the RNA-seq data of *H. brasiliensis* (PRJCA004986) were obtained from the National Genomics Data Center (NGDC, https://ngdc.cncb.ac.cn/).

References

Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281–97. https://doi.org/10.1016/s0092-8674(04)00045-5.

Olsen PH, Ambros V. The lin-4 Regulatory RNA Controls Developmental Timing in Caenorhabditis elegans by Blocking LIN-14 Protein Synthesis after the Initiation of Translation. Developmental Biology. 1999;216:671–80. https://doi.org/10.1006/dbio.1999.9523.

Llave C, Xie Z, Kasschau KD, Carrington JC. Cleavage of Scarecrow-like mRNA targets directed by a class of Arabidopsis MiRNA. Science. 2002;297:2053–6. https://doi.org/10.1126/science.1076311.

Ma X, Tang Z, Qin J, Meng Y. The use of high-throughput sequencing methods for plant MicroRNA research. RNA Biol. 2015;12:709–19. https://doi.org/10.1080/15476286.2015.1053686.

Liu W, Xu L, Wang Y, Shen H, Zhu X, Zhang K, et al. Transcriptome-wide analysis of chromium-stress responsive MicroRNAs to explore miRNA-mediated regulatory networks in radish (Raphanus sativus L). Sci Rep. 2015;5:14024. https://doi.org/10.1038/srep14024.

Jeyaraj A, Wang X, Wang S, Liu S, Zhang R, Wu A, et al. Identification of regulatory networks of MicroRNAs and their targets in response to Colletotrichum gloeosporioides in tea plant (Camellia sinensis L). Front Plant Sci. 2019;10. https://doi.org/10.3389/fpls.2019.01096.

Liang Y, Guan Y, Wang S, Li Y, Zhang Z, Li H. Identification and characterization of known and novel MicroRNAs in strawberry fruits induced by botrytis cinerea. Sci Rep. 2018;8:10921. https://doi.org/10.1038/s41598-018-29289-7.

Kuang W, Qin D, Huang Y, Liu Y, Cao X, Xu M. Analysis of the miR482 gene family in plants. Genes. 2024;15:1043. https://doi.org/10.3390/genes15081043.

Bao M, Bian H, Zha Y, Li F, Sun Y, Bai B, et al. miR396a-Mediated basic Helix–Loop–Helix transcription factor bHLH74 repression acts as a regulator for root growth in Arabidopsis seedlings. Plant Cell Physiol. 2014;55:1343–53. https://doi.org/10.1093/pcp/pcu058.

Baucher M, Moussawi J, Vandeputte OM, Monteyne D, Mol A, Pérez-Morga D, et al. A role for the miR396/GRF network in specification of organ type during flower development, as supported by ectopic expression of opulus trichocarpa miR396c in Transgenic tobacco. Plant Biol. 2013;15:892–8. https://doi.org/10.1111/j.1438-8677.2012.00696.x.

Baulies JL, Bresso EG, Goldy C, Palatnik JF, Schommer C. Potent Inhibition of TCP transcription factors by miR319 ensures proper root growth in Arabidopsis. Plant Mol Biol. 2022;108:93–103. https://doi.org/10.1007/s11103-021-01227-8.

Liu H, Guo S, Xu Y, Li C, Zhang Z, Zhang D, et al. OsmiR396d-Regulated OsGRFs function in floral organogenesis in rice through binding to their targets OsJMJ706 and OsCR4. Plant Physiol. 2014;165:160–74. https://doi.org/10.1104/pp.114.235564.

Liu W, Zhou Y, Li X, Wang X, Dong Y, Wang N, et al. Tissue-Specific regulation of Gma-miR396 family on coordinating development and low water availability responses. Front Plant Sci. 2017;8:1112. https://doi.org/10.3389/fpls.2017.01112.

Jones-Rhoades MW, Bartel DP. Computational identification of plant MicroRNAs and their Targets, including a Stress-Induced MiRNA. Mol Cell. 2004;14:787–99. https://doi.org/10.1016/j.molcel.2004.05.027.

Ding B, Yue Y, Chen X, Long X, Zhou Z. Identification and expression analysis of miR396 and its target genes in Jerusalem artichoke under temperature stress. Gene. 2024;893:147908. https://doi.org/10.1016/j.gene.2023.147908.

Wang W, Cheng M, Wei X, Wang R, Fan F, Wang Z, et al. Comprehensive evolutionary analysis of growth-regulating factor gene family revealing the potential molecular basis under multiple hormonal stress in gramineae crops. Front Plant Sci. 2023;14. https://doi.org/10.3389/fpls.2023.1174955.

Wang H, Zhang Y, Liang D, Zhang X, Fan X, Guo Q, et al. Genome–wide identification and characterization of miR396 family members and their target genes GRF in sorghum (Sorghum bicolor (L.) moench). PLoS ONE. 2023;18:e0285494. https://doi.org/10.1371/journal.pone.0285494.

Gao Z, Ma C, Zheng C, Yao Y, Du Y. Advances in the regulation of plant salt-stress tolerance by MiRNA. Mol Biol Rep. 2022;49:5041–55. https://doi.org/10.1007/s11033-022-07179-6.

Liu D, Song Y, Chen Z, Yu D. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol Plant. 2009;136:223–36. https://doi.org/10.1111/j.1399-3054.2009.01229.x.

Zhang B, Tong Y, Luo K, Zhai Z, Liu X, Shi Z, et al. Identification of GROWTH-REGULATING FACTOR transcription factors in lettuce (Lactuca sativa) genome and functional analysis of LsaGRF5 in leaf size regulation. BMC Plant Biol. 2021;21:485. https://doi.org/10.1186/s12870-021-03261-6.

Yuan S, Zhao J, Li Z, Hu Q, Yuan N, Zhou M, et al. MicroRNA396-mediated alteration in plant development and salinity stress response in creeping bentgrass. Hortic Res. 2019;6:1–13. https://doi.org/10.1038/s41438-019-0130-x.

Yuan H, Cheng M, Wang R, Wang Z, Fan F, Wang W, et al. miR396b/GRF6 module contributes to salt tolerance in rice. Plant Biotechnol J. 2024;22:2079–92. https://doi.org/10.1111/pbi.14326.

Yan J, Qiu R, Wang K, Liu Y, Zhang W. Enhancing alfalfa resistance to spodoptera herbivory by sequestering microRNA396 expression. Plant Cell Rep. 2023;42:805–19. https://doi.org/10.1007/s00299-023-02993-z.

Cheng H, Song X, Hu Y, Wu T, Yang Q, An Z, et al. Chromosome-level wild Hevea Brasiliensis genome provides new tools for genomic-assisted breeding and valuable loci to elevate rubber yield. Plant Biotechnol J. 2023;21:1058–72. https://doi.org/10.1111/pbi.14018.

Jahan MS, Haider MM, Rahman M, Biswas D, Misbahuddin M. Chemical pulping: evaluation of rubber wood (Hevea brasiliensis) as a Raw material for kraft pulping. Nord Pulp Pap Res J. 2011;26:258–62. https://doi.org/10.3183/npprj-2011-26-03-p258-262.

Teoh YP, Don MM, Ujang S. Assessment of the properties, utilization, and preservation of rubberwood (Hevea brasiliensis): a case study in Malaysia. J Wood Sci. 2011;57:255–66. https://doi.org/10.1007/s10086-011-1173-2.

Severo ETD, Calonego FW, Sansígolo CA, Bond B. Changes in the chemical composition and decay resistance of Thermally-Modified Hevea Brasiliensis wood. PLoS ONE. 2016;11:e0151353. https://doi.org/10.1371/journal.pone.0151353.

Meng X, Kong L, Zhang Y, Wu M, Wang Y, Li J, et al. Gene expression analysis revealed Hbr-miR396b as a key piece participating in reaction wood formation of Hevea Brasiliensis (rubber tree). Ind Crops Prod. 2022;177:114460. https://doi.org/10.1016/j.indcrop.2021.114460.

Chen J, Liu M, Meng X, Zhang Y, Wang Y, Jiao N, et al. Multiomics studies with co-transformation reveal MicroRNAs via miRNA-TF-mRNA network participating in wood formation in Hevea Brasiliensis. Front Plant Sci. 2023;14. https://doi.org/10.3389/fpls.2023.1068796.

Liebsch D, Palatnik JF. MicroRNA miR396, GRF transcription factors and GIF co-regulators: a conserved plant growth regulatory module with potential for breeding and biotechnology. Curr Opin Plant Biol. 2020;53:31–42. https://doi.org/10.1016/j.pbi.2019.09.008.

Kim JH, Tsukaya H. Regulation of plant growth and development by the GROWTH-REGULATING FACTOR and GRF-INTERACTING FACTOR duo. J Exp Bot. 2015;66:6093–107. https://doi.org/10.1093/jxb/erv349.

Bazin J, Khan GA, Combier J-P, Bustos-Sanmamed P, Debernardi JM, Rodriguez R, et al. miR396 affects mycorrhization and root meristem activity in the legume edicago truncatula. Plant J. 2013;74:920–34. https://doi.org/10.1111/tpj.12178.

Miao C, Wang D, He R, Liu S, Zhu J-K. Mutations in 396e and 396f increase grain size and modulate shoot architecture in rice. Plant Biotechnol J. 2020;18:491–501. https://doi.org/10.1111/pbi.13214.

Li A-L, Wen Z, Yang K, Wen X-P. Conserved miR396b-GRF regulation is involved in abiotic stress responses in Pitaya (Hylocereus polyrhizus). Int J Mol Sci. 2019;20:2501. https://doi.org/10.3390/ijms20102501.

Kim JH, Choi D, Kende H. The AtGRF family of putative transcription factors is involved in leaf and cotyledon growth in Arabidopsis. Plant J. 2003;36:94–104. https://doi.org/10.1046/j.1365-313X.2003.01862.x.

Ma J-Q, Jian H-J, Yang B, Lu K, Zhang A-X, Liu P, et al. Genome-wide analysis and expression profiling of the GRF gene family in oilseed rape (Brassica Napus L). Gene. 2017;620:36–45. https://doi.org/10.1016/j.gene.2017.03.030.

Jodder J. Regulation of pri-MIRNA processing: mechanistic insights into the MiRNA homeostasis in plant. Plant Cell Rep. 2021;40:783–98. https://doi.org/10.1007/s00299-020-02660-7.

Jayashree R, Rekha K, Venkatachalam P, Uratsu SL, Dandekar AM, Kumari Jayasree P, et al. Genetic transformation and regeneration of rubber tree (Hevea Brasiliensis Muell. Arg) Transgenic plants with a constitutive version of an anti-oxidative stress superoxide dismutase gene. Plant Cell Rep. 2003;22:201–9. https://doi.org/10.1007/s00299-003-0666-x.

Udayabhanu J, Huang T, Xin S, Cheng J, Hua Y, Huang H. Optimization of the transformation protocol for increased efficiency of genetic transformation in Hevea Brasiliensis. Plants (Basel). 2022;11:1067. https://doi.org/10.3390/plants11081067.

Guo Z, Kuang Z, Zhao Y, Deng Y, He H, Wan M, et al. PmiREN2.0: from data annotation to functional exploration of plant MicroRNAs. Nucleic Acids Res. 2022;50:D1475–82. https://doi.org/10.1093/nar/gkab811.

Crooks GE, Hon G, Chandonia J-M, Brenner SE. WebLogo: a sequence logo generator. Genome Res. 2004;14:1188–90. https://doi.org/10.1101/gr.849004.

Chen C, Wu Y, Li J, Wang X, Zeng Z, Xu J, et al. TBtools-II: A one for all, all for one bioinformatics platform for biological big-data mining. Mol Plant. 2023;16:1733–42. https://doi.org/10.1016/j.molp.2023.09.010.

Tian F, Yang D-C, Meng Y-Q, Jin J, Gao G. PlantRegMap: charting functional regulatory maps in plants. Nucleic Acids Res. 2020;48:D1104–13. https://doi.org/10.1093/nar/gkz1020.

Chou K-C, Shen H-B. Cell-PLoc: a package of web servers for predicting subcellular localization of proteins in various organisms. Nat Protoc. 2008;3:153–62. https://doi.org/10.1038/nprot.2007.494.

Thompson JD, Gibson TJ, Higgins DG. Multiple sequence alignment using ClustalW and ClustalX. Curr Protoc Bioinf. 2002. https://doi.org/10.1002/0471250953.bi0203s00. Chap. 2:Unit 2.3.

Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4. https://doi.org/10.1093/molbev/msw054.

Wang J, Chitsaz F, Derbyshire MK, Gonzales NR, Gwadz M, Lu S, et al. The conserved domain database in 2023. Nucleic Acids Res. 2023;51:D384–8. https://doi.org/10.1093/nar/gkac1096.

Chao J, Wu S, Shi M, Xu X, Gao Q, Du H, et al. Genomic insight into domestication of rubber tree. Nat Commun. 2023;14. https://doi.org/10.1038/s41467-023-40304-y.

Cheng H, Tang C, Huang H. The Reyan 7-33-97 rubber tree genome: insight into its Structure, composition and application. In: Matsui M, Chow K-S, editors. The rubber tree genome. Cham: Springer International Publishing; 2020. pp. 13–40. https://doi.org/10.1007/978-3-030-42258-5_2.