Jinhui Chen1,2*
Mingming Liu1,3
Xiangxu Meng1,2Yuanyuan Zhang4,5
Yue Wang1,2Nanbo Jiao1,2Jianmiao Chen1,3- 1Sanya Nanfan Research Institute of Hainan University, Hainan Yazhou Bay Seed Laboratory/Key Laboratory of Genetics and Germplasm Innovation of Tropical Special Forest Trees and Ornamental Plants, Ministry of Education, School of Forestry, Hainan University, Sanya, China
- 2Engineering Research Center of Rare and Precious Tree Species in Hainan Province, School of Forestry, Hainan University, Haikou, China
- 3School of Tropical Crops, Hainan University, Haikou, China
- 4Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences, Haikou, Hainan, China
- 5State Centre for Rubber Breeding, Haikou, Hainan, China
Introduction: MicroRNAs (miRNAs) are small endogenous non-coding RNAs that play an important role in wood formation in plants. However, the significance of the link between miRNAs and their target transcripts in wood formation remains unclear in rubber tree (Hevea brasiliensis).
Methods: In this study, we induced the formation of reaction wood by artificially bending rubber trees for 300 days and performed small RNA sequencing and transcriptome deep sequencing (RNA-seq) to describe the complement of miRNAs and their targets contributing to this process.
Results and discussion: We identified 5, 11, and 2 differentially abundant miRNAs in normal wood (NW) compared to tension wood (TW), in NW relative to opposite wood (OW), and between TW and OW, respectively. We also identified 12 novel miRNAs and 39 potential miRNA-mRNA pairs with different accumulation patterns in NW, TW, and OW. We noticed that many miRNAs targeted transcription factor genes, which were enriched in KEGG pathways associated with phenylpropanoid biosynthesis, phenylalanine metabolism, and pyruvate metabolism. Thus, miRNA-TF-mRNA network involved in wood formation via tension wood model were constructed. We validated the differential accumulation of miRNAs and their targets by RT-qPCR analysis and overexpressed miRNA in Nicotiana benthamiana with its potential target gene. These results will provide a reference for a deep exploration of growth and development in rubber tree.
Introduction
Rubber tree (Hevea brasiliensis) is the main source of natural rubber (NR). Once latex production is no longer economically viable, rubber trees are also used as timber. The wood from rubber trees has in fact become the main export commodity in southeast Asia (Nair, 2010). Its delicate color and outstanding physical performance make it an excellent option for flooring and home furnishings. Due to its high potential commercial value, increasing timber yield and quality has become a key point of biotechnology in the rubber tree industry (Priyadarshan, 2017). However, two major limitations hinder a more widespread use of rubber wood: (i) The extent of unlignified or partially lignified tension wood fiber, which is not easily digested by enzymes, is high, and the proportion of normal fibers is low (Mellerowicz and Gorshkova, 2012; Gritsch et al., 2015); and (ii) it has high sensitivity to biodegradation owing to low levels of phenolic compounds with biocidal activity (Pramod et al., 2019; Thaochan et al., 2020). The biosynthesis of lignin and polyphenolic derivatives in living trees, particularly in rapidly growing woody plants such as rubber trees, contributes to wood quality and durability (Kuyyogsuy et al., 2018; Pramod et al., 2019; Thaochan et al., 2020).
Lignin is an abundant biopolymer that is essential for plant cell wall integrity and stem strength (Shen et al., 2012). The biosynthesis of lignin monomers begins with phenylalanine deamination, leading to the production of three monolignin alcohols: coniferyl, sinapyl, and p-coumaryl alcohols (Eudes et al., 2012). Several genes that participate in lignin biosynthesis and modulate lignin levels have been identified in dicots (Sibout et al., 2005; Weng et al., 2010). For instance, knocking down 4-coumarate:CoA ligase (4CL) expression in hybrid poplar (Populus tremula × P. alba) resulted in a sharp decrease in lignin content and markedly altered wood chemical composition and wood metabolism (Voelker et al., 2010). The PtrWND (wood-associated NAC domain) genes were shown to induce the expression of wood biosynthetic genes, including associated structural genes and transcription factors, resulting in ectopic deposition of lignin in poplars (Zhong et al., 2010). Despite the above progress in our understanding of lignin biosynthesis, much remains to be investigated in terms of transcriptional and post-transcriptional regulation. Recently, regulation of wood formation via non-coding RNAs (ncRNAs) and microRNAs (miRNAs) has received increasing attention (Lu et al., 2013).
miRNAs are post-transcriptional modulators of gene function by promoting the cleavage of their complementary target messenger RNAs (mRNAs) and/or imposing translation repression (Bartel, 2004; Zhang et al., 2019). Several studies have illustrated the vital roles played by miRNAs in wood formation (Li et al., 2020; Yu et al., 2020). For instance, transcript levels for 17 of the 29 LACCASE (PtrLAC) genes in the black cottonwood (P. trichocarpa) genome decreased in P. trichocarpa trees overexpressing miR397a, in turn leading to a reduction of lignin levels (Lu et al., 2013). Similarly, the overexpression of miR319a in Populus tomentosa in seedlings resulted in delayed secondary growth and decreased xylem production (Hou et al., 2020). In particular, miR165b guides the development of pith secondary cytoderm is by restraining the AtHB15 (HOMEOBOX 15) expression domain (Du and Wang, 2015). However, how wood anatomical features like container shape and thickness are genetically governed is not very clear. These traits are important for cell wall composition and overall tree performance (Quan et al., 2018; Quan et al., 2019).
To explore the molecular basis of changes in transcript levels caused by miRNAs during wood formation, we sequenced small RNAs from three wood parts: normal wood (NW)which is from the stem of the tree at breast height, tension wood (TW) which is upper side of the bending trunk, and opposite wood (OW) which is lower side of the bending trunk. We then assigned predicted target genes to these wood-abundant miRNAs by identifying genes whose transcript levels were inversely correlated with miRNA abundance. We complemented this approach by using a bioinformatics method for miRNA target prediction and constructed the resulting miRNA–mRNA interaction network. The small RNAs found in this study are good candidates for miRNAs that are involved in wood formation. They may also help with the development of functional markers for molecular breeding in rubber trees and other tropical plants to help change the composition of lignin or physical characteristics.
Materials and methods
Plant materials and microscopy observations
The rubber trees (clone Reyan 7-33-97) used in this study were grown in the experimental greenhouse of Hainan University (Danzhou, Hainan, China; 109°29′25′′ E, 19°30′40′′ N) at the end of June 2016. To probe the genes involved in the formation of reaction wood, three rubber trees of similar age and with trunks of similar diameter (about 2 cm) were bent at a 30° angle for 300 days and were selected as experimental materials to force the formation of reaction wood (Figure S1). The bending was applied starting at 9 a.m. on August 17th, 2020, and ended at 9 a.m. on June 13rd, 2021, at which point the wood samples were rapidly processed. The wood quality from the collected samples was assessed by scanning electron microscopy (SEM, Phenom proX, the Netherlands) in Center for Analytical Instrumentation (Hainan University), and the resulting images were processed in ImageJ software to measure the area of the gelatinous (G) layer. The SEM images indicated that the TW (tension wood) reaction wood had an extremely dense cementitious layer (G-layer; Figure S2A) that was not present in NW (normal wood) or OW (opposite wood; Figures S2B, C). These observations were consistent with earlier findings (Sujan et al., 2015). Therefore, xylem samples from the collected trees were selected for further analysis.
Samples were collected for NW, TW, and OW from the same individuals to allow for direct comparison in an identical genetic background. Briefly, the bark above the sampling area was removed to expose the inner wood layers. A sharp razor blade was then used to collect TW (upper side) and OW (lower side) from the same branch (Li et al., 2013). The control for stem xylem tissue was NW and was collected about 1 m above the ground, before the bending point. All collected samples were about 2 cm × 1 cm and 4–5 mm in depth. Samples were harvested in the morning, quickly frozen in liquid nitrogen, and stored at –80°C until use.
RNA extraction and qualification
Total RNA was extracted from nine samples (NW1, NW2, NW3, OW1, OW2, OW3, TW1, TW2, and TW3) using a modified cetyl trimethyl ammonium bromide (CTAB) method (Chang et al., 1993). Each tree counted as one biological replicate (NW1-3, TW1-3, and OW1-3). Traces of genomic DNA were removed with RNase-free DNase I digestion (Takara, Beijing, China). RNA degradation and DNA contamination were assessed by electrophoresis on 1% (w/v) agarose gels. RNA concentration, quality, and integrity were estimated on a NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and an RNA Nano 6000 Chip on a Bioanalyzer 2100 (Agilent Technologies, CA, USA). Only RNA samples with an OD260/280 ratio of 1.9–2.2, an OD260/230 ratio ≥ 2.0, and RNA integrity number (RIN) values > 6.8 were processed for further experiments.
Small RNA library construction and sequencing
Three micrograms of total RNA per sample was used to construct sequencing libraries with a NEBNext multiplex small RNA library prep kit for Illumina (NEB, USA) following the manufacturer’s protocol. Briefly, the NEB 3′ SR adapter was ligated to the 3′ end of miRNAs, siRNAs (small interfering RNAs), and piRNAs (piwi-interacting RNAs). The SR real-time primer was then annealed to the 3′ SR adapter to initiate double-stranded DNA (dsDNA). The 5′ end adapter was connected to the 5′ ends of miRNAs. First-strand cDNA synthesis with M-MuLV Reverse Transcriptase (RNase H-). The libraries were amplified by 35 PCR cycles with index (X) primer, SR primer for Illumina, and LongAmp Taq 2X master mix. The PCR products were separated on 8% polyacrylamide gel for 80 min at 100 V. DNA fragments responding to 140-160 bp (the length of sRNAs with 5′ and 3′ adapters) were purified from the gel and eluted in 8 µL of elution buffer. Library quality and titer were evaluated using a DNA high sensitivity chip and Agilent Bioanalyzer 2100 instrument. A TruSeq SR Cluster Kit v3-cBot-HS (Illumina) was used for cluster formation on a cBot cluster generation system following the manufacturer’s instructions. Libraries were sequenced on an Illumina HiSeq 2,500 platform as 50-bp single-end reads.
Identification of known miRNAs and novel miRNAs
After adapter trimming from the raw reads, any clean reads shorter than 18 nt or low-quality reads (reads with N bases > 10% and reads with a 3′ end with Q < 20 [Q = –10Log10error_ratio]) were discarded. The remaining clean reads were mapped to the rubber tree reference genome (Liu et al., 2020) with Bowtie2 allowing no mismatch (Langmead et al., 2009). Reads derived from rRNAs, protein-coding genes, snRNAs, tRNAs, snoRNAs, and repeat sequences were removed by filtering the clean reads against the Rfam database and RepeatMasker (Tempel, 2012; Kalvari et al., 2018). Known miRNAs were identified with miRBase 20.0 (Meng et al., 2018); potential miRNAs and their secondary structures were determined using the miRDeep2 algorithm (Friedländer et al., 2012) and sRNA-tools-cli (http://srna-workbench.cmp.uea.ac.uk/). The formation of a typical hairpin structure was used as a criterion to identify novel miRNAs from the transcripts showing no match to known miRNAs. Unannotated small RNAs were assessed with miRDeep2 (Friedländer et al., 2012) and miREvo (Wen et al., 2012) using minimum free energy, possible Dicer cleavage sites, and the secondary structures of small RNA tags. Custom scripts were applied to count all candidate miRNAs and estimate the basic deviation between each position of all identified miRNAs and the first position of identified miRNAs of a given length. The secondary structure of the novel_28 mature sequence was predicted and compared to other miRNA families from rubber tree and other species.
Prediction of the target genes of miRNAs
The RNA-seq data were from our previous study (Meng et al., 2021). The online tool psRNATarget (https://www.zhaolab.org/psRNATarget/) was used to predict miRNA targets with default parameters with expectation ≤ 3 (Dai et al., 2018). A gene was deemed to be a putative target for a miRNA when its transcript levels showed a negative Pearson’s correlation coefficient with miRNA abundance (correlation < –0.8, P-value < 0.05). miRNA abundance was estimated with the formula (Zhou et al., 2010): Normalized abundance = mapped read count/total reads * 1,000,000.
KEGG enrichment analysis of co-expressed target genes
The predicted co-expressed target genes were subjected to KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway. KEGG (http://www.genome.jp/kegg/) pathway enrichment analysis was carried out with KOBAS software (Mao et al., 2005).
Construction for miRNA–mRNA interaction networks and tissue-specific expression analysis
The Pearson’s correlation coefficients between miRNAs and transcription factor genes were calculated using expression values from this study for miRNAs and from our previous study (Meng et al., 2021) for mRNAs or between transcription factor genes and other genes to identify miRNAs, transcription factor genes, and their co-expressed genes. The resulting interaction network was built in R (version 4.0.1) and visualized with Cytoscape (version 3.8.1).
Plasmid construction
According to the target sites (Figure 1A), to overexpress novel_28 (hbr-miR482c), the pre-miR482c sequence was amplified from genomic DNA isolated from normal wood tree with primers containing BamH I and Sac I restriction sites. The resulting PCR product was digested with BamH I and Sac I and ligated into the vector pBI121 downstream of the cauliflower mosaic virus (CaMV) 35S promoter (Shanghai Generay Biotech; Figure 1B). The full-length coding sequence of HbrCAD1 (GH714_013930) was also cloned into pBI121 (GH714_013930-pBI121; Figure 1B). Both constructs were transformed into Agrobacterium (Agrobacterium tumefaciens) strain GV3101. The primer sequences used for cloning are listed in Table S1.
Figure 1 Validation of the targeting of HbrCAD1 by its predicted miRNA novel_28 in Nicotiana benthamiana. β-actin was selected as internal reference; data are means ± standard error of three independent biological replicates. (A) Alignment of hbr-miR482c and HbrCAD1. (B) The plasmids pBI121-miR482c and GH714-013930-pBI121 used for the assay. (C) Principle of transient infiltration of N. benthamiana leaves with an Agrobacterium cell suspension. (D) Relative HbrCAD1 expression levels in different samples. Data are shown as Log2 (fold-change), with the expression of HbrCAD1 from the pair pBI121 + 35S:HbrCAD1 set to 1.
Verification of interaction between hbr-miR482c and its target
Agrobacterium-mediated transient expression in Nicotiana benthamiana leaves (English, 1997) was used to assess the targeting of HbrCAD1 transcripts by hbr-miR482c. Agrobacterium cultures harboring the hbr-miR482c or HbrCAD1 construct were resuspended in infiltration buffer, mixed, and infiltrated into six N. benthamiana leaves (Figure 1C) from 3-week-old plants. As negative controls, Agrobacterium cultures a mix of cultures harboring pB121 and the 35S::HbrCAD1 construct were infiltrated in N. benthamiana leaves.
Validation of miRNA expression and that of their target genes by RT-qPCR
A GoScript™ Reverse Transcription kit (Promega, USA) was used for first-strand cDNA synthesis from total RNA of the nine samples collected in this study. Six differentially expressed miRNAs (novel_47, novel_67, novel_28, novel_165, novel_177, and novel_101) and six target genes (ALDO1, PEPcase, CAD1, RF2, PKc_like, GT) of miRNA-mRNA correlation network were chosen (Figure S3). Gene-specific primers were designed for the target genes with Primer Premier v5 software (Table S1). qPCR was performed with TB Green Premix Ex Taq II (Tli RNase H Plus; Takara, Beijing, China) for the target genes according to the manufacturer’s instructions. PCR conditions were as follows: denaturation at 94°C for 2 min, then 40 cycles of 95°C for 5 s and 60°C for 30 s. Ubiquitin was used as internal reference for normalization of expression data of rubber tree samples (Meng et al., 2022), and β-actin was used for N. benthamiana (Nawaz et al., 2019). For miRNA, a miRNA RT-qPCR Detection Kit (Aidlab, Beijing, China) was used following the manufacturer’s instructions. PCR conditions were as above. U6 transcripts were used as internal reference for normalization (Zeng et al., 2009). A melting curve was performed from 60°C to 95°C to confirm the specificity of the amplicons. Relative expression levels of miRNAs and their target genes were estimated by the 2–△△Ct method (Schmittgen and Livak, 2008). Three technical replicates were analyzed per sample.
Results
Overview of small RNA sequencing from reaction wood
We constructed nine small RNA libraries from different wood tissues (NW, TW, and OW) to explore the role of miRNAs in wood development in rubber trees. We obtained between 10.15 and 14.50 million raw reads after sequencing. We removed low-quality reads and removed adapters to yield 9.86 to 14.35 million clean reads with a length ranging from 18 to 30 nucleotides (nt) (Table 1).
Table 1 Summary of reads from small RNA sequencing.
We aligned the clean reads against the rubber tree reference genome with Bowtie2 and used RSEM software to estimate read counts per gene model. We successfully mapped from 5,936,751 to 9,992,003 reads to the rubber tree genome, or 89.53%-93.80% of all clean reads across the nine small RNA libraries (Table 1). We then removed all reads mapping to tRNA (transfer RNA), snoRNA (small nucleolar RNA) and snRNA (small nuclear RNA) loci. We focused on the remaining unannotated reads. We observed a peak for 21-nt sRNAs, representing >14% of all candidates unannotated sRNAs in the libraries, with a second peak for 24-nt sRNAs (from 10% to 13%; Figure 2).
Figure 2 Length distribution of miRNA sequences sequenced from Hevea brasiliensis reaction wood.
Annotation of known and novel miRNAs expressed in reaction wood
We compared the unannotated sRNAs identified above to all miRNAs or their pre-miRNAs deposited in miRBase to identify known and novel miRNAs. We identified 22 (NW), 22 (TW), and 23 (OW) known miRNAs belonging to 21 miRNA families (Figure 3A). Of these, 19 miRNAs were shared across the NW, TW, and OW libraries (Figure 3B; Table S2). We determined that the first base of these known miRNAs tended to be a uracil (U) for 18–22-nt miRNAs (Figure S4).
Figure 3 Numbers of miRNAs identified in different xylem tissues. (A) Known and novel miRNAs identified from small RNA sequencing of NW, TW, and OW from 300-day rubber tree reaction wood. (B) Venn diagram showing the extent of overlap between known miRNAs in each tissue type (NW, TW, and OW) from 300-day rubber tree reaction wood.
The known miRNA hbr-miR166, hbr-miR396, hbr-miR408, and hbr-miR482 families were each represented by two members, while the remaining miRNA families only had one member (Figure S5). Furthermore, two miRNAs (hbr-miR6170 and hbr-miR6171) were specific to OW (Figure S5), suggesting that each wood tissue is associated with a slightly different set of miRNAs.
The sequences failed to match known miRNAs or pre-miRNAs were used to predict novel miRNA resulting in 89 novel miRNAs with typical hairpin structure (Figure 3, S6). The first base of these novel rubber tree miRNAs (18–30 nt in length) also largely started with a uracil (Figure S7). Furthermore, we cloned and assigned novel_28 to the miR482 family through BLAST search against other species in the miRBase database. We then predicted the secondary structure of novel_28 and compared its mature sequence to that of other hbr-miR482 family members, which revealed novel_28 is a new member of the hbr-miR482 family (Figures S6, S8). Thus, we renamed novel_28 as hbr-miR482c.
Differentially abundant miRNAs
To investigate how miRNAs might regulate wood formation, we quantified the abundance of all miRNAs as TPMs (transcripts per million), followed by a comparison of miRNA abundance between the tree wood types (|Fold change|≥2 and P ≤ 0.05). The OW vs. NW comparison yielded 11 differentially abundant miRNAs, 5 for the TW vs. NW comparison, and 2 for the TW vs. OW comparison. Of the five differentially abundant miRNAs between TW and NW samples, three were more abundant in TW and two were more abundant in NW tissues. Similarly, 5 miRNAs were more accumulated in OW and 6 were more accumulated in NW tissues. Finally, one miRNA was more abundant in each TW and OW tissue (Table S3).
Prediction of miRNA target genes
We predicted miRNA targets through psRNATarget (Dai et al., 2018), which identifies transcripts with complementary sequences to those of miRNA candidates. We independently calculated Pearson’s correlation coefficients between the abundance of each miRNA and the transcript levels of all rubber tree transcripts obtained from a previous study (Meng et al., 2021). We considered genes as candidate targets when their transcript levels were negatively correlated (<–0.8) with miRNA abundance. In the TW vs. NW comparison, we determined that 4 novel miRNAs have the potential to target 10 genes, but no clear target could be identified for the novel_86 miRNA. Similarly, we identified 31 targets for 8 of the differentially abundant miRNAs in the OW vs. NW comparison, with no predicted targets for the other three novel miRNAs (Table S4). Further, we identified 11 target genes for the 2 differentially abundant miRNAs from the comparison between TW and OW. We selected those potential targets with a negative correlation of –0.8 or below with their associated miRNA, resulting in 39 putative targets for 12 novel miRNAs (cor ≤ –0.8 and P ≤ 0.05) (Table S5).
Identification of genes co-expressed with transcription factor genes involved in rubber tree reaction wood formation
In the above list of target genes, we noticed that 53.85% (or 21 of 39) encode transcription factors (TFs; Figure S9). To define the putative targets of the encoded TFs, we measured the Pearson’s correlation coefficients (PCCs) between TF genes and all other rubber tree genes, using RNA-seq data available at the NCGC (National Genomics Data Center) under the accession numbers CRA004241 and CRA004243. A KEGG pathway enrichment analysis showed that these co-expressed targets are largely related to phenylpropanoid biosynthesis, phenylalanine metabolism, and fatty acid biosynthesis (Figure S10). These results suggest a central role for these TFs during reaction wood formation under prolonged mechanical stress.
Differentially expressed mRNA-miRNA pairs related to wood formation
We then explored how miRNA-mediated adjustment of transcripts affected reaction wood development. Here, we focused on target genes co-expressed with those TF genes that were associated with phenylalanine metabolism or carotenoid biosynthesis and constructed the underlying interaction network (Figure S10, 4). Our above analysis predicted that hbr-miR482c modulates the transcript levels of HbrCAD1, which is related to phenylpropanoid biosynthesis. Similarly, novel_67 might modulate the transcript levels of genes related to lignin biosynthesis by targeting C3H-type zinc finger TF family members (Figure 4); novel_93 was predicted to modulate the transcript levels of genes related to phenylalanine metabolism. In this network, hbr-miR482c and novel_76 each targeted two genes involved in wood formation. The C3H, FAR1 (FAR-RED IMPAIRED RESPONSE 1), bZIP, and MYB TF families possibly play vital functions in wood growth from our miRNA-TF-mRNA network. The target genes co-expressed with these TF genes, such as fatty acyl-CoA reductase 1 (FAR1, gene-GH714_004418), cinnamyl-alcohol dehydrogenase (CAD1), and p-coumarate 3-hydroxylase (C3H, gene-GH714_027993), were involved in carotenoid biosynthesis. We thus propose that the miRNA-TF-mRNA regulatory network described here may play a significant role in adjusting the molecular events necessary for reaction wood development in rubber tree. The expression levels of these genes and associated miRNAs are shown in Figure 5.
Figure 4 miRNA-transcription factor-mRNA networks associated with wood formation. Yellow octagons represent miRNAs, orange triangles represent transcription factor genes, and blue ellipses represent genes.
Figure 5 Expression profile of miRNAs and genes co-expressed with transcription factor genes. (A) Heatmap representation of the expression levels of genes co-expressed with transcription factor genes in different tissues. (B) Heatmap representation of the expression profile of miRNAs in TW, OW, and NW. Columns represent three different tissue types and rows represent different transcripts. Each square represents a transcript and the color indicates the level of expression; red represents up-regulation and blue represents down-regulation.
We wished to confirm the transcript levels of miRNAs and their target genes estimated from the RNA-seq datasets by a reverse transcription quantitative PCR (RT-qPCR) analysis of the same RNA samples. We observed a generally comparable trend in the accumulation of most miRNAs and their target genes between RT-qPCR and RNA-seq data in the OW, NW, and TW samples, although the fold-change (FC) values from RNA-seq did not precisely match the expression values obtained by RT-qPCR (Figure S3). These results indicate the dependability of the sequencing data.
HbrCAD1 transcripts are cleaved by hbr-miR482c
miRNAs can cause the cleavage of their target transcripts between the 10th and 11th nucleotides of their target site (Rhoades et al., 2002). In this study, we predicted from bioinformatics analysis that hbr-miR482c may influence lignin biosynthesis by targeting HbrCAD1 transcripts (Figure 4). To test this hypothesis experimentally, we co-expressed the precursor of hbr-miR482c and HbrCAD1, driven by the strong cauliflower mosaic virus (35S) promoter, into N. benthamiana leaves by Agrobacterium-mediated transient infiltration. As negative controls, we co-infiltrated the construct expressing the HbrCAD1 and the empty effector vector (pBI121). We observed that HCAD1 transcript levels decrease significantly when the hbr-miR482c precursor is co-expressed (Figure 1D). Thus we suspected that HbrCAD1 transcripts might be the target gene of hbr-miR482c.
Discussion
A model for reaction wood formation in rubber tree and other tree species
The scanning electron microscopy of the xylem samples from reaction wood had confirmed that the tension wood had the remarkably thick gelatinous layer (Figure S2). These results were supported by earlier findings and confirmed a tension wood model were constructed in our study (Sujan et al., 2015). We wished to identify genes important for wood formation based on their expression models in various wood tissues of rubber tree. We drew heatmaps and constructed an interaction network to reveal connections between genes related to phenylpropanoid biosynthesis.
Enriched KEGG terms underscored the involvement of phenylpropanoid biosynthesis, phenylalanine metabolism, and pyruvate metabolism in reaction wood formation. Indeed, our results showed that the genes co-expressed with miRNA-targeted TF genes were related to metabolism and physiological functions that all contribute to reaction wood development. Similar consequences have been reported in previous work, lending some support to our current study (Li et al., 2013; Lv et al., 2021). Future work should pay close attention to exploring the functions of noncoding RNAs and their candidate target genes predicted from our RNA-seq results. In particular, genes linked to wood development can now be identified from their expression patterns across wood growth regions.
Integrated analysis of the miRNA-TF-mRNA network
Important transcription factors associated with secondary growth have been identified, such as members of the C3H and MYB TF families (Demura and Fukuda, 2007; Zhong and Ye, 2009). Zhang et al. reported that miRNAs may be connected to tension wood development by regulating secondary cell wall biosynthesis in Moso bamboo (Phyllostachys edulis) (Zhang et al., 2018). In this study, we determined that the transcript levels of 21 TF genes from 13 families changed over the course of wood bending (Figure S9). Of these, TF genes from the C3H and MYB families may help adjust the expression of genes related to phenylpropanoid biosynthesis and possibly play a significant role in the wood development in rubber tree. Previous research revealed that overexpressing the MYB TF genes PtoMYB216, PtoMYB74, and PtoMYB92 from P. tomentosa induced the expression of genes related to lignin biosynthesis, resulting in thicker xylem cell walls, more xylem layers, ectopic lignin deposition, and enhanced lignin contents by 13–50% (Qiaoyan et al., 2013; Li et al., 2015; Li et al., 2018). Analogously, compared to the control, the expression levels of lignin biosynthetic genes, lignification ability, xylem volume, and lignin levels of C3H overexpressed in Arabidopsis all increased (Fornalé et al., 2015). We also observed that C3H-type zinc finger TF family members possibly take part in wood development by regulating the transcript levels of GLUCOSYLTRANSFERASE (GT) (gene-GH714_027993), which is associated with cellulose biosynthesis. This observation underscores the significance of C3H family members during plant development.
miRNA regulaged key genes in wood formation pathway
The proteins that are thought to catalyze glucan-chain elongation in cellulose and callose biosynthesis are processive GTs that belong to GT2 family. In Arabidopsis, 10 to 12 GT2 family members form CESA (cellulose synthase catalytic subunit) and callose synthase (Yang et al., 2016). Moreover, GTs are a vital component of wood development. For instance, in Arabidopsis, CesA8 participates in secondary cell wall formation, as its loss-of-function mutation produced plants with a delicate stem phenotype (Taylor et al., 1999; Taylor et al., 2000; Taylor et al., 2003). Similarly, a mutation in rice CesA8 caused a dramatic decrease in the cellulose content of secondary cell walls, resulting in a brittle culm phenotype (Zhang et al., 2009; Song et al., 2013). Meng et al. found that HbrCesA8, which is related to cellulose biosynthesis, may participate in reaction wood formation in rubber tree (Meng et al., 2021). These findings support our miRNA-mRNA model for rubber tree wood formation. Our multiomics-based approach produced a post transcription network of which several nodes are regulated by novel_67, which possibly affects some target genes during wood development.
The expression levels of CAD, 4CL, peroxidase (POD), and CAFFEATE O-METHYLTRANSFERASE (COMT) are specifically linked to lignin composition (Do et al., 2007; Wagner et al., 2009; Vanholme et al., 2010; Chanoca et al., 2019). Here, we showed that HbrCAD1 transcript levels are less abundant in OW tissues compared to NW, suggesting that HbrCAD1-catalyzed reactions might contribute less to wood formation in OW relative to NW tissues. Previous studies have shown that repressing PAL, CAD, or other enzymes of the lignin biosynthetic pathway may cause decreased lignin content (Chanoca et al., 2019). Furthermore, in Arabidopsis, studies have indicated that single or double knockout mutants in PODs representatively resulted in a small but significant decrease in lignin content and altered lignin biosynthesis in inflorescence stalks, indicating that modulating POD expression levels may affect lignin composition (Herrero et al., 2013; Barros et al., 2015). In general, we hypothesize that the expression pattern of HbrCAD1 contributes to establishing the difference in lignin levels across tissue types.
Conclusions
Based on the small RNA data obtained from wood tissues collected from rubber tree, we identified 114 miRNAs (25 known and 89 novel) present in 300-day reaction wood. We also established a network linking miRNAs, their putative TF target genes, and the genes that are co-expressed with these TF genes in the context of cellulose biosynthesis. Finally, we revealed the interaction landscape of these three regulatory layers in adjusting reaction wood growth and validated the network in wood formation of rubber trees. In summary, we described target genes associated with wood development in rubber tree and studied their post-transcriptional regulation. These results will provide the theoretical basis to clarify miRNA-mediated post-transcriptional mechanisms during wood growth and development in rubber trees.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://bigd.big.ac.cn/gsa, CRA007612 https://bigd.big.ac.cn/gsa, CRA004818.
Author contributions
JHC designed the research. MML, XM, YZ, YW, NJ and JMC performed the research. All authors analyzed and interpreted the data. JHC and MML wrote the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation of China (31960321), the PhD Scientific Research and Innovation Foundation of Sanya Yazhou Bay Science and Technology City (HSPHDSRF-2023-12-006) and the Scientific Research Fund Project of Hainan University (KYQD (ZR)1830).
Acknowledgments
We thanks Center for Analytical Instrumentation, Hainan University, for their help in our experiments. This research work is also supported by High-performance Computing Platform of YZBSTCACC.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1068796/full#supplementary-material
References
Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 116 (2), 281–297. doi: 10.1016/s0092-8674(04)00045-5
Barros, J., Serk, H., Granlund, I., Pesquet, E. (2015). The cell biology of lignification in higher plants. Ann. Bot. 115, 1053–1074. doi: 10.1093/aob/mcv046
Chang, S., Puryear, J., Cairney, J. (1993). A simple and efficient method for isolating RNA from pine trees. Plant Mol. Biol. Rep. 11 (2), 113–116. doi: 10.1007/BF02670468
Chanoca, A., de, V. L., Boerjan, W. (2019). Lignin engineering in forest trees. Front. Plant Sci. 10. doi: 10.3389/fpls.2019.00912
Dai, X., Zhuang, Z., Zhao, P. X. (2018). PsRNATarget: a plant small RNA target analysis server, (2017 Release). Nucleic Acids Res. 46 (W1), W49–W54. doi: 10.1093/nar/gky316
Demura, T., Fukuda, H. (2007). Transcriptional regulation in wood formation. Trends Plant Sci. 12, 64–70. doi: 10.1016/j.tplants
Do, C. T., Pollet, B., Thévenin, J., Sibout, R., Denoue, D., Barrière, Y., et al. (2007). Both caffeoyl coenzyme a 3-o-methyltransferase 1 and caffeic acid O-methyltransferase 1 are involved in redundant functions for lignin, flavonoids and sinapoyl malate biosynthesis in Arabidopsis. Planta. 226 (5), 1117–1129. doi: 10.1007/s00425-007-0558-3
Du, Q., Wang, H. (2015). The role of HD-ZIP III transcription factors and miR165/166 in vascular development and secondary cell wall formation. Plant Signal Behav. 10 (10), e1078955. doi: 10.1080/15592324.2015.1078955
English, J. J. (1997). Requirement of sense transcription for homology-dependent virus resistance and trans -inactivation. Plant J. 12, 597–603. doi: 10.1046/j.1365-313X.1997.d01-13.x
Eudes, A., George, A., Mukerjee, P., Kim, J. S., Pollet, B., Benke, P. I., et al. (2012). Biosynthesis and incorporation of side-chain-truncated lignin monomers to reduce lignin polymerization and enhance saccharification. Plant Biotechnol. J. 10 (5), 609–620. doi: 10.1111/j.1467-7652.2012.00692.x
Fornalé, S., Rencoret, J., Garcia-Calvo, L., Capellades, M., Encina, A., Santiago, R., et al. (2015). Cell wall modifications triggered by the down-regulation of coumarate 3-hydroxylase-1 in maize. Plant Sci. 236, 272–282. doi: 10.1016/j.plantsci.2015.04.007
Friedländer, M. R., MacKowiak, S. D., Li, N., Chen, W., Rajewsky, N. (2012). MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40 (1), 37–52. doi: 10.1093/nar/gkr688
Gritsch, C., Wan, Y., Mitchell, R. A. C., Shewry, P. R., Hanley, S. J., Karp, A. (2015). G-Fibre cell wall development in willow stems during tension wood induction. J. Exp. Bot. 66, 6447–6459. doi: 10.1093/jxb/err339
Herrero, J., Fernández-Pérez, F., Yebra, T., Novo-Uzal, E., Pomar, F., Pedreño, M. Á., et al. (2013). Bioinformatic and functional characterization of the basic peroxidase 72 from arabidopsis thaliana involved in lignin biosynthesis. Planta 237 (6), 1599–1612. doi: 10.1007/s00425-013-1865-5
Hou, J., Xu, H., Fan, D., Ran, L., Li, J., Wu, S., et al. (2020). Mir319a-targeted ptotcp20 regulates secondary growth via interactions with ptowox4 and ptownd6 in populus tomentosa. New Phytol. 228 (4), 1354–1368. doi: 10.1111/nph
Kalvari, I., Nawrocki, E. P., Argasinska, J., Quinones-Olvera, N., Finn, R. D., Bateman, A., et al. (2018). Non-coding RNA analysis using the rfam database. Curr. Protoc. Bioinf. 62 (1), e51. doi: 10.1002/cpbi.51
Kuyyogsuy, A., Deenamo, N., Khompatara, K., Ekchaweng, K. (2018). Chitosan enhances resistance in rubber tree (Hevea brasiliensis), through the induction of abscisic acid (ABA). Physiol. Mol. Plant P. 102, 67–78. doi: 10.1016/j.pmpp.2017.12.001
Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10 (3), R25. doi: 10.1186/gb-2009-10-3-r25
Li, H., Huang, X., Li, W., Lu, Y., Dai, X., Zhou, Z., et al. (2020). MicroRNA comparison between poplar and larch provides insight into the different mechanism of wood formation. Plant Cell Rep. 39 (9), 1199–1217. doi: 10.1007/s00299-020-02559-3
Li, X., Yang, X., Wu, H. X. (2013). Transcriptome profiling of radiata pine branches reveals new insights into reaction wood formation with implications in plant gravitropism. BMC Genomics 14, 768. doi: 10.1186/1471-2164-14-768
Li, C., Wang, X., Ran, L., Tian, Q., Fan, D., Luo, K. (2015). Ptomyb92 is a transcriptional activator of the lignin biosynthetic pathway during secondary cell wall formation in populus tomentosa. Plant Cell Physiol. 56 (12), 2436–2446. doi: 10.1093/pcp/pcv157
Li, C., Ma, X., Yu, H., Fu, Y., Luo, K. (2018). Ectopic expression of ptomyb74 in poplar and arabidopsis promotes secondary cell wall formation. Front. Plant Sci. 9, 1262 doi: 10.3389/fpls.2018.01262
Liu, J., Shi, C., Shi, C. C., Li, W., Zhang, Q. J., Zhang, Y., et al. (2020). The chromosome-based rubber tree genome provides new insights into spurge genome evolution and rubber biosynthesis. Mol. Plant 13 (2), 336–350. doi: 10.1016/j.molp.2019.10.017
Lu, S., Li, Q., Wei, H., Chang, M. J., Tunlaya-Anukit, S., Kim, H., et al. (2013). Ptr-miR397a is a negative regulator of laccase genes affecting lignin content in Populus trichocarpa. Proc. Natl. Acad. Sci. U.S.A. 110 (26), 10848–10853. doi: 10.1073/pnas.1308936110
Lv, C., Lu, W., Quan, M., Xiao, L., Li, L., Zhou, J., et al. (2021). Pyramiding superior haplotypes and epistatic alleles to accelerate wood quality andyield improvement in poplar breeding. Ind. Crop Prod. 171, 113891. doi: 10.1016/j.indcrop.2021.113891
Mao, X., Cai, T., Olyarchuk, J. G., Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics. 21, 3787–3793. doi: 10.1093/bioinformatics/bti430
Mellerowicz, E. J., Gorshkova, T. A. (2012). Tensional stress generation in gelatinous fibres: a review and possible mechanism based on cell-wall structure and composition. J. Exp. Bot. 63, 551–565. doi: 10.1093/jxb/err339
Meng, X. X., Kong, L. S., Zhang, Y. Y., Wu, M. J., Wang, Y., Li, J., et al. (2022). Gene expression analysis revealed hbr-miR396b as a key piece participating in reaction wood formation of Hevea brasiliensis (rubber tree). Ind. Crop Prod. 177, 114460. doi: 10.1016/j.indcrop.2021.114460
Meng, X. X., Wang, Y., Li, J., Jiao, N., Zhang, X., Zhang, Y., et al. (2021). RNA Sequencing reveals phenylpropanoid biosynthesis genes and transcription factors for Hevea brasiliensis reaction wood formation. Front. Genet. 29. doi: 10.3389/fgene.2021.763841
Meng, X., Zhang, X., Li, J., Liu, P. (2018). Identification and comparative profiling of ovarian and testicular microRNAs in the swimming crab Portunus trituberculatus. Gene. 640, 6–13. doi: 10.1016/j.gene
Nair, K. (2010). The agronomy and economy of important tree crops of the developing world. Agron. Economy Important Tree Crops Developing World. 30 (4), 313–351. doi: 10.1016/B978-0-12-384677-8.00008-4
Nawaz, Z., Kakar, K. U., Ullah, R., Yu, S., Zhang, J., Shu, Q. Y., et al. (2019). Genome-wide identification, evolution and expression analysis of cyclic nucleotide-gated channels in tobacco (Nicotiana tabacum l.). Genomics. 111 (2), 142–158. doi: 10.1016/j.ygeno.2018.01.010
Pramod, S., Reghu, C. P., Rao, K. S. (2019). Biochemical characterization of wood lignin of hevea brasiliensis. wood is good (Singapore: Springer), 199–209. doi: 10.1007/978-981-10-3115-1_19
Priyadarshan, P. M. (2017). Refinements to hevea rubber breeding. Tree Genet. Genomes. 13 (1), 1–17. doi: 10.1007/s11295-017-1101-8
Qiaoyan, T., Xianqiang, W., Chaofeng, L., Wanxiang, L., Li, Y., Yuanzhong, J., et al. (2013). Functional characterization of the poplar r2r3-myb transcription factor ptomyb216 involved in the regulation of lignin biosynthesis during wood formation. PloS One 8 (10), e76369–. doi: 10.1371/journal.pone.0076369
Quan, M., Du, Q., Xiao, L., Lu, W., Wang, L., Xie, J., et al. (2019). Genetic architecture underlying the lignin biosynthesis pathway involves noncoding RNAs and transcription factors for growth and wood properties in Populus. Plant Biotechnol. J. 17 (1), 302–315. doi: 10.1111/pbi.12978
Quan, M., Xiao, L., Lu, W., Liu, X., Song, F., Si, J., et al. (2018). Association genetics in Populus reveal the allelic interactions of Pto-MIR167a and its targets in wood formation. Front. Plant Sci. 9. doi: 10.3389/fpls.2018.00744
Rhoades, M. W., Reinhart, B. J., Lim, L. P., Burge, C. B., Bartel, B., Bartel, D. P. (2002). Prediction of plant microRNA targets. Cell. 110 (4), 513–520. doi: 10.1016/s0092-8674(02)00863-2
Schmittgen, T. D., Livak, K. J. (2008). Analyzing real-time PCR data by the comparative CT method. Nat. Protoc. 3, 1101–1108. doi: 10.1038/nprot.2008.73
Shen, H., He, X., Poovaiah, C. R., Wuddineh, W. A., Ma, J., Mann, D. G. J., et al. (2012). Functional characterization of the switchgrass (Panicum virgatum) R2R3-MYB transcription factor PvMYB4 for improvement of lignocellulosic feedstocks. New Phytol. 193 (1), 121–136. doi: 10.1111/j.1469-8137.2011.0392
Sibout, R., Eudes, A., Mouille, G., Pollet, B., Lapierre, C., Jouanin, L., et al. (2005). CINNAMYL ALCOHOL DEHYDROGENASE-C and -D are the primary genes involved in lignin biosynthesis in the floral stem of Arabidopsis. Plant Cell. 17, 2059–2076. doi: 10.1105/tpc.105.030767
Song, X. Q., Liu, L. F., Jiang, Y. J., Zhang, B. C., Gao, Y. P., Liu, X. L., et al. (2013). Disruption of secondary wall cellulose biosynthesis alters CADmium translocation and tolerance in rice plants. Mol. Plant 6, 768–780. doi: 10.1093/mp/sst025
Sujan, K. C., Yamamoto, H., Matsuo, M., Yoshida, M., Naito, K., Shirai, T. (2015). Continuum contraction of tension wood fiber induced by repetitive hygrothermal treatment. Wood Sci. Technol. 49, 1157–1169. doi: 10.1007/s00226-015-0762-4
Taylor, N. G., Howells, R. M., Huttly, A. K., Vickers, K., Turner, S. R. (2003). Interactions among three distinct CesA proteins essential for cellulose synthesis. Proc. Natl. Acad. Sci. U. S. A. 100, 1450–1455. doi: 10.1073/pnas.0337628100
Taylor, N. G., Laurie, S., Turner, S. R. (2000). Multiple cellulose synthase catalytic subunits are required for cellulose synthesis in Aarabidopsis. Plant Cell. 12, 2529–2539. doi: 10.1105/tpc.12.12.2529
Taylor, N. G., Scheible, W. R., Cutler, S., Somerville, C. R., Turner, S. R. (1999). The irregular xylem3 locus of Aarabidopsis encodes a cellulose synthase required for secondary cell wall synthesis. Plant Cell. 11, 769–779. doi: 10.1105/tpc.11.5.769
Tempel, S. (2012). Using and understanding RepeatMasker. Methods Mol. Biol. 859, 29–51. doi: 10.1007/978-1-61779-603-6_2
Thaochan, N., Pornsuriya, C., Chairin, T., Sunpapao, A. (2020). Roles of systemic fungicide in antifungal activity and induced defense responses in rubber tree (Hevea brasiliensis) against leaf fall disease caused by Neopestalotiopsis cubana. Physiol. Mol. Plant P 111, 101511. doi: 10.1016/j.pmpp.2020.101511
Vanholme, R., Ralph, J., Akiyama, T., Lu, F., Pazo, J. R., Kim, H., et al. (2010). Engineering traditional monolignols out of lignin by concomitant up-regulation of F5H1 and down-regulation of COMT in Arabidopsis. Plant J. 64 (6), 885–897. doi: 10.1111/j.1365-313X.2010.04353.x
Voelker, S. L., Lachenbruch, B., Meinzer, F. C., Jourdes, M., Ki, C., Patten, A. M., et al. (2010). Antisense down-regulation of 4CL expression alterslignification, tree growth, and saccharification potential of field-grownpoplar. Plant Physiol. 154 (2), 874–886. doi: 10.1104/pp.110.159269
Wagner, A., Donaldson, L., Kim, H., Phillips, L., Flint, H., Steward, D., et al. (2009). Suppression of 4-coumarate-CoA ligase in the coniferous gymnosperm Pinus radiata. Plant Physiol. 149 (1), 370–383. doi: 10.1104/pp.108.125765
Wen, M., Shen, Y., Shi, S., Tang, T. (2012). MiREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinforma. 13, 140. doi: 10.1186/1471-2105-13-140
Weng, J. K., Mo, H., Chapple, C. (2010). Over-expression of F5H in COMT-deficient Arabidopsis leads to enrichment of an unusual lignin and disruption of pollen wall formation. Plant J. 64 (6), 898–911. doi: 10.1111/j.1365-313X.2010.04391.x
Yang, W., Schuster, C., Beahan, C. T., Charoensawan, V., Peaucelle, A., Bacic, A., et al. (2016). Regulation of meristem morphogenesis by cell wall synthases in Arabidopsis. Curr. Biol. 26 (11), 1404–1415. doi: 10.1016/j.cub.2016.04.026
Yu, X., Gong, H., Cao, L., Hou, Y., Qu, S. (2020). MicroRNA397b negatively regulates resistance of Malus hupehensis to Botryosphaeria dothidea by modulating MhLAC7 involved in lignin biosynthesis. Plant Sci. 292, 110390. doi: 10.1016/j.plantsci.2019.110390
Zeng, C., Wang, W., Zheng, Y., Chen, X., Bo, W., Song, S., et al. (2009). Conservation and divergence of microRNAs and their functions in euphorbiaceous plants. Nucleic Acids Res. 38, 981–995. doi: 10.1093/nar/gkp1035
Zhang, B., Deng, L., Qian, Q., Xiong, G., Zeng, D., Li, R., et al. (2009). A missense mutation in the transmembrane domain of CESA4 affects protein abundance in the plasma membrane and results in abnormal cell wall biosynthesis in rice. Plant Mol. Biol. 71, 509–524. doi: 10.1007/s11103-009-9536-4
Zhang, X., Wang, W., Zhu, W., Dong, J., Cheng, Y., Yin, Z., Beahan, C., et al. (2019). Mechanisms and functions of long non-coding RNAs at multiple regulatory levels. Int. J. Mol. Sci. 20 (22), 5573. doi: 10.3390/ijms20225573
Zhang, H., Ying, Y. Q., Wang, J., Zhao, X. H., Zeng, W., Beahan, C., et al. (2018). Transcriptome analysis provides insights into xylogenesis formation in moso bamboo (phyllostachys edulis) shoot. Sci. Rep. 8 (1), 3951. doi: 10.1038/s41598-018-21766-3
Zhong, R., Ye, Z. H. (2009). Transcriptional regulation of lignin biosynthesis. Plant Signal. Behav. 4, 1028–1034. doi: 10.4161/psb.4.11.9875
Zhong, R., Lee, C., Ye, Z. H. (2010). Functional characterization of poplar wood-associated NAC domain transcription factors. Plant Physiol. 152, 1044–1055. doi: 10.1104/pp.109.148270
Keywords: reaction wood, phenylpropanoid biosynthesis pathway, lignin biosynthesis, Hevea brasiliensis, miRNA
Citation: Chen J, Liu M, Meng X, Zhang Y, Wang Y, Jiao N and Chen J (2023) Multiomics studies with co-transformation reveal microRNAs via miRNA-TF-mRNA network participating in wood formation in Hevea brasiliensis. Front. Plant Sci. 14:1068796. doi: 10.3389/fpls.2023.1068796
Received: 13 October 2022; Accepted: 17 May 2023;
Published: 14 August 2023.
Edited by:
Kai-Hua Jia, Shandong Academy of Agricultural Sciences, ChinaCopyright © 2023 Chen, Liu, Meng, Zhang, Wang, Jiao and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jinhui Chen, amluaHVpY2hlbkBoYWluYW51LmVkdS5jbg==