Gang Nie
Zongchao Liao†
Xia Wang
Xinquan Zhang- Department of Forage Science, College of Grassland Science and Technology, Sichuan Agricultural University, Chengdu, China
Chromium (Cr) is a heavy metal in nature, which poses a potential risk to toxicity to both animals and plants when releasing into the environment. However, the regulation of microRNA (miRNA)-mediated response to heavy metal Cr has not been studied in Miscanthus sinensis. In this study, based on high-throughput miRNA sequencing, a total of 104 conserved miRNAs and 158 nonconserved miRNAs were identified. Among them, there were 45 differentially expressed miRNAs in roots and 13 differentially expressed miRNAs in leaves. The hierarchical clustering analysis showed that these miRNAs were preferentially expressed in a certain tissue. There were 833 differentially expressed target genes of 45 miRNAs in roots and 280 differentially expressed target genes of 13 miRNA in leaves. After expression trend analysis, five significantly enriched modules were obtained in roots, and three significantly enriched trend blocks in leaves. Based on the candidate gene annotation and gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) function analysis, miR167a, novel_miR15, and novel_miR22 and their targets were potentially involved in Cr transportation and chelation. Besides, miR156a, miR164, miR396d, and novel_miR155 were identified as participating in the physiological and biochemical metabolisms and the detoxification of Cr of plants. The results demonstrated the critical role of miRNA-mediated responses to Cr treatment in M. sinensis, which involves ion uptake, transport, accumulation, and tolerance characteristics.
Introduction
Heavy metal pollution of soil could lead to losses in agricultural yields and potentially affects human health when it enters the food chain, which has become one of the major ecological problems worldwide (Islam et al., 2014). Chromium (Cr) is a heavy metal by nature with trivalent Cr (III) and hexavalent Cr (VI) as two major stable chemical forms, which poses a potential risk to toxicity to both animals and plants when releasing into the environment (Tiwari et al., 2013). Cr (VI) compounds are known to have stronger oxidizing activity than Cr (III), and excessive accumulation in plants can cause extremely deleterious effects on morphological, physiological, and biochemical processes of plants, such as inhibition of seed germination and plant growth, damage of cell ultrastructure, induction of oxidative stress, and uptake of mineral nutrition (Rout et al., 2000; Dube et al., 2003; Choudhury and Panda, 2005; Shanker et al., 2005; Panda, 2007; Diwan et al., 2010). Previous studies found that Glycine max had a significant translocation of Cr to the aboveground plant parts (Mei et al., 2002). Besides, Spirodela polyrhiza, Hydrodictyon reticulatum, and Ceratophyllum demersum could effectively reduce about 90% of the Cr concentration of water (Rai et al., 1995). Goupil et al. (2009) detected that heat shock proteins (HSPs), metallothioneins (MTs), and GR isoforms were upregulated in tomato (Lycopersicon esculentum) when subjected to Cr stress. In Cr-stressed Typha angustifolia, a significant induction in the expression of ATP synthase, Rubisco small subunit, and coproporphyrinogen III oxidase (Bah et al., 2010) was detected. In short, to manage Cr contamination and limit Cr accumulation in crops, it is necessary to gain a better understanding of the uptake, transportation, and sequestration of Cr and of the adaptive response of plants to the stress.
Miscanthus sinensis is a vigorous C4 perennial grass that originates from East Asia, which is widely cultivated as forage, ornamental grass, and energy crop worldwide (Clifton-Brown and Lewandowski, 2002; Heaton et al., 2004). Miscanthus sinensis could develop a dense net of fibrous roots that contribute to broad adaptation, high water-use efficiency, and great resistance to drought, salt, and heavy metal stress (Ye et al., 2005; Ezaki et al., 2008; Nie et al., 2017). Besides, it has been reported that M. sinensis could effectively absorb and accumulate the Cr (Zhu et al., 2010; He et al., 2013; Liang et al., 2014), and the contaminated biomass could be burnt for bioenergy production with great commercial values (Arduini et al., 2006). On the other hand, about 99% of the Cr taken up by the Miscanthus plant was retained by the hypogeal part (Arduini et al., 2006). Thirty-six proteins were identified to be differently expressed in M. sinensis in response to Cr stress, including oxidative stress-related proteins and metabolism-related proteins (Sharmin et al., 2012). However, the molecular regulation network of Cr detoxification and tolerance characteristics of M. sinensis remain unknown due to the dearth of omic information.
Generally, plant responses to metal toxicity are regulated at both transcriptional and posttranscriptional levels (Mendoza-Soto et al., 2012). MicroRNAs (miRNAs) are a class of noncoding small regulatory RNAs with the length of 20–22 nt, which bind to mRNA and participate in regulating gene expression via histone modification and DNA methylation (Diederichs and Sven, 2011; Khraiwesh et al., 2012; Tiwari et al., 2013), and they play a critical role in the whole growth period or fall into an environmental stress in recent studies (Noman and Aqeel, 2017). In rice (Oryza sativa), several miRNAs, such as miR-156, miR-169, miR-171, miR-199, miR-396, miR-398, miR-399, and miR-408, had been reported as responsive miRNAs under cadmium and arsenic (Tang et al., 2014; Sharma et al., 2015). A high-throughput small RNA sequencing approach revealed that miR159, miR160, miR319, miR396, and miR390 were downregulated in response to Al exposure in Medicago truncatula (Chen et al., 2012). In addition, M. truncatula Hg-responsive miRNAs, such as miR167, miR172, miR169, miR164, and miR395 families, were identified (Zhou et al., 2011). However, few studies on Cr-responsive miRNAs have been reported in recent years. In rice, 13 conserved miRNAs were expressed preferentially in response to Cr stress, namely, miR156, miR159, miR160, miR166, miR169, miR171, miR396, miR397, miR408, miR444, miR1883, miR2877, and miR5072 (Dubey et al., 2020). Fifty-four differentially expressed known miRNAs and 16 novel miRNAs (Raphanus sativus L.), such as miR156, miR164, miR167, miR171, miR396, and miR399, were identified in radish (Liu et al., 2015). In tobacco (Nicotiana tabacum), 53 known miRNAs and 29 unknown miRNAs, including miR156, miR159, miR166, miR167, miR171, miR396, and miR399, were identified (Bukhari et al., 2015).
Given the fact that miRNAs are key regulators, they serve as the core of gene regulatory networks response to heavy metal stress. Therefore, in this study, Illumina HiSeqTM2500 platform was employed for the miRNA sequencing on M. sinensis plant subjected to Cr treatment. The aim of this study is to identify Cr-responsive miRNA and primary potential targets and explore the possible pathway that is involved in tolerance to Cr stress. The results established in this study would help to further elucidate the molecular mechanisms of M. sinensis in response to Cr stress and provide valuable bioenergy resources for the development of soil phytoremediation strategies.
Materials and Methods
Plant Materials
The M. sinensis plant material (M20100819) was provided by Sichuan Agricultural University and was collected from Ya’an, Sichuan (N 29°57'30.1'', E 102°24'36.4', 1,388 m). The seedlings were planted in a round plastic pot (14.3 cm × 11 cm × 11 cm) with a matrix of quartz sand and were irrigated by a nutrient solution of 1× Hoagland under a constant temperature (28°C for light/25°C for darkness with cycles of 12/12 h) and humidity (75%) in a growth chamber. According to the preliminary experiment, the morphology of M. sinensis was completely established after 4 months (Nie et al., 2017), and the phenotypic change was most obvious under 200 mg/L of Cr stress based on the study by Arduini et al. (2006). Therefore, after 4 months of growth, 15 pots were treated with 200 mg/L Cr (VI; K2Cr2O7) solution. After Cr treatment, the samples of roots and leaves (signed for MR and ML) were collected at 0, 12, 24, and 72 h for miRNA sequencing (MR0 and ML0, MR12 and ML12, MR24 and ML24, and MR72 and ML72, respectively), and all of these samples were quickly deep-frozen in a container with liquid nitrogen and then stored at −80°C for later RNA extraction. Each sample has three replicates.
RNA Extraction and High-Throughput Sequencing
Total RNA was extracted from samples using the Direct-zol™ RNA MiniPrep Kit (Zymo Research Co., CA, United States), following the instructions of the manufacturer. Twenty-four RNA samples (MR0 and ML0: roots and leaves at 0 h, nontreated; MR12 and ML12: roots and leaves with Cr treated for 12 h; MR24 and ML24: roots and leaves with Cr treated for 24 h; and MR72 and ML72: roots and leaves with Cr treated for 72 h) were sent to the Tianjin Novogene Bioinformatic Technology Co., Ltd., (TianJing, China) for sequencing using the Illumina HiSeqTM2500.
Establishment of Small RNA Library and Bioinformatic Identification of miRNAs
An equal amount of RNA (RIN number > 7.0) per sample was used to construct small RNA libraries using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, United States), following the recommendations of the manufacturer. The 3' and 5' end adapters were ligated to the two ends of miRNAs. Then, the first strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H–). Long Amp Taq 2X Master Mix, SR Primer for Illumina, and index (X) primer were used for PCR amplification. The PCR products were further purified on 8% polyacrylamide gel (100 V, 80 min). The DNA fragments corresponding to 140–160 bp were recovered and dissolved in 8-μl elution buffer. Finally, the library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA high-sensitivity chips. Small RNA sequencing was performed on the Illumina HiSeqTM2500 platform to generate raw reads. The clean reads were obtained by removing reads containing poly N, with 5' adapter contaminants, without 3' adapter or the insert tag, containing poly A or T or G or C, and containing low quality and those smaller than 18 nt reads from the raw data. The small RNA tags were mapped to the Sorghum bicolor reference sequence (ftp://ftp.ensemblgenomes.org/pub/plants/release-41/gff3/sorghum bicolor; Langmead et al., 2009) without mismatch to analyze their expression and distribution on the reference. The mapped small RNA tags were used to identify known miRNA by aligning against the miRbase 20.0. The custom scripts were used to obtain the miRNA counts and base bias on the first position of an identified miRNA with certain length and on each position of all identified miRNAs, respectively. The characteristics of the hairpin structure of miRNA precursor can be used to predict novel miRNAs. The miREvo (Wen et al., 2012) and mirdeep2 (Friedlnder et al., 2012) software were integrated to predict novel miRNAs. The known miRNA used miFam.dat1 to identify families, and novel miRNA precursors were further submitted to Rfam2 to identify Rfam families.
Differential Expression Analysis of miRNAs
According to the methods studied by Zhou et al. (2010), the miRNA expression levels were estimated by transcript per million (TPM). TPM avoided the effect of quantitative accuracy and normalized the expression level of small RNAs in different sequencing amounts. TPM was calculated as follows:
(libsize: Sum of miRNA ReadCount of sample)
The data from the TPM normalization were used to identify the differentially expressed miRNAs based on the number of genomic tags in each sample and compare the miRNA abundance among the three sets of libraries. The fold change was calculated as follows:
A: a sample of A time.
B: a sample of B time.
The positive and negative values indicated upregulation and downregulation of small RNA, respectively. The Benjamini–Hochberg method was used to adjust the p value. Value of p < 0.05 was set as the threshold for significantly differential expression. To identify differentially expressed miRNAs among Cr treatment at different points in time, heat maps were generated using Novomagic tools.
Target Gene Prediction, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes Function Analysis
The target gene prediction of miRNAs was performed by psRobot_tar in psRobot (Wu et al., 2012). NovoMagic program3 was used to draw the Venn diagram. After finding the target genes of miRNAs, the same program was used for the GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The target genes were divided into three parts using the GO function package of R-Studio Team as follows: biological processes, molecular functions, and cellular components. The trend analysis was performed using the OmicShare tools.4 The small RNA sequencing data set was deposited into the NCBI database under accession number PRJNA595773.
Verification of miRNA by qRT-PCR Analysis
Six Cr-responsive miRNAs, i.e., miR156a, miR167a, miR396d, miR5564a, novel_miR149, and novel_miR32, were selected for quantitative real time PCR (qRT-PCR) analysis (Table 1). The miRNA tailing reverse transcription primer was used for the transcription of total RNA, which was designed according to the method of Shi and Chiang (Rui and Chiang, 2005). Reverse transcription was performed using the Mir-X miRNA First-Strand Synthesis Kit (Clontech Laboratories, Inc., Shanghai, China) according to the instructions of the manufacturer. The cDNA (i.e., samples were same with sequence) was then used for qRT-PCR using specific forward primers (Table 1) and mRQ 3' universal reverse primers (10 μM). The 18S rRNA was used as an internal reference for tailing qRT-PCR. The qRT-PCR was performed using a Bio-Rad CFX96 Touch (Bio-Rad, Hercules, CA, United States). Each 20 μl reaction contained 10 μl 2× miRNA qPCR master mix, 0.8 μl (10 μM) forward primer, 0.8 μl (10 μM) reverse primer, 6.4 μl RNase-free water, and 2 μl cDNA template. The PCR procedure was carried out at 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and at 60°C for 20 s. After the reaction was completed, a melting curve was set to evaluate primer specificity. The relative expression levels were calculated using the 2−△△CT method (Livak and Schmittgen, 2001).
Table 1. Primers information of qRT-PCR verification for selected six Miscanthus sinensis miRNAs.
Results
Data Analysis of Small RNA Sequences
After sequencing with previously cultured M. sinensis materials that were grown under 0 h Cr treatment (ML0 and MR0), 12 h Cr treatment (ML12 and MR12), 24 h Cr treatment (ML24 and MR24), and 72 h Cr treatment (ML72 and MR72), a total of 24 small RNA libraries were established to identify and predict the miRNAs associated with Cr treatment. The sequences of shortages, low quality tags, N > 10% reads, adaptors, and adaptor–adaptor ligation were filtered out, and 355,453,701 clean reads were obtained (Supplementary Table S1). About 0.52–1.04 Gb raw reads were obtained from small RNA sequencing, and the averages of Q20 and Q30 were more than 99.55 and 98.74%, respectively. The miRNA length ranging from 18 to 30 nt was further mapped to the reference sequence using Bowtie software, with the average mapping rate of 62.20%.
Identification of Cr-Responsive Differentially Expressed miRNAs in M. sinensis
A total of 104 known miRNAs (Supplementary Table S2) were identified in 24 libraries, and 158 novel miRNAs (Supplementary Table S3) were predicted, which belong to 35 miRNA families. Among them, 54 differentially expressed miRNAs were detected (in comparison group: ML12 vs. ML0, ML24 vs. ML12, ML72 vs. ML24, MR12 vs. MR0, MR24 vs. MR12, and MR72 vs. MR24), including nine known miRNAs and 45 novel miRNAs, and four miRNAs (i.e., novel_miR149, novel_miR15, novel_miR22, and novel_miR32) are commonly identified in two plant organizations (Table 2). Out of 45 novel miRNAs, 13 were identified in the three comparison groups of leaves. In the root, there were 45 differentially expressed miRNAs, i.e., nine known miRNAs and 36 novel miRNAs. After Cr stress, 31 miRNAs were upregulated, 38 miRNAs were downregulated, and the most dramatic change occurred at the roots (Figure 1A). The hierarchical clustering analysis showed that these miRNAs were preferentially expressed in a certain tissue, most of which were highly expressed in the root tissue (Figure 1B). Especially, the expression trend of novel_miR15 and novel_miR22 is consistent in roots and leaves, and these two miRNAs could cluster together in a heat map. In addition, novel_miR15 and novel_miR22 clustered with the miR164 family. According to the expression of known miRNA, the miRNAs from the same family had similar expression pattern, for example, both miRNAs from miR396 family and miR164 family were upregulated in the root after the Cr stress. Furthermore, the Venn diagrams were used to compare and analyze the miRNAs differentially expressed at six comparison groups (Figure 1C). Among them, no common miRNA was found in these three groups of leaves, while there were 11 differentially expressed miRNAs (i.e., novel_miR14, novel_miR118, novel_miR32, novel_miR15, miR156a, miR167a, novel_miR182, novel_miR140, novel_miR183, novel_miR33, and novel_miR8) in the roots shared by MR12 vs. MR0 and MR72 vs. MR24. Interestingly, only one differentially expressed miRNA was identified in MR24 vs. MR12 comparable group (novel_miR155), indicating rare response at posttranscriptional levels to Cr stress during this stage.
Table 2. Numbers of differentially expressed miRNA targets in leaf and root organizations.
Figure 1. The identified differentially expressed microRNA (miRNA) response to Chromium (Cr) stress. (A) The histogram showed the number of upregulated and downregulated differentially expressed miRNAs caused by Cr treatment of six comparison groups. (B) Heat map of the expression levels of 54 differentially expressed miRNAs. Different colors indicate different expression levels, as shown in the scale. (C) Venn diagram showed the number of differentially expressed miRNAs in different comparison groups.
Prediction of Target Genes of the Differentially Expressed miRNAs
A total of 5,451 target genes of 54 differentially expressed miRNAs were predicted, and each miRNA was annotated to 123 genes on average (Table 2). The most abundant novel_miR118 corresponded to 222 target genes. However, novel_miR14 and novel_miR194 were not annotated to any target genes. Among them, 3,215 target genes were also differentially expressed in the transcriptome sequencing analysis after Cr treatment (|Fold change| > 1, FDR < 0.05, the data set was deposited into the NCBI database with accession number PRJNA594536, but not yet published).
A total of 833 differentially expressed target genes of 45 miRNAs were identified in roots. After trend analysis, five significantly enriched modules were obtained showing various expression trends; including continuous up- and downregulating in roots (Figure 2A; Supplementary Table S4). The number of genes in profile 10 was 188, which were the target genes of 37 miRNAs. Profile 10 enriched the most differentially expressed target genes in roots under Cr stress. Among them, there is a sulfate transmembrane transporter (Cluster-61280.0), which was targeted by miR167a. Besides, genes in profile 17 were 74, which were the target genes of 40 miRNAs, and a Mg2+ transporter protein (i.e., Cluster-42176.0, targeted by novel_miR22) and two ABC transporter proteins (i.e., Cluster-86748.1466 was the target gene of novel_miR15, and Cluster-25283.0 was the target gene of novel_miR3 and novel_miR8) were detected in this study. Interestingly, a M. sinensis MT-like protein gene (i.e., Cluster-86748.41948 in profile 6) was also the target gene of novel_miR15.
Figure 2. Trend analysis of the differentially expressed target genes of specific miRNAs in (A) roots and (B) leaves, and profiles are ordered based on the number of genes assigned.
In addition, the genes in profile 19 were 112, which were the target genes of 41 miRNAs, and genes in profile 19 showed a continuously upregulated expression trend. An activated mitogen-activated protein kinase (Cluster-15974.0) was targeted by miR396d, a probable leucine-rich repeat receptor-like serine/threonine-protein kinase (Cluster-86748.44739) was targeted by miR396e, and a pyruvate dehydrogenase kinase (Cluster-86748.49811) was targeted by novel_miR44, which can participate in the signal transduction in plants after Cr stress through phosphorylation. In profile 0, the number of genes was 80, which were the target genes of 32 miRNAs. A cytochrome P450 gene (Cluster-86748.26078 was the target gene of miR156a) and an auxin-responsive protein (Cluster-86748.61059, target by novel_miR124) were detected, which potentially participated in the plant oxidation–reduction process and showed continuous downregulating expression trend in the block. There were 71 target genes in profile 15, which were the target genes of 32 miRNAs, and a glutaminyl-peptide cyclotransferase (i.e., Cluster-56602.0 was the target gene of miR156a) was detected in this study.
In leaves, the number of differentially expressed target genes of 13 miRNAs was 280, and three significantly enriched trend blocks were obtained after the trend analysis (Figure 2B; Supplementary Table S5). The number of genes in profile 14 was 51, including a MT protein gene (i.e., Cluster-86748.38482 was targeted by novel_miR155) participating in metal ion binding, one transcription factor of MYB (i.e., Cluster-86748.56774 was targeted by novel_miR17), and one transcription factor of bHLH129 isoform X2 (i.e., Cluster-86748.36885 was targeted by novel_miR149).
The number of genes in profile 7 was 37, which were the target genes of 11 differentially expressed miRNAs, and the number of genes in profile 5 was 30, which were the target genes of 12 miRNAs. In profile 7, three chloroplast OS genes (i.e., Cluster-86748.59528 was targeted by novel_miR1, and Cluster-86748.60212 and Cluster-86748.60213 were targeted by novel_miR155) and a structural maintenance of chromosome protein (Cluster-97477.0 was targeted by novel_miR1) were detected. In profile 5, Cluster-86748.43155 targeted by novel_miR13 plays a role of protein kinase activity at cell membrane.
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Function Analysis of Differentially Expressed Target Genes
In GO function analysis, two different organs of M. sinensis showed different enrichment trends. The function of the target genes of differentially expressed miRNAs in root was annotated with three categories based on their functions, namely, biological processes, cellular components, and molecular functions; among them, biological processes were mainly annotated (Figure 3). For cellular components, these target genes were annotated with cytoplasmic part, intracellular organelle, and organelle. In molecular functions, the gene function is mainly concentrated in oxidoreductase activity, structural constituent of ribosome and transferase activity, and the transferring glycosyl group. In GO function analysis of leaves, molecular functions play a major role in three function analysis, namely, ion binding, ATP binding, and phosphotransferase activity. However, biological processes and cellular components enriched fewer genes. For KEGG analysis, unigenes were searched through the KEGG database to determine the biological processes in which they are involved. Roots were significantly enriched in starch and sucrose metabolism, and leaves were enriched in the ribosome and plant hormone signal transduction (Figure 4).
Figure 3. Gene ontology (GO) diagrams of target genes of differentially expressed miRNAs in (A) roots and (B) leaves. Green, blue, and red colors represent biological processes, cellular components, and molecular functions, respectively.
Figure 4. Kyoto Encyclopedia of Genes and Genomes (KEGG) diagrams of target genes of differentially expressed miRNAs in (A) roots and (B) leaves.
qRT-PCR Analysis of Differentially Expressed miRNA
In order to verify the accuracy of high-throughput sequencing results and validate the expression pattern of miRNAs detected in this study, six differentially expressed miRNAs were randomly selected for qRT-PCR validation. The results showed that the relative expression obtained by qRT-PCR was basically consistent with the TPM values, and most of them have specific expression in a certain tissue (Figure 5). However, there are still some differences between transcriptomic sequencing and qRT-PCR. The expression of miR156a increased continuously in TPM expression, but the qRT-PCR expression of miR156a decreased sharply in MR72, which might be due to the sensitivity and specificity of the two techniques.
Figure 5. Quantitative real time PCR (qRT-PCR) validation of differentially expressed miRNAs detected by miRNA sequencing. Vertical bars indicate mean values ± SD from three biological replicates.
Discussion
Due to the rapid development of industry and the irregular discharge of industrial wastes, heavy metal pollution has become one of the most urgent problems in ecological remediation (He et al., 2005). Previous studies indicated that miRNAs are involved in different heavy metal stress, participating in regulation of gene expression. It has been reported that Miscanthus is resistant to heavy metals and can remove heavy metal soil activity (Zhang et al., 2009). However, there are few studies on the extensive identification of Cr-reactive miRNA and its target genes in M. sinensis. In this study, 104 conservative miRNAs and 158 nonconservative miRNAs were identified from root and leaf library. Among them, compared with nonconservative miRNAs, most conservative miRNAs had relatively high reads. In addition, the mean number of members of the conservative miRNA family is smaller than that of the nonconservative miRNA family, which is different from previous studies in other species, such as Fragaria × ananassa (Li et al., 2013), M. truncatula (Zhou et al., 2011), and Brassica napus (Chen et al., 2012). In underlying radish plants, 54 conservative miRNAs and 16 novel differentially expressed miRNAs were identified under Cr stress (Liu et al., 2015). In tobacco (N. tabacum) roots, 29 conserved miRNAs and 14 novel miRNAs expressed differently under Cr stress (Bukhari et al., 2015). A total of 512 and 568 known miRNAs were identified from Cr treatment and control rice plants, respectively, and 13 conserved miRNAs depicted preferential expression under Cr treatment (Dubey et al., 2020). In total, nine known and 45 novel miRNAs were significantly differentially regulated under Cr stress in this study. Furthermore, there were no known miRNA expression differences in the leaves, and all of them depicted differentially regulated in roots.
Plant response to heavy metal stress was a complex process, including dynamic changes of a large number of miRNAs and target genes. Novel_miR15 and novel_miR22 were two differentially expressed miRNAs in both roots and leaves in this study, which were clustered with miR164a and miR164b, showing an increasing tendency after the treatment. Previous studies indicated that miRNA164 was widely involved not only in plant seed germination, cell division, cell secondary wall formation, meristems formation, flowering, and senescence (Wang et al., 2014), but also in plant response to abiotic stress, such as heavy metal stress (Huang et al., 2009; Zhou et al., 2011). The target genes of miRNA164 family were NAC transcription factors (Xu et al., 2016), and it was found that overexpression of miR164a and miR164b resulted in the reduced expression of NAC gene family members, CUC1 and CUC2, in Arabidopsis, which were associated with the development of meristematic tissue and the establishment of organ primordium boundaries (Mallory et al., 2004; Nikovics et al., 2006). Accordingly, it was assumed that Cr increased the expression level of these miRNAs, inhibiting the expression of NAC genes related to plant growth and development; thus, the distribution of nutrient and energy could be a turn to plant resistance.
Sulfate transporter (ST) protein had been identified as a major transporter for Cr (VI) uptake. In yeast, the STs (i.e., Sul1 and Sul2) are responsible for Cr (VI) uptake (Chang et al., 2003). Furthermore, a recent study showed that high-affinity STs (i.e., Sultr1 and 2) were involved in Cr (VI) uptake by the roots of Arabidopsis thaliana (Xu et al., 2021). An increasing number of studies have demonstrated that miR395 plays an important role in the regulation of STs, and it is crucial for the regulation of sulfate homeostasis (He et al., 2020). In this study, a target gene “Cluster-61280.0” (targeted by miR167a) was identified to regulate the synthesis of Sultr2 transporter. Therefore, miR167a may be a new miRNA involved in the regulation of sulfate proteins under Cr stress. MT protein has metal binding ability and excessive metal ion detoxification (García-Hernández and Taiz, 1998; Jain et al., 2016), and it has high affinity for various heavy metals due to its rich cysteine content (Cobbett and Goldsbrough, 2002). Jain et al. (2016) found that MT gene expression was significantly increased in sugarcane (Saccharum spp. hybrid) under Cr stress. In the study involving the root of M. sinensis, a MT-like protein gene (Cluster-86748.41948) was identified as an important predicted target gene of novel_miR15, which occupied an important position in ion transport function and participated in the detoxification of heavy metal Cr. Besides; novel_miR15 also targeted an ABC transporter protein (Cluster-86748.1466), which was differentially expressed after Cr treatment. The ABC transporter family is one of the largest membrane protein families with a wide variety of species distributed throughout the organism and a variety of transportable substrates, including ions, carbohydrates, lipids, and heavy metals (i.e., As, Cd, As, Hg, and Zn; Long et al., 2011; Wang et al., 2017). Accordingly, it was speculated that novel_miR15 and miR164 had a similar heavy metal detoxification effect based on the clustered results in heat map and predicted target genes of them.
In M. sinensis, most differentially expressed miRNAs were identified in roots after Cr treatment, including known miR156 and miR167. For miR156, it was generally identified as a critical factor to help plants to coordinate the relationship between development and stress tolerance and targeted Squamosa promoter binding protein-like genes family. Overexpression of miR156 in A. thaliana could increase plant stress tolerance (Cui et al., 2015). In addition, a target of miR156 in B. napus encodes glutamyl transferase, an enzyme involved in the formation of chelate states in plant, through which the plant formed chelate states or glutathione tripeptides chelated the heavy metals and eventually sequestered them into vacuoles (Cobbett, 2000; Zhou et al., 2011). Besides, the highly toxic hexavalent Cr ions are redox to the less toxic trivalent Cr ions, which are then transported in a possible detoxification pathway. Cytochrome P450 proteins are important in the oxidative, peroxidative, and reductive metabolisms of numerous endogenous compounds (Joo, 1997). In rice, CYP90D2 encodes C-3 dehydrogenase and has been considered to be involved in brassinosteroid biosynthesis and catabolism (Hong et al., 2003). In this study, cytochrome P450 (Cluster-86748.26078) was also identified as an important target gene of miR156 under Cr stress, which would provide a new insight for further explanation of the detoxification mechanism. MiR167 targeting auxin response factor (ARF) gene family also play an important role in regulating root extension and plant shape control (Goetz et al., 2006; Gutierrez et al., 2009). Under salt stress, miR167 was upregulated in both A. thaliana (Liu et al., 2008) and Caragana intermedia (Zhu et al., 2013), possibly by regulating the expression of ARFs, affecting the transcription of auxin response-related genes, thus slowing down the growth and development of plants.
In conclusion, a total of 104 conserved miRNAs and 158 nonconserved miRNAs were identified. Among them, 45 miRNAs were differentially expressed miRNAs in roots and 13 miRNAs were differentially expressed miRNAs in leaves. There were 833 differentially expressed target genes of 45 miRNAs in roots and 280 differentially expressed target genes in leaves. Based on candidate gene annotation and GO and KEGG function analysis, the miR167a, novel_miR15, and novel_miR22 and their targets were potentially involved in Cr transportation and chelation. Besides, miR156a, miR164, miR396d, and novel_miR155 were identified by participating in the physiological and biochemical metabolisms and the detoxification of Cr of plants. These results provided a valuable reference original insight into the role of miRNAs on the metabolic fingerprinting of Cr, regarding the targeted miRNAs. In future, the identified miRNAs and target genes will provide a fundamental process toward understanding miRNA-mediated regulatory networks of M. sinensis response under Cr (VI) stress.
Data Availability Statement
The original contributions presented in the study are publicly available. This data can be found here: miRNA-seq data were deposited into NCBI database under accession number PRJNA595773, and the transcriptome sequence data were deposited into NCBI database under accession number PRJNA594536.
Author Contributions
XZ and GN contributed to the experiment design and funding acquisition. MZ, ZL, JZ, and AL performed the material preparation and data collection. GN, MZ, XW, and JC analyzed the data. GN, ZL, MZ, and AL wrote the first draft of the manuscript. All authors commented on previous versions of the manuscript, and read and approved the final manuscript.
Funding
This work was supported by the National Project on Sci-Tec Foundation Resources Survey (2017FY100602), China Agriculture Research System of MOF and MARA, the National Natural Science Foundation of China (NSFC 31801432), and the Ministry of Education of the People’s Republic of China “Chunhui” Program Cooperation Project (20191383-101).
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.666117/full#supplementary-material
Footnotes
1. ^http://www.mirbase.org/ftp.shtml
2. ^http://rfam.sanger.ac.uk/search/
References
Arduini, I., Masoni, A., and Ercoli, L. (2006). Effects of high chromium applications on Miscanthus during the period of maximum growth. Environ. Exp. Bot. 58, 234–243. doi: 10.1016/j.envexpbot.2005.09.004
Bah, A. M., Sun, H., Fei, C., Jing, Z., Dai, H., Zhang, G., et al. (2010). Comparative proteomic analysis of Typha angustifolia leaf under chromium, cadmium and lead stress. J. Hazard. Mater. 184, 191–203. doi: 10.1016/j.jhazmat.2010.08.023
Bukhari, S. A. H., Shang, S., Zhang, M., Zheng, W., Zhang, G., Wang, T. Z., et al. (2015). Genome-wide identification of chromium stress-responsive micro RNAs and their target genes in tobacco (Nicotiana tabacum) roots. Environ. Toxicol. Chem. 34, 2573–2582. doi: 10.1002/etc.3097
Chang, K. S., Won, J. I., Lee, M. R., Lee, C. E., and Hwang, S. (2003). The putative transcriptional activator MSN1 promotes chromium accumulation in Saccharomyces cerevisiae. Mol. Cell 16, 291–296.
Chen, L., Wang, T., Zhao, M., Tian, Q., and Zhang, W. H. (2012). Identification of aluminum-responsive microRNAs in Medicago truncatula by genome-wide high-throughput sequencing. Planta 235, 375–386. doi: 10.1007/s00425-011-1514-9
Choudhury, S., and Panda, S. K. (2005). Toxic effects, oxidative stress and ultrastructural changes in Moss Taxithelium Nepalense (Schwaegr.) broth. Under chromium and lead phytotoxicity. Water Air Soil Pollut. 167, 73–90. doi: 10.1007/s11270-005-8682-9
Clifton-Brown, J. C., and Lewandowski, I. (2002). Screening Miscanthus genotypes in field trials to optimise biomass yield and quality in southern Germany. Eur. J. Agron. 16, 97–110. doi: 10.1016/S1161-0301(01)00120-4
Cobbett, C. S. (2000). Phytochelatin biosynthesis and function in heavy-metal detoxification. Curr. Opin. Plant Biol. 3, 211–216. doi: 10.1016/S1369-5266(00)00066-2
Cobbett, C. S., and Goldsbrough, P. (2002). Phytochelatins and metallothioneins: roles in heavy metal detoxification and homeostasis. Annu. Rev. Plant Biol. 53, 159–182. doi: 10.1146/annurev.arplant.53.100301.135154
Cui, L. G., Shan, J. X., Shi, M., Gao, J. P., and Lin, H. X. (2015). The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants. Plant J. 80, 1108–1117. doi: 10.1111/tpj.12712
Diederichs, J. W., and Sven (2011). Argonaute proteins regulate microRNA stability: increased microRNA abundance by Argonaute proteins is due to microRNA stabilization. RNA Biol. 8, 1149–1157. doi: 10.4161/rna.8.6.17665
Diwan, H., Khan, I., Ahmad, A., and Iqbal, M. (2010). Induction of phytochelatins and antioxidant defence system in Brassica juncea and Vigna radiata in response to chromium treatments. Plant Growth Regul. 61, 97–107. doi: 10.1007/s10725-010-9454-0
Dube, B. K., Tewari, K., Chatterjee, J., and Chatterjee, C. (2003). Excess chromium alters uptake and translocation of certain nutrients in citrullus. Chemosphere 53, 1147–1153. doi: 10.1016/S0045-6535(03)00570-8
Dubey, S., Saxena, S., Chauhan, A. S., Mathur, P., and Chakrabaroty, D. (2020). Identification and expression analysis of conserved microRNAs during short and prolonged chromium stress in rice (Oryza sativa). Environ. Sci. Pollut. Res. Int. 27, 380–390. doi: 10.1007/s11356-019-06760-0
Ezaki, B., Nagao, E., Yamamoto, Y., Nakashima, S., and Enomoto, T. (2008). Wild plants, Andropogon virginicus L. and Miscanthus sinensis Anders, are tolerant to multiple stresses including aluminum, heavy metals and oxidative stresses. Plant Cell Rep. 27, 951–961. doi: 10.1007/s00299-007-0503-8
Friedlnder, M. R., Mackowiak, S. D., Na, L., Wei, C., and Nikolaus, R. R. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids 40, 37–52. doi: 10.1093/nar/gkr688
García-Hernández, M., and Taiz, M. L. (1998). Metallothioneins 1 and 2 have distinct but overlapping expression patterns in Arabidopsis. Plant Physiol. 118, 387–397. doi: 10.1104/pp.118.2.387
Goetz, M., Vivian-Smith, A., Johnson, S. D., and Koltunow, A. M. (2006). AUXIN RESPONSE FACTOR8 is a negative regulator of fruit initiation in Arabidopsis. Plant Cell 18, 1873–1886. doi: 10.1105/tpc.105.037192
Goupil, P., Souguir, D., Ferjani, E., Faure, O., Hitmi, A., and Ledoigt, G. (2009). Expression of stress-related genes in tomato plants exposed to arsenic and chromium in nutrient solution. J. Plant Physiol. 166, 1446–1452. doi: 10.1016/j.jplph.2009.01.015
Gutierrez, L., Bussell, J. D., Pacurar, D. I., Schwambach, J., Pacurar, M., and Bellini, C. (2009). Phenotypic plasticity of adventitious rooting in Arabidopsis is controlled by complex regulation of auxin response factor transcripts and microRNA abundance. Plant Cell 21, 3119–3132. doi: 10.1105/tpc.108.064758
He, J., Li, H., Luo, J., Ma, C., Li, S., Qu, L., et al. (2013). A transcriptomic network underlies microstructural and physiological responses to cadmium in Populus x canescens. Plant Physiol. 162, 424–439. doi: 10.1104/pp.113.215681
He, L., Wang, H., Zhao, Q., Cheng, Z., and Liu, W. (2020). Tomato grafting onto Torubamu (Solanum melongena): miR166a and miR395b reduce scion cd accumulation by regulating sulfur transport. Soil Sci. Plant Nutr. 452, 267–279. doi: 10.1007/s11104-020-04564-7
He, Z. L., Yang, X. E., and Stoffella, P. J. (2005). Trace elements in agroecosystems and impacts on the environment. Trace Elem. Med. Biol. 19, 125–140. doi: 10.1016/j.jtemb.2005.02.010
Heaton, E. A., Long, S. P., Voigt, T. B., Jones, M. B., and Clifton-Brown, J. (2004). Miscanthus for renewable energy generation: European Union expierence and projections for Illinois. Mitig. Adapt. Strateg. Glob. Chang. 9, 433–451. doi: 10.1023/B:MITI.0000038848.94134.be
Hong, Z., Ueguchi-Tanaka, M., Umemura, K., Uozu, S., Fujioka, S., Takatsuto, S., et al. (2003). A rice brassinosteroid-deficient mutant, ebisu dwarf (d2), is caused by a loss of function of a new member of cytochrome P450. Plant Cell 15, 2900–2910. doi: 10.1105/tpc.014712
Huang, S. Q., Peng, J., Qiu, C. X., and Yang, Z. M. (2009). Heavy metal-regulated new microRNAs from rice. J. Inorg. Biochem. 103, 282–287. doi: 10.1016/j.jinorgbio.2008.10.019
Islam, M. K., Alam, I., Khanam, S., Lee, S. Y., and Huh, M. R. (2014). Accumulation and tolerance characteristics of chromium in nine jute varieties (Corchorus spp. and Hibiscus spp.). Plant Omics 7, 392–402.
Jain, R., Singh, A., Singh, P., Solomon, S., Tripathi, P., Chandra, A., et al. (2016). Study on physio-biochemical attributes and metallothionein gene expression affected by chromium (VI) in sugarcane (Saccharum spp. hybrid). J. Environ. Biol. 37:375.
Joo, E. Y. (1997). Induction of cytochrome P-450 protein in wounded avocado roots (Persea americana, mill). Prev. Nutr. Food Sci. 2, 328–332.
Khraiwesh, B., Zhu, J. K., and Zhu, J. (2012). Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochim. Biophys. Acta 1819, 137–148. doi: 10.1016/j.bbagrm.2011.05.001
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
Li, H., Mao, W., Liu, W., Dai, H., Liu, Y., Ma, Y., et al. (2013). Deep sequencing discovery of novel and conserved microRNAs in wild type and a white-flesh mutant strawberry. Planta 238, 695–713. doi: 10.1007/s00425-013-1917-x
Liang, X., Liu, J., Shi, S., Yuan, W., Chang, E., Ming, G., et al. (2014). Uptake of heavy metals by native herbaceous plants in an antimony mine (Hunan, China). CLEAN–Soil Air Water 42, 81–87. doi: 10.1002/clen.201200490
Liu, H. H., Tian, X., Li, Y. J., Wu, C. A., and Zheng, C. C. (2008). Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA 14, 836–843. doi: 10.1261/rna.895308
Liu, W., Xu, L., Wang, Y., Shen, H., Zhu, X., Zhang, K., et al. (2015). Transcriptome-wide analysis of chromium-stress responsive microRNAs to explore miRNA-mediated regulatory networks in radish (Raphanus sativus L.). Sci. Rep. 5:14024. doi: 10.1038/srep14024
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Long, Y., Li, Q., and Cui, Z. (2011). Molecular analysis and heavy metal detoxification of ABCC1/MRP1 in zebrafish. Mol. Biol. Rep. 38, 1703–1711. doi: 10.1007/s11033-010-0283-z
Mallory, A. C., Dugas, D. V., Bartel, D. P., and Bartel, B. (2004). MicroRNA regulation of NAC-domain targets is required for proper formation and separation of adjacent embryonic, vegetative, and floral organs. Curr. Biol. 14, 1035–1046. doi: 10.1016/j.cub.2004.06.022
Mei, B., Puryear, J. D., and Newton, R. J. (2002). Assessment of Cr tolerance and accumulation in selected plant species. Plant Soil 247, 223–231. doi: 10.1023/A:1021509115343
Mendoza-Soto, A. B., Sánchez, F., and Hernández, G. (2012). MicroRNAs as regulators in plant metal toxicity response. Front. Plant Sci. 3:105. doi: 10.3389/fpls.2012.00105
Nie, G., Tang, L., Zhang, Y., Huang, L., Ma, X., Cao, X., et al. (2017). Development of SSR markers based on transcriptome sequencing and association analysis with drought tolerance in perennial grass Miscanthus from China. Front. Plant Sci. 8:801. doi: 10.3389/fpls.2017.00801
Nikovics, K., Blein, T., Peaucelle, A., Ishida, T., Morin, H., Aida, M., et al. (2006). The balance between the MIR164A and CUC2 genes controls leaf margin serration in Arabidopsis. Plant Cell 18, 2929–2945. doi: 10.1105/tpc.106.045617
Noman, A., and Aqeel, M. (2017). miRNA-based heavy metal homeostasis and plant growth. Environ. Sci. Pollut. Res. Int. 24, 10068–10082. doi: 10.1007/s11356-017-8593-5
Panda, S. K. (2007). Chromium-mediated oxidative stress and ultrastructural changes in root cells of developing rice seedlings. J. Plant Physiol. 164, 1419–1428. doi: 10.1016/j.jplph.2007.01.012
Rai, U. N., Sinha, S., Tripathi, R. D., and Chandra, P. (1995). Wastewater treatability potential of some aquatic macrophytes: removal of heavy metals. Ecol. Eng. 5, 5–12. doi: 10.1016/0925-8574(95)00011-7
Rout, G. R., Samantaray, S., and Das, P. (2000). Effects of chromium and nickel on germination and growth in tolerant and non-tolerant populations of Echinochloa colona (L.) link. Chemosphere 40, 855–859. doi: 10.1016/S0045-6535(99)00303-3
Rui, R. S., and Chiang, V. L. V. (2005). Facile means for quantifying microRNA expression by real-time PCR. BioTechniques 39, 519–525. doi: 10.2144/000112010
Shanker, A. K., Cervantes, C., Loza-Tavera, H., and Avudainayagam, S. (2005). Chromium toxicity in plants. Environ. Int. 31, 739–753. doi: 10.1016/j.envint.2005.02.003
Sharma, D., Tiwari, M., Lakhwani, D., Tripathi, R. D., and Trivedi, P. K. (2015). Differential expression of microRNAs by arsenate and arsenite stress in natural accessions of rice. Metallomics 7, 174–187. doi: 10.1039/C4MT00264D
Sharmin, S. A., Alam, I., Kim, K. H., Kim, Y. G., Kim, P. J., Bahk, J. D., et al. (2012). Chromium-induced physiological and proteomic alterations in roots of Miscanthus sinensis. Plant Sci. 187, 113–126. doi: 10.1016/j.plantsci.2012.02.002
Tang, M., Mao, D., Xu, L., Li, D., Song, S., and Chen, C. (2014). Integrated analysis of miRNA and mRNA expression profiles in response to cd exposure in rice seedlings. BMC Genomics 15:835. doi: 10.1186/1471-2164-15-835
Tiwari, K. K., Singh, N. K., and Rai, U. N. (2013). Chromium Phytotoxicity in radish (Raphanus sativus): effects on metabolism and nutrient uptake. Bull. Environ. Contam. Toxicol. 91, 339–344. doi: 10.1007/s00128-013-1047-y
Wang, R. F., Hu, Y. S., Gao, W. R., Zhang, Y. X., and Song, X. S. (2014). Functions of NAC transcription factors family in stress response in plants. Plant Physiol. J. 50, 1494–1500. doi: 10.13592/j.cnki.ppj.2014.0363
Wang, X. Z., Sun, W. M., Yi Feng, M. A., Han, E. Q., Han, L., Sun, L. P., et al. (2017). Research progress of ABC transporters in Arabidopsis thaliana. Plant Physiol. J. 02, 4–15. doi: 10.13592/j.cnki.ppj.2016.0404
Wen, M., Shen, Y., Shi, S., and Tang, T. (2012). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140. doi: 10.1186/1471-2105-13-140
Wu, H. J., Ma, Y. K., Tong, C., Meng, W., and Wang, X. J. R. (2012). PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids 40, W22–W28. doi: 10.1093/nar/gks554
Xu, Z. R., Cai, M. L., Chen, S. H., Huang, X. Y., and Wang, P. (2021). High-affinity sulfate transporter Sultr1;2 is a major transporter for Cr(VI) uptake in plants. Environ. Sci. Technol. 55, 1576–1584. doi: 10.1021/acs.est.0c04384
Xu, S., Liu, N., Mao, W., Hu, Q., Wang, G., and Gong, Y. (2016). Identification of chilling-responsive microRNAs and their targets in vegetable soybean (Glycine max L.). Sci. Rep. 6:26619. doi: 10.1038/srep26619
Ye, B., Saito, A., and Minamisawa, K. (2005). Effect of inoculation with anaerobic nitrogen-fixing consortium on salt tolerance of Miscanthus sinensis. Soil Sci. Plant Nutr. 51, 243–249. doi: 10.1111/j.1747-0765.2005.tb00028.x
Zhang, C. B., Wang, J., Ke, S. S., and Jin, Z. X. (2009). Effects of natural inhabitation by Miscanthus floridulus on heavy metal speciations and function and diversity of microbial community in mine tailing sand. J. Plant Ecol. 33, 629–637. doi: 10.3773/j.issn.1005-264x.2009.04.001
Zhou, L., Chen, J., Li, Z., Li, X., Hu, X., Huang, Y., et al. (2010). Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One 5:e15224. doi: 10.1371/journal.pone.0015224
Zhou, Z. S., Zeng, H. Q., Liu, Z. P., and Yang, Z. M. (2011). Genome-wide identification of Medicago truncatula microRNAs and their targets reveals their differential regulation by heavy metal. Plant Cell Environ. Int. 35, 86–99. doi: 10.1111/j.1365-3040.2011.02418.x
Zhu, J., Li, W., Yang, W., Qi, L., and Han, S. (2013). Identification of microRNAs in Caragana intermedia by high-throughput sequencing and expression analysis of 12 microRNAs and their targets under salt stress. Plant Cell Rep. 32, 1339–1349. doi: 10.1007/s00299-013-1446-x
Keywords: chromium, high-throughput sequencing, microRNA, abiotic stress, Miscanthus sinensis
Citation: Nie G, Liao Z, Zhong M, Zhou J, Cai J, Liu A, Wang X and Zhang X (2021) MicroRNA-Mediated Responses to Chromium Stress Provide Insight Into Tolerance Characteristics of Miscanthus sinensis. Front. Plant Sci. 12:666117. doi: 10.3389/fpls.2021.666117
Edited by:
Eduardo Gusmão Pereira, Universidade Federal de Viçosa, BrazilReviewed by:
Milan Borisev, University of Novi Sad, SerbiaVivek Dogra, Institute of Himalayan Bioresource Technology (CSIR), India
Copyright © 2021 Nie, Liao, Zhong, Zhou, Cai, Liu, Wang and Zhang. 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: Xinquan Zhang, zhangxq@sicau.edu.cn
†These authors have contributed equally to this work