Abstract
Tea tree [Camellia sinensis (L.) O. Kuntze] is an important leaf
(sometimes tender stem)-using commercial plant with many medicinal
uses. The development of newly sprouts would directly affect the yield
and quality of tea product, especially significant for Pingyang
Tezaocha (PYTZ) which takes up a large percent in the early spring tea
market. MicroRNA (miRNA), particularly the conserved miRNAs, often
position in the center of subtle and complex gene regulatory systems,
precisely control the biological processes together with other factors
in a spatio-temporal pattern. Here, quality-determined metabolites
catechins, theanine and caffeine in PYTZ sprouts including buds (sBud),
different development stages of leaves (sL1, sL2) and stems (sS1, sS2)
were quantified. A total of 15 miRNA libraries of the same tissue with
three repetitions for each were constructed to explore vital miRNAs
during the biological processes of development and quality formation.
We analyzed the whole miRNA profiles during the sprout development and
defined conserved miRNA families in the tea plant. The differentially
expressed miRNAs related to the expression profiles buds, leaves, and
stems development stages were described. Twenty one miRNAs and eight
miRNA-TF pairs that most likely to participate in regulating
development, and at least two miRNA-TF-metabolite triplets that
participate in both development and quality formation had been
filtered. Our results indicated that conserved miRNA act boldly during
important biological processes, they are (i) more likely to be linked
with morphological function in primary metabolism during sprout
development, and (ii) hold an important position in secondary
metabolism during quality formation in tea plant, also (iii) coordinate
with transcription factors in forming networks of complex multicellular
organism regulation.
Keywords: conserved miRNA, sprouts development, quality formation,
transcription factors, Camellia sinensis (L.) O. Kuntze
Introduction
Originally produced in China, green tea nowadays is the star among the
top list beverages, attribute to its good taste, health benefits, and
mysterious process, which brings considerable economic benefit in
planting and exporting countries such as China, India, Kenya, and Sri
Lanka. Based on reports from the China Tea Marketing Association 2017
([32]http://www.stats.gov.cn/), approximately 10.3 million tons of
fresh tea leaves were harvested to produce various tea products. As
young leaves and tender stems from the tea tree are processed to
prepare “tea”, the developmental characters are supposed to have a
direct and significant bearing on the yield and the quality of tea
product.
MicroRNAs (miRNAs), are endogenous single-stranded non-coding small
RNAs, that could both regulate their target messenger RNAs (mRNAs) at
chromatin state and could also perfectly or imperfectly bind to their
targets for further translation suppression by cleaving at some
complementary site (Rubio-Somoza and Weigel, [33]2011; Zheng et al.,
[34]2015). The plant miRNA families were thus placed to be at the
central position within gene expression programs, always with small
numbers per cell and large amounts of transcripts (Voinnet, [35]2009),
yet have powerful effect in developmental regulation, morphogenesis,
stress responses (Axtell and Bowman, [36]2008; De Lima et al.,
[37]2012; Jones-Rhoades, [38]2012; Yang et al., [39]2013). Still,
seldom gene regulation could be completed without consideration of the
transcription factors. Many transcription factors in the plant kingdoms
are highly conserved even stride over large evolutionary distances, and
for some of them, they could still share similar developmental roles in
diverse species (Zhao et al., [40]2013; Xu et al., [41]2016).
Hypothesis demonstrated that the miRNA binding sites evolve faster than
the transcription factor binding sites, as the ways to repress a gene
is relatively much more than to activate one (Chen and Rajewsky,
[42]2007).
However, some miRNA seems extremely well conserved (Lu et al.,
[43]2006). A large portion of the conserved miRNAs and their
conventional target transcription factors as well as F-box proteins
play pivotal roles in governing plastic behavior during development,
such as phase change and plant architecture (Kidner, [44]2010;
Rubio-Somoza and Weigel, [45]2011), making miRNA-TF mRNA pairs more
fascinating. At least 7 kinds of miRNAs were widely reported to
regulate in the three stages of leaf development and leaf morphology.
At the initiation stage, a division of leaf primordia are commonly
considered to be the key stage in the leaf development process, which
comes from a group of cells localized on the flanks of the shoot apical
meristem (SAM) loses their indeterminacy (Micol and Hake, [46]2003).
During this stage, miR390/ARF pathway has been described in the
regulation of leaf polarity (Braybrook and Kuhlemeier, [47]2010).
miR165/166 regulates the leaf polarity by targeting the HD-ZIP genes
and thus control the adaxial cell fate (Rubio-Somoza and Weigel,
[48]2011; Sun, [49]2012). Recent discovery revealed that the leaf
dorsoventral polarity (adaxial-abaxial) signals which may cause
mechanical heterogeneity of the cell wall, is linking to the
methyl-esterification of cell-wall pectins in tomato and Arabidopsis
(Qi et al., [50]2017). The shape and architecture of leaf need the
orchestration of auxin, KNOX genes and miRNA regulation. KNOX genes
could be down-regulated by CUC transcriptional regulators, which are
important for organ boundaries building (Takada and Tasaka, [51]2002;
Chen, [52]2009), floral patterning, and leaf morphogenesis (Micol and
Hake, [53]2003; Engstrom et al., [54]2004). NAC (NAM, CUC1/2-like) is
one branch of CUC gene family regulated by miR164. MiR164/GOB (a CUC2
ortholog gene), well-studied in tomato, is necessary for controlling
leaf polarity and determining the serration or smooth of the leaf
boundaries (Berger et al., [55]2009). MiR319, encoding by three loci
including miR-JAW (miR319a) in Arabidopsis, regulates five TEOSINTE
BRANCHED/CYCLOIDEA/PCF (TCP) family members (Palatnik et al., [56]2003,
[57]2007), which could also lead to the regulation of CUC genes.
Overexpression of miR319 or loss function of these five TCP genes would
result in crinkly leaves (Palatnik et al., [58]2003; Liu et al.,
[59]2018). TCP regulated growth and senescence via jasmonic acid
synthesis pathway (Schommer et al., [60]2008). The cell number and cell
size, which reported to be precisely spatial and temporal controlled
(Usami et al., [61]2009), are mainly regulated by an SQUAMOSA PROMOTER
BINDING PROTEIN PROTEIN-LIKE (SPL)-dependent pathway (Ferreira e Silva
et al., [62]2014; Xu et al., [63]2016). In Arabidopsis, miR156 targets
11 of the 17 SPL genes, among which SPL3, 4, and 5 accelerates the
juvenile-to-adult phase change, SPL9 and SPL15 regulate plastochron
length (Wang et al., [64]2008; Wu et al., [65]2009; Xu et al.,
[66]2016). MiR396 plays an important role in plant leaf growth and
development, most likely by repressing Growth-Regulating Factor (GRF)
genes in Arabidopsis. Transgenic miR396-overexpressing plants have
narrow-leaf phenotypes due to a reduction in cell number (Liu et al.,
[67]2009).
During the long cultivation history for more than 2,000 years in China,
numerous elite tea varieties have been bred for different
characteristics like early germination, high yield, good performance
under environmental stress, and distinctive aroma or flavor. Camellia
sinensis (L.) O. Kuntze “Pingyang Tezaocha” (PYTZ), an elite cultivar
with short internodes selected in Zhejiang Province in the late
century, is now popularized in the north tea area in China attributes
to its high yield for about 3tons green tea products per hectare [Data
from e-China tea from Tea Research Institute of China Academy of
Agriculture Sciences AS (TRI, CAAS)]
([68]http://www.e-chinatea.cn/other_shujuku.aspx) (Zhao et al.,
[69]2017). What's more, its early germination in April helps taking up
a large percent in the early spring tea market annually (Yang,
[70]2015). In tea plant, phenolic compounds is one of the most
important secondary metabolites, accounting for 18% to 36% dry weight
in the fresh leaves and tender stem (Jiang et al., [71]2013), is also
the main flavor components and functional ingredients that had been
intensely studied in the past decades for its effective and extensive
pharmacological activities (Zhao et al., [72]2013). Accordingly, the
accumulation of some phenolic-like nongalloylated catechins like
epigallocatechin (EGC) and epicatechin (EC) (Zhao et al., [73]2017),
quinic acid and flavonol glycosides are gradually increasing along with
the developing stages (Jiang et al., [74]2013). What's more, some key
enzyme genes involved in the biosynthetic pathway of phenolic compounds
in different organs and leaves at different developmental stages also
have the similar expression patterns, such like CsDHQ/DHS2 (DHQ/DHS,
3-dehydroquinate synthase), CsCHS1 (CHS, chalcone synthase), CsUGT78E1
(UGT, uridine diphosphate Glycosyltransferase), Cs4CL1
(4CL,4-coumaroyl-CoA ligase), CsF3′H1 (F3′H, flavonoid 3′-hydroxylase),
and some TF genes like Sg4 of CsMYB family (Jiang et al., [75]2013; Li
et al., [76]2017a), CsMYB5-1 and bHLH24-3 (Jiang et al., [77]2013).
Beyond catechins in tea plant, theanine and caffeine are the other two
characteristic constituents determine tea quality (Xia et al.,
[78]2017). No matter it is primary or secondary, there is no doubt that
metabolisms synchronize along with plant growth and development, are
under the precise entire spatio-temporal network control (Chen and
Rajewsky, [79]2007).
Here, the content of three kinds of main taste compounds catechins,
theanine, and caffeine in different tissues of spring sprouts including
bud, two stages of leaves, and stems of PYTZ were quantified by
High-Performance Liquid Chromatography-Mass Spectrometry (HPLC-MS/MS).
MiRNA libraries of the same tissues were constructed by Illumina HiSeq
technology in order to explore how miRNA works between development and
quality formation. Key miRNAs involved in regulating sprout development
have been speculated based on computational expression. Conserved miRNA
families in tea plant were obtained and mainly studied. To what extent
the conserved miRNAs might be linked with morphogenesis function during
sprout development was further investigated through other six
morphologically-different tea cultivars. Regulations in metabolic
pathways of conserved miRNA together with their target genes,
especially transcription factor genes that would finally determine tea
quality have been studied and discussed. The consistency of performance
between development and quality need more cross understanding and
balance in the subsequent process of screening tea cultivars.
Materials and Methods
Plant Material and RNA Isolation
The four-year-old tea plant cultivar Camellia sinensis (L.) O. Kuntze
“Pingyang Tezaocha” (PYTZ) were planted in the Germplasm of Qingdao Tea
Repository at the Tea Research Institute located in Qingdao (35°N119°E,
Qingdao city, China) under natural light condition. To ensure the
successiveness of gene expression during the development of newborn
branch, we collected the samples of bud, leaves, and stems orderly
downwards from the top ([80]Figure 1A) (Shen et al., [81]2019). For
normalization, the buds about 3 cm long (sBud), and the first leaf
below the bud (sL1) about 3.5cm long, the second leaf with higher
maturity below the bud (sL2) 4.5cm long, the stem between the first
leaf and the second leaf (sS1) with about 1.5cm long, and the more
mature stem between the second leaf and the fish leaf (sS2) with 2.1cm
long were measured and collected. For collecting samples of RNA,
healthy buds, leaves and stems at different developmental stages were
collected and frozen immediately in liquid nitrogen and stored in −80°C
freezers before use (Fan et al., [82]2015). Three biological replicates
were collected and pooled from at least five individuals, and each
biological replicate contained more than five buds, leaves and stems.
The total RNA for each sample was extracted using TRIzol reagent
(Invitrogen, Burlington, ON, Canada). The quality, purity,
concentration, and integrity of the total RNA was checked using 1%
agarose gel electrophoresis, NanoDrop Photometer Spectrophotometer
(IMPLEN, Westlake Village, CA, USA), Qubit RNA Assay Kit in Qubit 2.0
Fluorometer (Life Technologies, Carlsbad, CA, USA), and RNA Nano 6000
Assay Kit of the Bioanalyser 2100 system (Agilent Technologies, Santa
Clara, CA, USA), respectively. RNA samples with a 260/280 ratio between
1.8 and 2.0, 260/230 ratio between 2.0 and 2.5, and RNA integrity
number more than 8.0, were used for sequencing and quantitative PCR
analysis described below.
Figure 1.
[83]Figure 1
[84]Open in a new tab
The contents of characteristic metabolites. (A) Samples of bud, leaves,
and stems were collected orderly downwards from the top, named by sBud,
sL1, sL2, sS1, and sS2, from buds, leaves, and stems, correspondingly.
(B) The content of total catechins, caffeine, and theanine. (C–J) The
content of GC, EGC, C, EGCG, EC, ECG, caffeine, and theanine in each
tissue detected by HPLC-MS/MS. Standard catechins were
(-)-gallocatechin (GC), (-)-epigallocatechin (EGC), (+)-catechin (C),
(-)-epigallocatechin gallate (EGCG), (-)-epicatechin (EC), and
(-)-epicatechin gallate (ECG). (Above each column, the different low
case letters mean significant difference, p < 0.05, while labeled with
the same ones mean insignificance, p > 0.05. p-values were calculated
with Student's t-test).
Extraction and Quantification of Catechins, Caffeine, and Theanine
The extraction of catechins and caffeine followed by a previously
described method with minor modifications (Jiang et al., [85]2013; Wang
et al., [86]2018): 0.2 g of each fresh samples (sBud, sL1, sL2, sS1,
and sS2) were ground in liquid nitrogen and extracted with an
extraction solution (80% methanol and 20% water), followed by vortexing
and sonicating for 30 min at a low temperature. Then, the samples were
centrifuged at 3,500 g for 15 min, and the residues were re-extracted
twice as mentioned above until the final volume of the pooled
supernatants was 2 mL. The supernatants were then extracted three times
with chloroform and three times with ethyl acetate. The pooled
supernatant was concentrated to remove the ethyl acetate at a low
temperature with a vacuum pump. Finally, the product was dissolved in
200 μL methanol for quantification. The theanine was extracted as
reported by the method of Jeon et al. ([87]2017) with some
modifications. One gram of each finely ground sample was mixed with 100
mL boiling distilled water and brewing for 10 min (with the help of
magnetic stirrer). All the obtained extract were filtered by 0.45 μ
nylon membrane (after cooling down) and approximately 1 ml of the
sample solution were centrifuged at 13,000 rpm for 10 min prior to HPLC
analysis.
The isolation and detection of quality-related metabolites catechins,
caffeine, and theanine in the sprouts of PYTZ were performed by high
performance liquid chromatography-mass spectrometry (HPLC-MS/MS). HPLC
analyses were performed on an Agilent 1298 LC system (Agilent, Santa
Clara, CA, USA), and MS/MS detection was carried out using an Agilent
6460 Series Triple Quadrupole instrument (Agilent). Caffeine, theanine
and six major tea standards: catechins, (-)-epigallocatechin gallate
(EGCG), (-)-epigallocatechin (EGC), (-)-epicatechin gallate (ECG),
(-)-epicatechin (EC), (-)-gallocatechin (GC), and (+)-catechin (C) were
purchased from Sigma (St Louis, MO, USA). An Agilent 20RBAX RRHD
Eclipse Plus C18 column (particle size: 1.8 mm, length: 100 mm, and
internal diameter: 2.1 mm) was used at a flow rate of 1 mL min^−1. For
catechins and caffeine, the mobile phase consisted of 0.4% acetic acid
in water and 100% acetonitrile; and the gradient of latter increased
linearly from 0 to 10% (v/v) within 5 min, and to 35% at 20 min, to 10%
at 21 min, to 1% at 25 min. For theanine, the mobile phase consisted of
HPLC water and acetonitrile; and the gradient of former remained at
100% within 10 min, and decreased linearly from 100% to 20% (v/v) to 12
min, and kept at 20% to 20 min, to 100% at 22 min, and kept at 100% to
40 min. Mass spectra were acquired simultaneously using electrospray
ionization in the positive and negative ionization modes over the range
of m/z 100 to 2000. A drying gas flow of 6 L min^−1, drying gas
temperature of 350°C, nebulizer pressure of 45 psi, and capillary
voltages of 3,500 V were used. The compounds were identified
qualitatively using LC-MS by comparing the retention times (t[R]),
wavelengths of maximum absorbance (λ[max]), protonated/deprotonated
molecules ([M+H]^+/[M–H]^−), and major fragment ions with those of the
authentic standards and published literature (Jiang et al., [88]2013;
Jeon et al., [89]2017; Wang et al., [90]2018).
Library Construction and Small RNA Sequencing
For sRNA library construction, 3 μg of total RNA per sample was used
for the RNA sample preparations. Sequencing libraries were generated
using NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB,
USA). The library preparations were sequenced on an Illumina HiSeq™
2500 sequencer, by Gene Denovo Biotechnology Co. (Guangzhou, China).
The generated 50 bp single-end reads were then filtered out the impure
sequences (adaptor sequences and the low quality reads) and removed
cellular structural RNAs such as rRNA, snoRNA, snRNA, and tRNA based on
the alignment with small RNAs in GeneBank database (Release 209.0) and
Rfam database (11.0). The clean reads were mapped to the tea tree
genome without mismatch to analyze their expression and distribution
(NCBI Sequence Read Archive Database under accession [91]PRJNA381277).
Tags that mapped to exons or introns and repeat sequences were also
removed.
Identification of Known miRNAs and Novel miRNA
Since tea miRNA dataset was not included in the miRBase, the clean tags
were subjected to a Blastn search against miRBase 21.0, to identify and
annotate known miRNAs from all other plant miRNAs, allowing two
mismatches. All the known miRNAs were further checked for the existence
through 72 plant species, to figure out their conservative property.
The unannotated tags were aligned with tea tree genome to identify
novel miRNA candidates according to their genome positions and hairpin
structures predicted by software Mireap
([92]https://github.com/liqb/mireap, version 0.20).
miRNA Expression Profiles and Prediction of Target mRNAs
The expression levels of both known miRNA and novel miRNA from each
sample were calculated and normalized to transcripts per million (TPM)
(Wu et al., [93]2017b). The formula is TPM = Actual miRNA counts /
Total counts of clean tags^*10^6. Meanwhile, the correlation
coefficient between every two replicas was calculated to evaluate
repeatability between samples. Differential expression analysis across
samples was performed using the DEGseq (2010) R package. miRNAs with p
< 0.05 and log2-fold change ≥ 2 in comparison were set as the threshold
for significantly differentially expressed miRNAs (DEM). Candidate
target genes were predicted by using software PatMatch (Version 1.2)
blasting against tea tree genome, abiding by some rigorous parameters
as follows: No more than four mismatches between sRNA/target (G-U bases
count as 0.5 mismatches); For the miRNA/target duplex (5′ of miRNA),
(a) no more than two adjacent mismatches, (b) no adjacent mismatches in
positions 2–12, (c) no mismatches in positions 10–11, (d) no more than
2.5 mismatches in positions 1–12, and the minimum free energy (MFE) of
the miRNA/target duplex should be no < 60% compared to the MFE of the
miRNA bound to its perfect complement (Yan et al., [94]2005; Wu et al.,
[95]2017a).
Functional Enrichment Analysis of Target mRNAs
Gene Ontology (GO) enrichment analysis and KEGG pathway analysis were
performed to the target mRNAs of DEM in order to comprehensively figure
out their biological functions. All DEM target genes were mapped to GO
terms in the Gene Ontology database ([96]http://www.geneontology.org/),
then the enriched significant GO terms (taking FDR ≤ 0.05 as a
threshold, derived from calculated p-value) comparing to tea tree
genome background were categorized into three levels, “biological
process,” “cellular component” and “molecular function”. KEGG is the
major public pathway-related database (Kanehisa et al., [97]2008) for
further understand how genes interact with each other to play roles in
certain biological functions. The calculating formula is the same as
that in GO analysis. KEGG pathway enrichment analysis identify
significantly enriched metabolic pathways or signal transduction
pathways (Liu et al., [98]2014). Some online platforms or commercial
services that based on same or different mathematical algorithms could
help us to reconstruct gene networks, for example STRING
([99]https://string-db.org/), Pathway Commons
([100]https://www.pathwaycommons.org/) (Luna et al., [101]2016),
ANDSystem (Ivanisenko et al., [102]2019), and so on (Saik et al.,
[103]2018). Functional enrichment of both target genes of miRNAs in
single samples and DEM in a compare group were carried out in our
analysis. Here, STRING was used to show the enrichment networks.
Trend analysis was aiming at the expression of all miRNAs performing in
continuous tissues samples (Bud/L1/L2 and Bud/S1/S2) to cluster genes
with similar expression patterns. Trend analysis was carried out by
software Short Time-series Expression Miner (Ernst and Bar-Joseph,
[104]2006) (under parameters -pro 20 -ratio 1.0). GO and KEGG pathway
enrichment analysis was then be done to target genes of miRNAs in each
trend, and the p-value was obtained by hypothesis testing. Those GO
term and KEGG pathway were defined as significant ones satisfying Q
value ≤ 0.05. Q value was that p-value corrected by FDR (Benjamin and
Hochberg, [105]1995).
Quantitative PCR for miRNAs and mRNAs
The expression profiles of mature miRNAs and the potential target mRNAs
were further validated by quantitative PCR. Synthesis of the first
strand cDNA was performed with Mir-X™ miRNA First-Strand Synthesis Kit
(Cat. No. 638313, Clontech Laboratories, Inc., CA, USA), with 5.8S rRNA
served as an internal control. The first strand cDNA of mRNA were
synthesized by using PrimeScript™ RT reagent Kit with gDNA Eraser
(Perfect Real Time) (Code No. RR047A, Takara, Tokyo, Japan), with
glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene for
normalization. The primers of miRNAs and mRNA were listed in
[106]Supplementary Table 1. Quantitative PCR was carried out with SYBR
Premix Ex Taq™ IIKit (Tli RNase H Plus) (Code No. RR820A, Takara,
Tokyo, Japan), on a LightCycler 480 instrument (Roche Molecular
Systems, Inc., Indianapolis, IN, USA). The amplification program of
miRNA was performed under the following parameters: 95°C for 10 min, 40
cycles at 95°C for 15 s, 60°C for 1 min (Zheng et al., [107]2015). The
amplification program of mRNA was performed at 94°C for 10 s, 58°C for
10 s and 72°C for 10 s (Li et al., [108]2017b). Triplicates of each
reaction were performed, and 5.8S rRNA and GAPDH were used as
endogenous control separately. CT values obtained through quantitative
PCR were analyzed using 2^−ΔΔCT methods to calculate relative fold
change values.
In addition, fresh spring sprouts were plucked from seven tea varieties
from one experimental tea garden of Tea Research Institute, Fujian
Academy of Agricultural Sciences in Fu'an, China (27°10′N, 119°35′E):
Camellia sinensis “Jinfenghuang” (JFH), Camellia sinensis “Pingyang
Tezaocha” (PYTZ), Camellia sinensis “Zhengdayin” (ZDY), Camellia
sinensis “Dayewulong” (DYWL), Camellia sinensis “Huangdan” (HD),
Camellia sinensis “Jiukeng 6” (JK), Camellia sinensis “Queshe” (QS).
The bud, leaves, and stems were also sampled at the same position
mentioned above for quantitative PCR analysis. The result of relative
expressions was presented as a heatmap by using TBtools (Chen et al.,
[109]2018).
Transcription Factor Prediction and miRNA-mRNA-Metabolite Network
Construction
All the mRNA genes target by DEM predicted above were blasted against
the transcription factor (TF) database from the plant
([110]http://planttfdb.cbi.pku.edu.cn/, version 4.0) to annotate
potential TFs. The resulting target TF genes were classified in each
tissue for analysis. The resulting TF genes were further blasted
against the tea tree genome (Xia et al., [111]2017) (NCBI Sequence Read
Archive Database No. PRJNA381277). The methods for TF blasting and
expression analysis were followed by Zhao et al. ([112]2017). Different
expression profiles were finally grouped into profile 1, 2, and so on.
For the network analysis was based on Savoi's method (Savoi et al.,
[113]2016), the mRNA and metabolite association were obtained based on
Pearson correlation coefficient between the contents of each metabolite
and the expression levels of mRNA, and filtered the pairs when the
absolute value of cor was larger than 0.9 and p-value was smaller than
0.05. The miRNA and mRNA association were obtained by Spearman's
correlation coefficient according to their expression levels calculated
by TPM and filtered the pairs when cor-value was no larger than −0.5,
and the p-value was smaller than 0.05. The network was visualized by
Cytoscape (V3.6.0) (Praneenararat et al., [114]2012).
Availability of Supporting Data
Clean Illumina sequencing reads of 15 small RNA of PYTZ sprout have
been deposited in the NCBI Sequence Read Archive Database under
accession [115]PRJNA510482.
Results
The Contents of Quality-Related Metabolites in Tea Sprouts
The buds, the first leaves, the first stems, sometimes with the second
leaves together, are the most common raw materials for producing green
tea. The contents of the most quality-related metabolites in tea
sprouts that contributing flavors and health-promoting functions were
quantified by HPLC-MS/MS ([116]Figures 1B–J). Among the three kinds of
characteristic metabolites, the content of total catechin accounts for
the most majority proportion, while distributing significantly
different among tender leaves and stems ([117]Figure 1B). The
galloylated catechins such like EGCG ([118]Figure 1F) and ECG
([119]Figure 1H), take up 67.4% in total catechin concentration
obtained by summation of the individual components ranged from 37.3 to
119.2 mg/g ([120]Figure 1B). Both EGCG and ECG had a relatively low
concentration in sS2. The level of GC ([121]Figure 1C), C ([122]Figure
1E), and EC ([123]Figure 1G), were found to be low in bud and leaves,
especially for C. The concentrations of caffeine ranging from 4.7 to
9.7 mg/g, accumulated to its highest levels in sBud and lowest in sS2
([124]Figure 1I). Theanine, however, showed its highest accumulation in
stems in PYTZ sprouts ([125]Figure 1J).
Overview of microRNA Profile and Its Mapping to Tea Tree Genome
To figure out what role the miRNA play during the formation of
characteristic metabolites along with the development of the PYTZ
sprout, 15 sRNA-Seq libraries including buds, leaves, and stems were
separately sequenced on Illumina HiSeq™ 2500 platform generating a
total of 292,653,360 raw reads. After removing dirty reads containing
adapters and low quality bases, in average, the clean tags of
14,080,519 for sBud (bud), 12,776,739 for sL1 (the 1st leaf),
11,830,623 for sL2 (the 2nd leaf), 11,484,341 for sS1 (the younger
stem), and 9,822,915 for sS2 (the older stem) were retained. The
filtering data of each procedure were listed in [126]Supplementary
Table 2. Most clean tags had the length of 21-24nt, in which the 24 nt
sRNAs were the most abundant ([127]Figure 2A). The proportion of
different length tags has no obvious difference among the five sample
groups, and generally showed the trends of increased and then decreased
bounded by the 24 nt sRNAs. Notably, the number of 24 nt sRNAs in sS1
was the lowest, while in sL1 was the highest.
Figure 2.
[128]Figure 2
[129]Open in a new tab
Overview of PYTZ microRNA libraries. (A) Size-class distribution of
sRNAs in the sBud, sL1, sL2, sS1, and sS2 libraries. Designations the
same as on previous figure. (B) Number statics of tissue-specific
miRNAs. (C) Number statics of differentially expressed miRNAs (DEM)
between each two sample groups. (D–G) The number of up/down-regulated
DEM between different tissues in G1/G2.
About 76.97% clean tags were perfectly mapped to tea tree genome (NCBI
Sequence Read Archive Database No. PRJNA381277), which indicated a
credible quality of sequencing, and the rates of a genomic match were
similar across these samples. We removed the tags mapped to exons
located in positive-sense strands which might be fragments from mRNA
degradation. (For statistics of mapping to tea tree genome, see
[130]Supplementary Table 3). The tags mapped to repeat sequences were
also excluded. The clean tags were then aligned with the Rfam database
(11.0) and the percentage of annotation was summarized in
[131]Supplementary Table 4. The average of rRNA, snRNA, snoRNA, and
tRNA in the samples took up of 18.11, 0.19, 0.56, and 1.01%.
Known miRNAs Identification and Novel microRNAs Prediction in PYTZ
To identify known miRNAs in tea, all the unannotated unique tags were
blast-searched against plant miRNAs in miRBase (Release 21.0, June
2014). Overall, a total of 1,928,678 miRNA clean tags were identified
from 15 libraries ([132]Table 1), and 156 known miRNAs were identified
([133]Supplementary Table 5). Most of the identified known miRNAs
(81.97%) belonged to the 21 nt length miRNA families, the remaining
ones belonged to 18–24 nt miRNAs families ([134]Supplementary Table 5).
For the 156 known miRNAs, 122 precursors and 99 kinds of characteristic
hairpin structures of the known miRNAs were identified. The length of
the precursors varied from 71 to 288 nt, with an average of 144 nt, and
the average minimum free energy (MFE) was −57.23 kcal/mol, ranging from
−22.3 to −90.2 kcal/mol ([135]Supplementary Table 6). The abundance of
miRNA including the novel miRNA showed a high distribution in leaf than
in stem (ratio value in [136]Table 1), with that in S1 the lowest.
Table 1.
Numbers and ration of clean tags for conserved miRNA and novel miRNA.
Sample Total Known miRNA Novel miRNA
Mirna num Tags uniq Tags total/ratio Mirna num Tags uniq Tags
total/ratio
sBud-1 13025496 90 2281 145547(1.12%) 986 1205 35371(0.27%)
sBud-2 13672142 87 2150 96863(0.71%) 780 1003 22599(0.17%)
sBud-3 15543919 86 2172 79234(0.51%) 669 876 23593(0.15%)
sL1-1 13217393 100 2571 237953(1.80%) 1000 1244 47097(0.36%)
sL1-2 12086276 98 3140 205017(1.70%) 1008 1276 36970(0.31%)
sL1-3 13026547 91 1967 127242(0.98%) 854 1066 30010(0.23%)
sL2-1 11146588 94 2523 151824(1.36%) 820 1037 26207(0.24%)
sL2-2 13556636 90 2518 162349(1.20%) 737 872 21904(0.16%)
sL2-3 10788646 90 2485 201288(1.87%) 783 1005 27762(0.26%)
sS1-1 10821241 73 1763 48215(0.45%) 508 693 11804(0.11%)
sS1-2 11025165 76 1753 48019(0.44%) 483 669 10145(0.09%)
sS1-3 12606617 81 1902 74064(0.59%) 659 852 18301(0.15%)
sS2-1 9125465 81 2079 114250(1.25%) 745 982 25010(0.27%)
sS2-2 11303334 90 2451 130149(1.15%) 847 1025 26705(0.24%)
sS2-3 9039947 79 2144 106664(1.18%) 781 945 21290(0.24%)
total 179985412 1306 33899 1928678(16.31%) 11660 14750 384768(3.25%)
[137]Open in a new tab
For the four nucleic acids, the frequency of cytosine (C) (32.09%) and
uracil (U) (29.65%) is higher than guanine (G) (19.26%) and adenosine
(A) (19.00%). In the five samples, U had a high appearance at the 1st,
17th, 22th, and 23rd positions, with an average of 84.34, 61.80, 54.57,
and 52.48%, respectively ([138]Supplementary Figure 1). C occupied a
very high percentage (87.44%) at 19th position. The analysis showed
that A had a relatively high proportion at 9th position in the stem
(sS1 and sS2) and 17th position in leaf (sL1 and sL2), in contrast, A
was seldom present at 2nd, 13th, 18th, and 20th positions in the five
tissues. For the first nucleotide bias analysis, U had the absolute
predominance in miRNAs with the length of 20, 21, and 22 nt
([139]Supplementary Figure 2).
The remaining reads which couldn't get mapped to known miRNAs were used
to identify novel miRNAs. 384,768 novel miRNA tags were identified from
15 libraries ([140]Table 1), and 1186 novel miRNA tags were identified
by predicting the hairpin structures of their precursor sequences
([141]Supplementary Table 7). The length of the novel miRNAs ranged
from 18 to 27 nt, different with known miRNAs, the 22 nt length miRNA
families were the most abundant (46.71%), followed by 21nt (41.23%).
These novel miRNA were involved in 1130 hairpin miRNA precursors. The
length of these precursors varied from 65 to 373 nt, with an average of
178 nt. The average minimum free energy (MFE) was −56.26 kcal/mol,
ranging from −18.1 to −292.3 kcal/mol ([142]Supplementary Table 7). The
numbers of novel miRNAs were most in sL1 and lowest in sS1, the trend
was the same with known miRNAs.
Whole miRNA Expression Characters in PYTZ
To figure out the whole miRNAs express patterns in the spring sprouts
of PYTZ, we need to evaluate the reliability of parallel experimental
results as well as operational stability. The expression level of all
miRNAs including known miRNAs and novel miRNAs from 15 libraries was
normalized to generate TPM, which further used to compute the related
coefficients. The strong correlation between every two biological
replicates for interlibrary of all five sample groups brought out that
the sequencing results are highly reliable ([143]Supplementary Figure
3). The correlations between sL1 and sL2, sS1, and sS2 were
substantially higher than other inter-groups, suggesting closely
associated integral processes in the separate development of leaf and
stem.
To lock the target miRNAs which might be responsible for the tea shoot
development, we firstly define sBud/ sL1/sL2 as Group 1 (G1), sBud/
sS1/sS2 as Group 2 (G2) to see the whole miRNA change characters. In
G1, a total of 226 miRNAs showed different expression among the three
tissues and classified into 8 profiles according to their trends.
Overall, the down-expression trend miRNAs took up a larger percentage
(69.5%), which belonged to profile 3 (75 miRNAs), profile 0 (67
miRNAs), and profile 1 (15 miRNAs). The up-expression trend miRNAs
which had the most abundant expression in sL2 belonged to profile 6 (40
miRNAs), profile 7 (12 miRNAs), and profile 4 (1 miRNA). There were
other 9 and 7 miRNAs which had the highest and lowest expression level
in sL1, separately ([144]Figure 3A).
Figure 3.
[145]Figure 3
[146]Open in a new tab
Expression profile trends of the whole miRNA in G1 (A) and G2 (B).
In G2, a total of 273 miRNAs showed different expression among the
three tissues and had been classified into 8 profiles according to
their trends. Overall, 91 miRNAs showed the down-expression trend
belonging to profile 3 (42 miRNAs), profile 1 (34 miRNAs), and profile
0 (15 miRNAs). 64 miRNAs showed the up-expression trend belonging to
profile 6 (35 miRNAs), profile 4 (25 miRNAs), and profile 7 (4 miRNA).
It is noteworthy that there were 95 miRNAs and 23 miRNAs showed the
lowest and highest expression level in sS1, separately ([147]Figure
3B).
Conserved miRNAs Families and Tissue-Specific miRNAs
The 156 known miRNAs belonging to 125 families, among which 27 families
were well-conserved that present in more than 10 plant species out of
72 plant species ([148]Table 2). miR156 was the most popular one, which
was found in 51 plant species, followed by miR396 and miR166, which
were conserved in 47 and 45 plant species, respectively.
Table 2.
27 conserved miRNA families from PYTZ searched in other 72 plant
species.
Species M IR 156 M IR 396 M IR 166 M IR 171 M IR 160 M IR 167 M IR 164
M IR 172 M IR 319 M IR 159 M IR 169 M IR 408 M IR 390 M IR 395 M IR 398
M IR 168 M IR 399 M IR 162 M IR 393 M IR 394 M IR 482 M IR 403 M IR 530
M IR 2111 M IR 477 M IR 535 M IR 2118
Brachypodium distachyon + + + + + + + + + + + + + + + + + + + + – – + –
– – +
Glycine max + + + + + + + + + + + + + + + + + + + + + + + + – – +
Arabidopsis thaliana + + + + + + + + + + + + + + + + + + + + – + – + –
– –
Oryza sativa + + + + + + + + + + + + + + + + + + + + – – + – – + +
Arabidopsis lyrata + + + + + + + + + + + + + + + + + + + + – + – + – –
–
Malus domestica + + + + + + + + + + + + + + + + + + + + + + – + + + +
Prunus persica + + + + + + + + + + + – + + + + + + + + + + + + + + –
Vitis vinifera + + + + + + + + + + + + + + + + + + + + + + – + + + –
Populus trichocarpa + + + + + + + + + + + + + + + + + + + + + + + + + –
–
Solanum tuberosum + + + + + + + + + – + + + + + – + + + – + – + – + – –
Citrus sinensis + + + + + + + + + + + + + + + – + + + + + + + – + + –
Cucumis melo + + + + + + + + + + + + + + + + + + + + – – + + + – –
Manihot esculenta + + + + + + + + + + + + + + – + + + + + + + + + + + –
Amborella trichopoda + + + + + + + + + + + – + + + + – – + + – – – + +
+ –
Medicago truncatula + + + + + + + + + + + + + + + + + + + – + – + + – –
+
Nicotiana tabacum + + + + + + + + + + + + + + + + + + – + + – – – + – –
Aegilops tauschii + + + + + + + + + – + + + + + + + – + + – – – – – – +
Zea mays + + + + + + + + + + + + + + + + + + + + + – – – – – +
Sorghum bicolor + + + + + + + + + + + + + + + + + + + + – – – – – – +
Theobroma cacao + + + + + + + + + – + – + + + + + + + + – + + + – + –
Carica papaya + + + + + + + + + + + + + + + – – + + + – – – – + + –
Linum usitatissimum + + + + + + + + + + + + + + + + + + + + – – + – – –
–
Ricinus communis + + + + + + + + + + + + + + + + + + + – – + – – – + –
Solanum lycopersicum + + + + + + + + + + + – + + – + + + – + + + – – +
– –
Triticum aestivum + + – + + + + – + + + + – + + – + – – – – – + – – – –
Aquilegia caerulea + + + + + + – + + + + + – + + + + – – – + – + – + +
–
Brassica napus + + + + + + + + – + + – + + – + + + + + – + – + – – –
Brassica rapa + + – + + + + + + + – + + + + + – + – – – + – + – – –
Cynara cardunculus + + – + + + + + + – + + + + + + + – + + – – – – – –
–
Physcomitrella patens + – + + + + – – + – – + + + – – – – – – – – – – +
+ –
Gossypium hirsutum + + + – + + + + – – + – + – + – + + + + + – – – – –
–
Hevea brasiliensis + + + – – – – – + + – + – – + – – – – – + – – – – –
+
Vigna unguiculata + – – – + – + + + – + + – + – + + + – – + – – – – – +
Gossypium raimondii – – + – – + + + – – – – – – + – + – – – + – + – + –
–
Salvia sclarea + + + + – – + + – – + – – + + – + – – + – – – – – – –
Pinus taeda + + + + – – – – + + – + + – + – – – – – + – – – – – –
Festuca arundinacea + + + + + – + – – + + – – – – – – – – – – – – – – –
–
Saccharum sp. + + + – – + – – – + + + – – – + – – – – – – – – – – –
Arachis hypogaea + – – – + + – – – + – + – – + – – – – + – – – – – – –
Helianthus tuberosus + – – + + – – – – + – – – – – – – + + – – + + – –
– –
Hordeum vulgare + – + + – – – – – + + – – – – + + – – – – – – – – – –
Lotus japonicus – + – + – + – + – – – + + – – – – – – – – – – + – – –
Pinus densata – + + + – – – – – + + – + – – – – + – – + – – – – – –
Selaginella moellendorffii + + + + + – – – + + – + – – – – – – – – – –
– – – – –
Digitalis purpurea + + + – + + – + – – – + – – – – – – – – – – – – – –
–
Picea abies – + + – + – – – – – – – – + – – – – – – + – – – – + –
Acacia auriculiformis – + – – + – – + + – – – – – – + – + – – – – – – –
– –
Citrus trifoliata + – + + – + + – + – – – – – – – – – – – – – – – – – –
Phaseolus vulgaris – – + – – – – – + + – – – – – – + – – – + – – – – –
+
Saccharum officinarum + + – – – + – – – + – + – – – + – – – – – – – – –
– –
Brassica oleracea – – – + – – – + – – – – – – + – – – – – – – – – – – –
Citrus clementina – + – + – + – – – – – – – – – + – – – – – – – – – – –
Citrus reticulata – – + + – – – – – – – – – – – + – – – – – – – – – – –
Cunninghamia lanceolata – – + – – – + – – – – – – – – – – + – – – – – –
– – –
Glycine soja – – – – – + – – – – – – – – – – – – – – + – – – – – –
Helianthus annuus + – – – + – – – – – – – – – – – – – – – – – – – – – –
Helianthus paradoxus + – + + – – – – – – – – – – – – – – – – – – – – –
– –
Helianthus petiolaris – – + – – – – – – – – – – – – – – + – – – + – – –
– –
Rehmannia glutinosa – – – – – – + – – – – – – – – – – – – – – – – – – –
–
Acacia mangium – + – – – – – – + – – – – – – – – – – – – – – – – – –
Avicennia marina + + – – – – – – – – – – – – – – – – – – – – – – – – –
Bruguiera cylindrica + + – – – – – – – – – – – – – – – – – – – – – – –
– –
Bruguiera gymnorhiza + + – – – – – – – – – – – – – – – – – – – – – – –
– –
Helianthus argophyllus + – – – – – – – – – – – – – – – – – – – – + – –
– – –
Helianthus ciliaris + – – – – – + – – – – – – – – – – – – – – – – – – –
–
Panax ginseng – – – – – – – – – – – – – – – – – – – – + – – – – – +
Chlamydomonas reinhardtii – – – – – – – – – – – – – – – – – – – – – – –
– – – –
Elaeis guineensis – – – – – – – + – – – – – – – – – – – – – – – – – – –
Gossypium herbaceum – – – – – – – – – – + – – – – – – – – – – – – – – –
–
Helianthus exilis – – – – – – – – – – – – + – – – – – – – – – – – – – –
Populus euphratica – – – – – – – – – – – – – – – – – – – – – – – – – –
–
Triticum turgidum – – – – + – – – – – – – – – – – – – – – – – – – – – –
Sum 51 47 45 43 41 40 37 37 37 36 36 34 33 33 33 32 32 31 26 26 23 17
16 15 15 13 12
[149]Open in a new tab
The expression of miRNA usually tells more about the code of regulating
new shoot elongation and development. In order to filter the
tissue-specific miRNAs in PYTZ, we merged the DEM from sL1 and sL2 into
sL, sS1 and sS2 into sS, and removed duplicates, separately. The Venn
diagram ([150]Figure 2B) showed that most miRNA were existed in all
tissues (82.49%) or at least in one tissue, regardless of their
relatively high or medium expression abundance. Interestingly, some
miRNA could only have their expression in specific tissues. Five miRNAs
were bud-specific that could only be expressed in the bud, 31 were
leaf-specific, and 28 were stem-specific, which had been summarized in
[151]Table 3.
Table 3.
Tissue specific miRNAs of PYTZ.
Tissue specific miRNAs of PYTZ
Bud miR1128-x, miR781-y, miR8590-y, miR870-y, miR9759-y
Leaf miR1042-x, miR1045-x, miR1057-y, miR1063-x, miR1515-x, miR164-y,
miR1866-y, miR2083-y, miR2111-y, miR2275-y, miR393-x, miR5061-y,
miR5181-y, miR5523-y, miR5538-x, miR5653-y, miR6173-y, miR6281-x,
miR7528-y, miR7711-x, miR7725-x, miR845-z, miR858-y, miR8665-y,
miR9569-x, novel-m0703-5p, novel-m0722-5p, novel-m0945-3p,
novel-m0953-5p, novel-m1089-3p, novel-m1125-5p
Stem miR1127-x, miR1865-x, miR2275-x, miR3512-y, miR3630-x, miR4388-y,
miR474-x, miR474-y, miR5021-x, miR5385-x, miR5834-x, miR6485-x,
miR7713-x, miR7717-x, miR7762-y, miR8007-y, miR8681-y, miR9863-x,
novel-m0024-3p, novel-m0095-3p, novel-m0155-3p, novel-m0331-5p,
novel-m0392-5p, novel-m0828-3p, novel-m0906-3p, novel-m1037-3p,
novel-m1040-3p, novel-m1041-5p
[152]Open in a new tab
Differentially Expressed miRNAs (DEM)
The differentially expressed miRNAs (DEM) were pairwise compared among
sBud, sL and sS, with their expression values higher than a 2-fold
change and p ≤ 0.05, aiming to find out key miRNAs during the
development. The numbers of DEM between tissues were summarized in
[153]Figures 2C–G. It's worth noting that according to the
developmental order, there was a sharp down-trend numbers of DEM (235)
and 79 up-trend DEM in sL1 than in sS1 (referred to as sL1 vs. sS1),
compared to 67 down-trend and 5 up-trend DEM in sL1 vs. sL2, reminding
that there are distinct regulatory changes and thus a metabolites
accumulate differences bounded between L1 and S1 ([154]Figure 2C). DEM
with ups and downs both in G1 ([155]Figures 2D,F) and G2 ([156]Figures
2E,G) were also classified. The continuously changing ones were seemed
to be possibly interesting regulators, such as 5 DEM (miR390-x,
novel-m0578-5p, novel-m0634-5p, novel-m0503-3p, novel-m0531-5p) that
had downtrend expression in G1 ([157]Figure 2F), miR396-x uptrend in G1
([158]Figure 2D) and novel-m0331-5p ([159]Figure 2E) uptrend in G2.
GO Enrichment and KEGG Pathway Analyses of DEM
miRNA sequences were searched against tea tree genomic sequences using
the plant miRNA potential target finder to predict target mRNAs. The
annotation of the target unigene of DEMs was conducted based on GO
enrichment and KEGG analyses. In this study, a total of 5501 potential
unigenes were predicted to be targeted by 934 miRNAs, including 138
conserved and 796 novel miRNAs. Among the miRNA, miR5385-x targeted the
most unigenes (770), followed by miR5658-x (533), and miR8577-x (223).
280 miRNAs targeted one unigene, while most miRNAs could target
multiple sites. Similarly, one unigene was also targeted by several
miRNAs, and there were 1351 unigenes could be regulated by more than
one miRNA, 58 of which were targeted by no <10 conserved miRNA. (The
complete list of target genes of all miRNA were listed in
[160]Supplementary Table 8).
Gene Ontology (GO) enrichment analysis offered a strictly defined
concept to describe properties of the target genes and recognize the
main biological functions in a dynamic-updated controlled vocabulary.
Within biological process categories, represented GO terms associated
with these target genes in all tissues ([161]Supplementary Figure 4)
were related to “metabolic process” the most, followed by “cellular
process” and “single-organismal process”. Within cellular components
categories, the unigenes were similarly represented, mainly in “cell,”
“cell part,” “membrane,” “organelle” and their parts. Within the
molecular function categories, the top two GO terms were “catalytic
activity” and “binding”.
Key DEM Involved in Growth and Development
We focused on the expressions of miRNAs to filter the possible ones
that participate in growth and development, and we found that all of
the reported growth and development associate miRNAs belong to the up-
or down- trends pattern, except miR172. So we firstly narrowed down
miRNAs with similar expression trends in the fore-mentioned G1 and G2
(in which the down-trends including profile 3, 0 and 1; the up-trends
including profile 6, 7 and 4), and then further screened by GO and KEGG
pathway analysis. Twenty one miRNAs, including 6 novel miRNAs were
screened out to be potential developmentally important miRNAs in PYTZ
([162]Table 4). Mature sequences of these miRNAs and their target genes
in tea genome were also listed in [163]Table 4. The heat map of the 21
miRNAs which represent their transcription levels calculated by TPM in
the samples were displayed in [164]Figure 4. Not each miRNA has the
same expression pattern in G1 compared with that in G2: 3 miRNAs only
changed in G1, with 1 uptrend and 2 downtrends; 4 miRNAs only changed
in G2, with 2 uptrends and 2 downtrends; 11 miRNAs had the same trends
in both G1 and G2, with 3 uptrends and 8 downtrends. Interestingly,
miR319-y had a downtrend expression pattern in G1 and uptrend in G2.
Some other miRNAs hadn't been included in these trends might not
because of their expression trends or levels, but the difference of
expressions among development stages was not significant (P < 0.05,
log2-fold change ≥2).
Table 4.
Potential developmentally important miRNAs in Pingyang tezaocha.
miRNA Expression pattern Mature miRNA sequence (5^′-3^′) Target mRNA ID
in tea genome
miR156-x G2 up CUGACAGAAGAGAGUGAGCAC CSA007311 CSA009012
[MATH: CSA011373_ :MATH]
CSA013149 CSA017838
[MATH: CSA019508_ :MATH]
CSA020439
[MATH: CSA023442_ :MATH]
[MATH: CSA031667_ :MATH]
miR160-x G1/G2 down UGCCUGGCUCCCUGUAUGCCA CSA021765 CSA026035 CSA026847
miR164-x G1/G2 up UGGAGAAGCAGGGCACGUGCA CSA013362 CSA027185 CSA033493
miR165-y G1 up UCGGACCAGGCUUCAUCCCCU
[MATH: CSA023057_ :MATH]
[MATH: CSA030874_ :MATH]
miR166-x G1/G2 down GGAAUGUUGGCUGGCUCGAUG CSA001544 CSA017082 CSA028252
CSA028940 CSA031362 CSA031363
miR166-y - UCGGACCAGGCUUCAUUCCCC CSA030874
miR166-z - UCGGACCAGGCUUCAUUCCCU CSA030874
miR319-x G1down GAGCUUCCUUCUGUCCACUU CSA003974 CSA005021 CSA007451
CSA014848 CSA017867 CSA018024 CSA018464 CSA021711 CSA036427
miR319-y G1down
G2 up UUGGACUGAAGGGAGCUCCCU CSA002826 CSA013833 CSA031080
[MATH: CSA036087_ :MATH]
miR390-x G1/G2down AAGCUCAGGAGGGAUAGCGCC CSA001776 CSA001956 CSA006028
CSA007412 CSA011969 CSA016385 CSA022179 CSA026282 CSA030630 CSA036196
miR396-x G1/G2 up UUCCACAGCUUUCUUGAACUU CSA003330
miR396-y G2 up GCUCAAGAAAGCUGUGGGAAG CSA008398
miR5083-y G1/G2 down CUACAAUUAUCUGAUCAAA CSA036373
miR8175-y G2 down UCCCCGGCAACGGCGCCA CSA013921 CSA018936
miR8577-x G1 down UGAGAUGAUGAUCAUGAU CSA008171
[MATH: CSA030921_ :MATH]
novel-m0675-3p G1/G2 down GAAUAUGAUGAAUUUGAAUG CSA006098
novel-m0297-3p G1/G2 down ACCCCUAACCCCAACACCCAAUC CSA012411
novel-m0243-3p G1/G2 down GAUCAGGAUGAAGCAACAUU CSA027591 CSA034664
novel-m0284-3p G1/G2 down UUGGGCUGGGCAGAAAUUGGGC CSA016454
novel-m0800-3p G1/G2 up GUUCAGUGAAGCUGUGGAAAG CSA010612 CSA028516
novel-m0187-3p G2 down AUUUCCCUUUCCAAAUUCCUU CSA012066 CSA024775
[165]Open in a new tab
The underlined genes are the targeted mRNA belonged to TF genes. The
gray shaded ones are the genes have reciprocal expression profiles with
responding miRNAs.
Figure 4.
[166]Figure 4
[167]Open in a new tab
Heat map of the 21 miRNAs by the value of TPM in the 15 samples. The
three major clusters on the left showed evolutionary relationships of
the miRNAs, with closely related miRNAs neighboring clustered. Similar
expression patterns could be found within one cluster.
Potential Transcription Factor Target Genes
As transcription factors were intensely studied in their numerous
important roles during plant growth and development in many species
(Ramachandran et al., [168]1994; Zhang et al., [169]2009; Chen et al.,
[170]2010), we are here supposed to analysis transcription factor genes
for identifying key TFs performing this function. All the 5,501
predicted target genes were blasted against Plant Transcription Factor
Database ([171]http://planttfdb.cbi.pku.edu.cn/), resulting in a total
of 46 kinds of transcription factors involving 352 mRNAs were detected
(The full list of identified TFs was provided in [172]Supplementary
Table 9). Then, types and numbers of TFs genes targeted by miRNAs in
each tissue (sBud, sL1, sL2, sS1, and sS2) were analyzed and summarized
in [173]Figure 5. On the whole, the numbers of transcription factors
genes targeted by miRNAs in sL2 were the most, and that in sS1 was the
least. MYB showed the widest involvement, followed by HD-ZIP and bHLH
(basic Helix-Loop-Helix) transcription factors. Some kinds of
transcription factors couldn't be targeted in all tissues, such as that
the CAMTA family of calmodulin binding TF genes in sS2 and GRF (Growth
Regulating Factor) in sL1.
Figure 5.
[174]Figure 5
[175]Open in a new tab
Types and numbers of transcription factors (TFs) genes targeted by
miRNAs in each tissue.
Quantitative PCR was further performed to validate the mentioned 21
growth and development associate miRNAs ([176]Figure 6) and their
predicted target mRNAs ([177]Figure 7), among which 14 of them were TF
genes (underlined in [178]Table 4). Universally, plant miRNAs might be
involved in many complicated and diverse functions in the complex
regulatory networks, the fundamental role of miRNAs is to suppress the
expression of target genes (Tang and Chu, [179]2017). Herein, we got
eight miRNA-TF genes pairs with reciprocal expression profiles (gray
shaded mRNAs in [180]Table 4), and the complementary correspondence of
miRNA toward the target sites were shown in [181]Supplementary Figure
5. The miRNA-TF genes pairs were miR156-x-CSA011373 (SBP),
miR156x-CSA019508 (SBP), miR156x-CSA023442 (SBP), miR156x-CSA031667
(SBP), miR165y-CSA023057 (Class III HD-Zip), miR165y-CSA030874
(HD-Zip), miR319y-CSA036087 (MYB), and miR8577x-CSA030921 (bHLH).
Figure 6.
[182]Figure 6
[183]Open in a new tab
Relative expressions of the 21 miRNA by quantitative PCR.
Figure 7.
[184]Figure 7
[185]Open in a new tab
Relative expressions of the predicted target mRNA genes by quantitative
PCR. Nodes on axis X for each inset image from left to the right are:
sBud, sL1, sL2, sS1, and sS2. Designations the same as on previous
figure.
Expression Profiles of Reported Morphological miRNA in Different Tea
Varieties
Originated in Yunnan and Tibet region of China, the tea tree has been
evolved over thousands of years and now at least 246 cultivars have
been selected breeding with significant differences in morphology and
physiology. Such as the focused cultivar in this study, PYTZ, has
oblong leaf shape, blunt tip, tight dentate margins, and shorter
internode. These characteristics are one of the indicators of screening
and distinguishing cultivars. Plant growth and development are
accompanied by morphogenesis that some regulators including small RNAs
and TFs may participate in both biological processes simultaneously. In
order to do some basal research linking developmental associate miRNA
toward morphology, the expression profiles of reported morphological
miRNAs were performed in other six tea varieties with obvious
differences in leaf morphology. The focused 21 miRNAs were again
checked by quantitative PCR for their expression levels in several
representative cultivars. Bud, the 1st leaf, the 2nd leaf, the younger
stem, and the older stem were also sampled from each variety in each
group for quantitative PCR analysis ([186]Figure 8). For each miRNA,
the relative expression in JFH was set as the reference so as to get
better understanding of the expression levels among varieties. In
general, most of the miRNAs had similar expression patterns in ZYD,
PYTZ, DYWL, and QS, with high expression levels in the 1st leaf and
then the 2nd leaf. And for other miRNAs were likely to have high
expression levels in the bud, such like miR160x in HD and miR156y in
JFH.
Figure 8.
[187]Figure 8
[188]Open in a new tab
Heat map of the relative expression levels of 21 miRNAs in different
tissues by Real-time PCR. (A) Samples from tea tree cultivars
corresponding to (B). Heat map of the relative expression levels of 21
miRNAs in different tissues from the Camellia sinensis “Jinfenghuang”
(JFH), Camellia sinensis “Pingyang Tezaocha” (PYTZ), Camellia sinensis
“Zhengdayin” (ZDY), Camellia sinensis “Dayewulong” (DYWL), Camellia
sinensis “Huangdan” (HD), Camellia sinensis “Jiukeng 6” (JK), Camellia
sinensis “Queshe” (QS).
Network Analysis on miRNA, Target mRNA, and Quality-Related Metabolites
Correlation analyses were conducted to figure out the extent of the 21
conserved miRNA mentioned above participated in the quality formation
progress ([189]Figure 9). Surprisingly, compared to catechins and
caffeine, theanine had a quite strong relationship and relatively large
numbers with target mRNAs, which could be definitely set as the center
regulated metabolite, at least during the development of PYTZ sprouts.
For the galloylated catechins, ECG was regulated by much more multiple
mRNAs than EGCG. The kinds of mRNA participated in regulating the
nongalloylated catechins were approximately equal. As for the conserved
miRNA, miR8577-x, miR160-x, and novel-m0187-3p were the most social
ones. However, the neighborhood relationships for miR156-x and miR319-y
seemed rather simple.
Figure 9.
[190]Figure 9
[191]Open in a new tab
The miRNA-mRNA-metabolites association networks. The round rectangle
placed in the center were several metabolites, the ellipses at the
interlayer were potential target mRNAs, and the diamond at the
outermost layer were miRNAs. Line thickness represented the strength of
the relationship.
Discussion
Different miRNAs Were Involved in Different Tissues and Stages During the
Sprout Development in PYTZ
How a plant builds leaves from a few cells that grow, divide, and
differentiate to form into the complex organ has been well-studied, the
same with the research of mechanical regulation (Braybrook and
Kuhlemeier, [192]2010; Chen, [193]2012; Qi et al., [194]2017). miRNA,
typically multigene families, allowing for subtlety and complexity of
control in different regulatory processes, are described as factors in
many aspects of plant development (Kidner, [195]2010).
As tea is the leaf-using crash plant, the performance of leaf
development has more practical significance. Thus, we took the
perspective of looking at the function of miRNA in different tissues
and development stages separately. In the comparison of G1 and G2
mentioned above, miRNAs that participate in leaf and stem development
are different, hinting that different miRNAs need coordinate working in
leaf and stem developmental process, separately. For example, both
novel-m0155-3p and novel-m0331-5p were stem-specific miRNAs and have
high expression levels in sS2 than in sS1, which had strong
possibilities responsible for stem elongation. And likewise, one kind
of miRNA may function differently in different tissues development. In
the 78 miRNAs that both exist in leaf and stem ([196]Figure 2B), 28 of
them have different including miR5049y and other 27 novel miRNAs could
be found both in leaf and stem, in which 23 of them share one similar
expression pattern ([197]Supplementary Figure 6) (P < 0.05, log2-fold
change ≥ 1.5), 27 of them with the expression levels in sL1 higher than
in sL2 ([198]Supplementary Figures 6A–C), 24 of them with the
expression levels in sS2 higher than in sS1 ([199]Supplementary Figures
6A,D). The non-conserved miR5049, included in profile1
([200]Supplementary Figure 6A), had been reported to be drought stress
response miRNA in the root of drought-tolerant cultivar wheat (Akdogan
et al., [201]2016).
Interestingly, the number of down-trend miRNAs (157 from profile 3,
profile 0 and profile 1) are much more than up-trend ones in G1 (53
from profile 6, profile 7 and profile 4) ([202]Figure 3A) and in G2 (91
vs. 64) ([203]Figure 3B), this is the same case for the abundance of
these miRNAs. Discarding the ones with TPM lower than 100, the
abundance of down-trend miRNAs occupied 79.04% and 54.30% in G1 and G2,
separately ([204]Supplementary Table 10). The percentage indicated that
during the development, especially in leaf, mRNAs regulated by miRNA
have a large percentage in increasing tendency. This trend is
consistent with the expression levels from tea leaf transcriptome that
72% of genes were up-regulated in the second leaf stage compared to the
first leaf stage (Guo et al., [205]2017). The coherence of expression
levels of the regulator and the content of secondary metabolites is
particularly impressive, a point we return to below.
Evolutionarily Conserved miRNAs Were Closely Connected to Morphogenesis
Functions During the Sprout Development
miRNA is usually be concerned whether to be conserved or not, which
typically depends on their degree of presentation in all or at least
most of the species, and thus the division could be influenced by
sampling and the phylogenetic diversity of available species that
miRNAs have been characterized and annotated (Baldrich et al.,
[206]2018). In the 21 developmental associated miRNAs filtered out in
this study ([207]Table 4), except for 6 novel miRNAs, 11 of the known
miRNAs are conserved miRNAs that could be found in at least 33 plant
species out of 72 plant species ([208]Table 2). Notably, none of the
miRNAs from tea could be found in Chlamydomonas reinhardtii, which is a
conventionally model for a photosynthetic cell in studying
photosynthesis (Funes et al., [209]2007), abiotic stress (Hema et al.,
[210]2007), circadian clock (Ral et al., [211]2006) and so on. This
result to some extent is an echo of the conservation of miRNA within
one kingdom, and no miRNA had been found conserved in green algae and
land plants (Baldrich et al., [212]2018). Interestingly, there was no
tea miRNA found in Populus euphratica either, which is an ideal model
system of woody plants for research into the abiotic stress resistance
(Li et al., [213]2009), such like drought (Li et al., [214]2011) and
salt (Li et al., [215]2013). Previously studies on Populus euphratica
had been reported that only 9 out of 21 miRNAs families (miR156;
miR163; miR172; miR398; miR393; miR171; miR408; miR169; miR472) were
conserved in other plants, with other 12 miRNA family candidates show
none homologies in Populus, Arabidopsis, and Oryza (Li et al.,
[216]2009) and can thus be considered as quite ancient and independent
evolution species.
The seven tea cultivars studied above were famous cultivars that
frequently used in producing fermented or non-fermented tea in China.
Typically in the tea processing industry, tea-processing suitability
and tea quality are basically determined by the main characteristic
metabolic compounds, which directly linked up with the development and
morphogenesis of the tea sprouts (Xia et al., [217]2017). In the 11
known conserved tea miRNA, miR156 was the most popular one, which was
found in 51 plant species, followed by miR396 and miR166, which were
conserved in 47 and 45 plant species, respectively. Similar with their
high abundance (miR166, miR319, miR396, miR160, and miR390 were listed
in the top five kinds of miRNAs), which means that conservation is not
only represent low sequence variation across diverse plant species, but
also to be the large and older miRNA families with abundant copy and
target number (Chavez Montes et al., [218]2014), in order to grantee
their tightly constrained roles in function and less gene loss in the
regulatory network (Shi et al., [219]2017). Mature miRNAs in plant
often have multiple target genes with similar complementary sequences,
among which these evolutionarily conserved miRNAs and their
predominantly target genes characteristically play essential roles in
developmental regulation, morphogenesis, stress responses (Axtell and
Bowman, [220]2008; Yang et al., [221]2013). Even more noteworthy is,
the tissue-specific miRNAs ([222]Table 2) may contribute to the
development of the specific tissue, which doesn't mean they are the
dominated ones and conserved ones, either. For example, miR319 had been
reported to increase the number of longitudinal small veins thus might
account for the leaf blade width (Yang et al., [223]2013) and miR159
was involved in stem elongation (Tsuji et al., [224]2006), but they are
both conserved miRNAs and have expression in bud, leaves, and stems.
Only miR164, miR393 and miR2111 were leaf-specific miRNAs and conserved
miRNAs as well. Not all of the conserved miRNAs have similar expression
patterns in the investigated cultivars ([225]Figure 8), which might be
the result of the flexible of the “fine-tuners,” to enhance the ability
of a fast response to evolution (Muleo, [226]2012). The conservation in
sequence doesn't always represent functional conservation (Ason et al.,
[227]2006). Though the evidence of miRNA in the plant is less than that
in the animal, it is widely accepted that plants miRNA genes are
evolved independently as they do in the animal kingdoms. It is thus
believed that the larger the miRNA family is, which means the more
multiple paralogous copies of one miRNA in plants, brings more flexible
of evolution rate and more possibility in diversification though, the
more essential function in development it may be involved. More
research combined the contents of metabolites with the function of
conserved miRNAs in species-level phenotypic differences needs to be
further studied.
Development Associate miRNAs Might Play Crucial Roles in the Quality
Formation Together With Their Potential Target Transcription Factor Genes
Higher plants evolved precise and robust spatio-temporal patterns of
gene regulatory systems, among which transcription factors and miRNAs
are two of the best studied regulatory mechanisms separately at
transcriptional and post-transcriptional level (Chen and Rajewsky,
[228]2007). TFs and miRNAs generally do not work isolation, but
instead, together with co-regulators in the same layer or not, to form
large networks of cooperating and interacting in complex multicellular
organisms (Dawid, [229]2006). But they are usually positioned at the
center of regulating many aspects of developmental plasticity along
with the life cycle (Rubio-Somoza and Weigel, [230]2011). We perform
network and enrichment analysis to the 352 TF genes targeted by miRNA
in STRING (Szklarczyk et al., [231]2017), and marked the developmental
relative process ([232]Supplementary Figure 7). Most of them were
involved were clustered in some certain pathway in regulating gene
expression and primary metabolic process, which could be easily
understood that TFs were involved in the primary metabolic process
because of the fundamental maintenance of living for plant themselves.
In shrinking the research objectives, we further perform GO analysis to
the DEM in G1 and G2 above, many pathways were enriched in
developmental and morphological process ([233]Supplementary Figure 8),
especially in the G1 up trend expression profiles, which also confirmed
us the effective way of filtering key miRNAs. The eight miRNA-TF mRNA
pairs verified by their reciprocal expression relationships were much
more likely to participate in development regulation, which doesn't
imply only this eight pairs of miRNA-TFs were involved.
Metabolism was along with the development. It has been reported that
plant miRNA were widely involved in quality formation regulation (Wu et
al., [234]2014; Liu et al., [235]2017b). As the three kinds of
characteristic metabolites which finally determine the quality of tea
(Xia et al., [236]2017; Wei et al., [237]2018), catechins mainly confer
astringent taste, theanine contributes to the umami and sweet tastes,
and caffeine offers a bitter taste (Wei et al., [238]2018). We did GO
pathway enrichment analysis to the target mRNA genes of DEM between
every two samples from Bud, sL1, sL2, sS1, and sS2, with the evolved
participants showed in [239]Table 5. The correlation of the chemical
analysis on catechins, theanine, caffeine and the soluble matter would
finally affect the sensory evaluation of green tea taste. In our study,
we found that theanine turned out to be even more active in the network
([240]Figure 9). Target mRNAs which belonging to TF genes were further
picked out and constructed a more metabolic directivity one due to
their correlation ([241]Supplementary Figure 9). TF genes like
CSA013022 (HD ZIP) and CSA029222 (ARF) had a strong positive
relationship with the biosynthesis of theanine, while the later also
positively regulated EC, referring to miR165-y, miR166-z, and miR160-x,
respectively. CSA031667, a TF gene belonging to SBP family, had a
positive correlation with C, and be controlled by miR156-x. When taking
the eight miRNA-TF pairs mentioned above into consideration together,
there were at least two triplets that participate in both development
and quality formation: miR156-x-CSA031667 (SBP)-C and
miR319-y-CSA036087 (MYB)-theanine. Molecular mechanisms of sprout
development and accumulation of metabolites would be gradually
uncovered after the release of tea tree genome (Xia et al., [242]2017;
Wei et al., [243]2018) and tea organic transcriptomes (Zheng et al.,
[244]2015, [245]2016; Guo et al., [246]2017; Liu et al., [247]2017a).
More connections would be further studied toward small RNAs to improve
breeding efficiency of developing better cultivars with higher quality.
Table 5.
GO pathway enrichment analysis to the target mRNA genes of DEM between
each two samples from Bud, sL1, sL2, sS1, and sS2.
Characteristic pathway Metabolic pathway Annotated genes Related mRNA
miRNA
Flavonoid pathway Phenylalanine metabolism Amidase CSA014307 miR8577-x
Omega-amidase CSA016772 novel-m0263-3p
Phenylpropanoid biosynthesis Cinnamyl-alcohol dehydrogenase (CAD)
CSA019604 novel-m0018-3p, novel-m0119-3p, novel-m0279-3p,
novel-m0498-3p, novel-m0608-3p, novel-m0674-3p, novel-m0888-3p,
novel-m0904-5p, novel-m0977-3p, novel-m1032-3p, novel-m1080-3p,
novel-m1106-3p
CSA005610 miR1511-y, miR5139-x, miR8155-y, novel-m0999-3p
CSA009196 miR6118-y
Peroxidase CSA026001 miR5658-x, novel-m0675-3p
Beta-glucosidase CSA014637 novel-m0349-3p, novel-m0360-3p,
novel-m0428-3p
Shikimate O-hydroxycinnamoyltransferase (SHT) CSA011696 miR8577-x
Flavonoid biosynthesis Naringenin 3-dioxygenase (F3H) CSA004930
miR8577-x
Chalcone isomerase (CHI) CSA006623 novel-m0466-3p, novel-m0809-3p
Anthocyanin biosynthesis Anthocyanidin 3-O-glucoside
5-O-glucosyltransferase (UGT75C1) CSA036671 miR1865-x, novel-m0674-3p,
novel-m1032-3p
UGAT CSA013643 miR4995-x, miR6483-y
Caffeine Pathway Caffeine Metabolism Urate oxidase (UOX) CSA020658
novel-m0021-5p, novel-m0508-3p, novel-m0698-3p, novel-m0834-3p,
novel-m0964-3p, novel-m0964-3p
Xanthine dehydrogenase/oxidase CSA006612 miR5059-x
Amino acid Pathway Citrate cycle (TCA cycle) Isocitrate dehydrogenase
(NAD+) CSA018911 miR5083-y
Pyruvate dehydrogenase E2 component (dihydrolipoamide
acetyltransferase) CSA034852 miR8175-y
2-oxoglutarate dehydrogenase E1 component CSA032631 novel-m0381-3p
2-oxoglutarate dehydrogenase E2 component (dihydrolipoamide
succinyltransferase) CSA019997 novel-m0184-5p
ATP citrate (pro-S)-lyase CSA033096 CSA027464 CSA034373 miR5059-x
Succinyl-CoA synthetase beta subunit CSA027414 novel-m0018-3p,
novel-m0119-3p, novel-m0279-3p, novel-m0498-3p, novel-m0608-3p,
novel-m0674-3p, novel-m0888-3p, novel-m0904-5p, novel-m0977-3p,
novel-m1032-3p, novel-m1080-3p, novel-m1106-3p
Aconitate hydratase CSA024268 CSA025056 CSA020729 novel-m0945-3p
Fumarate hydratase, class II CSA004158 novel-m0574-5p
CSA036624 miR5054-y
Succinyl-CoA synthetase beta subunit CSA027414 novel-m0018-3p,
novel-m0119-3p, novel-m0279-3p, novel-m0498-3p, novel-m0608-3p,
novel-m0674-3p, novel-m0888-3p, novel-m0904-5p, novel-m0977-3p,
novel-m1032-3p, novel-m1080-3p, novel-m1106-3p
Alanine, aspartate and glutamate metabolism Glutamate dehydrogenase
(NAD(P)+) CSA030852 CSA002178 CSA019055 novel-m0541-3p
Alanine transaminase CSA034511 novel-m0410-3p, novel-m0435-3p,
novel-m0549-3p, novel-m0553-5p, novel-m0554-5p, novel-m0752-5p,
novel-m0769-3p, novel-m0863-3p, novel-m0938-3p, novel-m1097-3p,
novel-m1123-5p
4-aminobutyrate—pyruvate transaminase CSA023854 novel-m0074-3p,
novel-m0075-3p
Calcium-binding protein CML CSA036463 novel-m1108-5p
CSA009111
Amidophosphoribosyl transferase CSA032114 miR5059-x, miR6281-x
3-deoxy-7-phosphoheptulonate synthase CSA021087 novel-m1042-3p
CSA011007 novel-m0657-3p, novel-m1042-3p
Arogenate/prephenate dehydratase CSA006958 miR5054-y, novel-m0388-3p
[248]Open in a new tab
Data Availability
The datasets generated for this study can be found in Sprouts
development of tea plants, PRJNA510482.
Ethics Statement
The authors declare that we have complied with all relevant ethical
regulations.
Author Contributions
ZD conceived the study. LZ, CC, JS, and YW performed the experiment and
analyzed the data. LZ wrote the paper. All authors read and approved
the final manuscript.
Conflict of Interest Statement
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.
Acknowledgments