MicroRNAs and their targeted genes associated with phase changes of stem explants during tissue culture of tea plant

Elucidation of the molecular mechanism related to the dedifferentiation and redifferentiation during tissue culture will be useful for optimizing regeneration system of tea plant. In this study, an integrated sRNAome and transcriptome analyses were carried out during phase changes of the stem explant culture. Among 198 miRNAs and 8001 predicted target genes, 178 differentially expressed miRNAs and 4264 potential targets were screened out from explants, primary calli, as well as regenerated roots and shoots. According to KEGG analysis of the potential targets, pathway of “aminoacyl-tRNA biosynthesis”, “proteasome” and “glutathione metabolism” was of great significance during the dedifferentiation, and pathway of “porphyrin and chlorophyll metabolism”, “mRNA surveillance pathway”, “nucleotide excision repair” was indispensable for redifferentiation of the calli. Expression pattern of 12 miRNAs, including csn-micR390e, csn-miR156b-5p, csn-miR157d-5p, csn-miR156, csn-miR166a-3p, csn-miR166e, csn-miR167d, csn-miR393c-3p, csn-miR394, csn-miR396a-3p, csn-miR396 and csn-miR396e-3p, was validated by qRT-PCR among 57 differentially expressed phase-specific miRNAs. Validation also confirmed that regulatory module of csn-miR167d/ERF3, csn-miR156/SPB1, csn-miR166a-3p/ATHB15, csn-miR396/AIP15A, csn-miR157d-5p/GST and csn-miR393c-3p/ATG18b might play important roles in regulating the phase changes during tissue culture of stem explants.

and that miR156 and miR160 could modulate the shoot regeneration 20,21 , while miR171 could influence the shoot branching 22 . In light of these evidences, it would be important to identify the miRNAs present in tea plants.
As one of the most important cash crops worldwide, tea plant (Camellia sinensis (L.) O. Kuntzes) is famous for its leaves production which is non-alcoholic and good for human health. Unfortunately, it is difficult for many researchers to study on the development mechanism of tea plant because of lacking a perfect plant regeneration system. Nowadays, although increasing attention is being paid on tea plantlets regeneration through organogenesis and somatic embryogenesis, a significant difference in regeneration frequency was observed from various explants, and very low frequency was usually witnessed during induction of many explants [23][24][25][26][27][28] . Therefore, investigation on mechanism of the dedifferentiation and redifferentiation during tissue culture might be useful for optimizing efficient regeneration system of tea plant.
Numerous conserved miRNAs and their targets have recently been identified in C. sinensis. Some of these miRNAs are responsible for cold stress 29 , Ectropis oblique feeding 30 and drought stress 31 . Using small RNA (sRNA) sequencing tech, Sun and colleagues identified 69 conserved and 47 novel miRNAs related to catechins biosynthesis from the EGCG-enriched tea plant line No. 1005 32 . It was found that 175 conserved and 83 novel miRNAs were mainly present in one bud and two tender leaves of the tea plant 33 . Although there is a large amount of information on miRNA expression related to dedifferentiation and redifferentiation in some other plants, relatively little information is available for tea plant. In the present study, Illumina HiSeq 2500 technology was occupied to sequence the putative miRNAs and mRNAs for investigating their expression profiles during the dedifferentiation and redifferentiation of tea plant tissue culture. The regulatory effect of miRNAs on the targeted genes during the phase transition of tissue culture was confirmed by qPCR.

Results
Transcriptome sequencing and assembly. Average of the obtained clean data for each sample exceeded 6GB after sequencing. A total of 194,014 transcripts was achieved from the stem explants, primary calli, redifferentiated roots and shoots through de novo assembly. 100,099 unigenes were obtained from the transcription data of these samples, with 789 bp in average length and 1,514 bp in N50 length; and 22,540 unigenes had a length of above 1000 bp, accounting for 22.52% of the total sequence number (Table 1).
High-throughput sequencing of sRNAs. After sRNA sequencing, more than 21 M raw reads and 15 M clean reads were obtained from each biological sample, around 1/10 of clean reads could mapped onto the reference transcriptome unigenes, and the mapped reads were 2.76 M and 2.62 M, 2.13 M and 1.15 M in stem explants, primary-calli, regenerated roots and shoots, respectively (Table 2). Comparison showed that a certain proportion of rRNA and tRNA, as well as a small amount of snoRNAs and Repbase repeat sequences were contaminated in clean reads, but sRNA proportion exceeded 75% except in regenerated shoots (Table 3).
Quantity proportion of the common reads and type proportion of the unique reads in the sRNA sequences were evaluated among the stem explant and its derivatives. Only a few unique reads (5.66-5.88%) shared their common types among these four samples, however, these common types represented 38.51-55.82% of the total reads ( Fig. 1), indicating that the sRNAs varied significantly during dedifferentiation and redifferentiation, and a few proportion of the sRNAs with abundant expression might play important roles during phase transition. miRnAs and their expression patterns during phase change. After blasting with the miRNA database miRBase (v21), a total of 59 conserved plant miRNAs were identified, in which 56, 54, 50 and 50 ones were observed in the explant, primary callus, regenerated root and shoot, respectively (  novel miRNAs were also predicted by the software miRDeep2 34 on Bayesian model, based on the distribution information of the reads on the precursor sequence and precursor structure energy information (RNAfold randfold), in which 133, 117, 126 and 113 ones were observed in these four types of biological samples, respectively. Thus, 198 miRNAs were identified in this study (Supplementary Table S1). These miRNAs belonged to 53 families, such as miR166, miR396, miR159, miR535, miR160, miR482, miR171_1 and so on.  Table 3. Sequence count of the small RNAs annotated in different database*. * Data in parentheses represent the percentage (%). According to calculation of TPM values of the miRNAs and the threshold of q value < 0.005 & |log2(Fold change) > 1| between different samples, 178 differentially expressed miRNAs, including 52 conserved and 126 novel miRNAs, were screened out from the detected 198 miRNAs. 94 miRNAs (24 conserved and 70 novel) were up-regulated during dedifferentiation of the stem explant. Among them, 25 miRNAs (6 conserved and 19 novel) were up-regulated furthermore while 69 miRNAs (18 conserved and 51 novel) were down-regulated during root redifferentiation; similarly, 36 miRNAs (8 conserved and 28 novel) were up-regulated while 58 miRNAs (16 conserved and 42 novel) were down-regulated during shoot redifferentiation. Meanwhile, 84 miRNAs (28 conserved and 56 novel) were down-regulated during phase change from stem explant to primary callus, of which 58 (20 conserved, 38 novel) and 52 (18 conserved, 34 novel) ones were up-regulated during regeneration of root and shoot, while 26 (8 conserved, 18 novel) and 32 (10 conserved, 22 novel) ones were down-regulated during these regenerations. Cluster analysis visually revealed that expression levels of these miRNAs were significantly different among explant, primary callus, regenerated root and shoot, especially between explant and its derivatives (Fig. 2), indicating that these differentially expressed miRNAs might be involved in regulation of phase change through its effect called posttranscriptional gene silencing. potential target genes of the miRnAs. 196 miRNAs (59 conserved and 137 novel) were predicted to be able to target 8001 potential function genes. As for differentially expressed 178 miRNAs, 1811 genes were potentially targeted by 52 conserved miRNAs, while 6025 genes were done by 126 novel miRNAs. The targeted genes were blasted with COG, GO, KEGG, KOG, Pfam, SwissProt, eggNOG, and NR databases. A total of 4264 targets were annotated, including 2964 sequences with a length greater than 1000 bp and 1300 sequences with a length between 330 bp and 1000 bp ( Table 5). Most of these targets were mainly related to the function of "transcription" (such as transcription factors and growth-regulating factors), "signal transduction mechanisms" (such as receptor protein kinase) and "posttranslational modification, protein turnover, chaperones" (Supplementary Table S2), and significantly enriched in the GO items of "growth" and "signaling" in the biological process, and "nucleoid" and "cell junction" in cellular component, as well as "electron carrier activity" and "enzyme regulator activity" in the molecular function ( Supplementary Fig. S1).  Table 4. Count result of the miRNAs and targeted mRNAs*. *There were 2 novel miRNAs without predicted target genes in 198 miRNAs. www.nature.com/scientificreports www.nature.com/scientificreports/ A total of 1518 targets were annotated in various pathways after being blasted with KEGG database (https:// www.kegg.jp/kegg/kegg1.html) 35 . Among them, 710 targets were observed during the dedifferentiation stage and enriched in KEGG pathway of "aminoacyl-tRNA biosynthesis", "proteasome", "terpenoid backbone biosynthesis", "phagosome", "cutin, suberine and wax biosynthesis", "glycerolipid metabolism", "phosphatidylinositol signaling system" and "glutathione metabolism" (Fig. 3A). During the roots regeneration, 613 targets were mainly enriched in the KEGG pathway of "porphyrin and chlorophyll metabolism", "mRNA surveillance pathway", "nucleotide excision repair", "protein export", "proteasome", "plant hormone signal transduction", "2-oxocarboxylic acid metabolism", "pyruvate metabolism" and "riboflavin metabolism" (Fig. 3B); during shoots regeneration, 624 targeted genes were mostly enriched in pathway of "folate biosynthesis", "protein processing in endoplasmic reticulum", "pyrimidine metabolism", "mRNA surveillance pathway", "nucleotide excision repair", "protein export", "2-oxocarboxylic acid metabolism", "brassinosteroid biosynthesis", "pyruvate metabolism" and "porphyrin and chlorophyll metabolism" (Fig. 3C). Obviously, genes related to the pathway of "aminoacyl-tRNA biosynthesis", "proteasome" and "glutathione metabolism" might be necessary for cell division and replication during dedifferentiation, which play an important role in the stress and defense of phase change 36,37 ; meanwhile, genes associated with pathway of "porphyrin and chlorophyll metabolism", "mRNA surveillance pathway", "nucleotide excision repair", "protein export", "2-oxocarboxylic acid metabolism" and "pyruvate metabolism" might play important roles in plastid rebuilt 37 , cellular stress 38 , DNA repair 39 and secondary metabolism during the redifferentiation.
During primary callus formation, significantly low expression of csnmiR167d was observed, while expression of its target ERF003 was specifically up-regulated. ERF003, one of the ethylene-responsive transcription factors, acts as a transcriptional activator through binding to the GCC-box promoter element and regulates expression of the genes involved in the response to stress factors and components of stress signal transduction pathways 40 . According to our study, miR167d might be an important regulator of the signaling pathway through modulating the ERF003 during callus formation of stem explant. In Arabidopsis, over-expression of miR167 would inhibit the somatic embryos formation via negatively regulating ARF6 and ARF8 13 , which implied that low expression of miR167 might be in favor of plant embryonic callus formation. Many previous studies revealed that ethylene could induce callus formation and plant regeneration in citrus 41 , apple 42 , barley 43 and Arabidopsis 44 ; and that up-regulated ethylene biosynthesis would promote the process of embryogenesis through activating the ethylene-responsive transcription factor in soybean 45 . Thus, miR167d/ERF003 module might play an important influence in dedifferentiation of stem explant.  www.nature.com/scientificreports www.nature.com/scientificreports/ As well known, the highly conserved miR156/SPL (Squamosa promoter-binding-like protein) module was reported to regulate stress response 46 , leaf and fruit development 47,48 , grain size 49,50 , shoot regeneration 21 , and callus induction 51 . In Arabidopsis, miR156 could respond to auxin signaling and modulate the target SPL 3/9/10, consequently control the quantity of lateral roots 52 . In present study, low expression of miR156 and up-regulated expression of SBP1 were observed during root and shoot regenerations. Transcriptional factor SBP1, an ortholog of SPL, can bind to the AP1 promoter and regulate the expression of SQUAMOSA involved in development 53 . Thus, the effect of miR156-SPB1 module on the redifferentiation is worthy to be further studied.
Low expression of the miR396, miR166a-3p, miR157d-5p and miR393c-3p, coupled with up-regulated expression of AIP15A, ATHB15, GST and ATG18b were observed during shoot regeneration, indicating module of miR396/AIP15A, miR166a-3p/ATHB15, miR157d-5p/GST and miR393c-3p/ATG18b might be responsible for shoot redifferentiation. In Cymbidium, promoted auxin biosynthesis was required for efficient shoot regeneration 54 . AIP15A belongs to the early auxin-responsive SAUR family gene. Up-regulated expression of this gene suggested that auxin and auxin signaling pathway were both essential for shoot regeneration of tea plant 55,56 . Previous research also showed that ATHB15 can exert its regulation at the early stage of shoot induction 57 . GST is involved in the conjugation of reduced glutathione to a wide number of exogenous and endogenous hydrophobic electrophiles and plays a detoxification role under stress conditions including pathogens infection, cold, drought and wounding treatments [58][59][60] . The extremely up-regulated expression of GST indicated that tissues and cells might have to face stress under light condition during shoot and regeneration of the callus. ATGs encode autophagy-related proteins which are required for autophagy process, and cooperate with jasmonate-and WRKY33-mediated signaling pathways in the regulation of plant defense responses to biotic and abiotic stresses 61 . Autophagy-defective mutant atg2-2 exhibited powdery mildew resistance and mildew-induced cell death in Arabidopsis 62 . In apple, overexpression of ATG18a could promote drought tolerance through activating autophagy and reducing oxidation damage 63 . Thus, up-regulated expression of GST and ATG18b reflected the stresses and defense responses took place during redifferentiation of tea plant tissue, which might be modulated by miR157d-5p and miR393c-3p.
These miRNA targeted genes were mainly involved in transcription regulation and signaling transduction. Expression change of these genes would cascadedly influence the expression of downstream genes associated with cell division, proliferation and specialization, consequently alter the metabolism pathway related to the dedifferentiation and redifferentiation, such as the genes associated with "aminoacyl-tRNA biosynthesis", "proteasome" and "glutathione metabolism", "porphyrin and chlorophyll metabolism", "mRNA surveillance pathway", "nucleotide excision repair", "protein export", "2-oxocarboxylic acid metabolism" and "pyruvate metabolism", finally lead to phase transition in tissue culture of stem explant.

conclusions
One hundred and ninety eight miRNAs were detected from the stem explant and its tissue culture derivatives, of which 178 miRNAs exhibited differential expression in various samples, 57 miRNAs exhibited phase-specific expression patterns. Expression of 12 phase-specific miRNAs and interaction of 6 miRNA/target gene modules were validated by qPCR. The miR167d, miR156, miR396, miR166a-3p, miR157d-5p and miR393c-3p might involve in regulation of dedifferentiation and redifferentiation of stem explant by targeting the ERF003, SPB1, AIP15A, ATHB15, GST and ATG18b, respectively.

Methods and Materials
preparation of tea plant samples. Tissue culture seedlings of C. sinensis cultivar 'Jinxuan' were micropropagated and maintained on the Murashige and Skoog (MS) medium with addition of 2 mg/L 6-benzylaminopurine (BAP), 0.1 mg/L naphthalene acetic acid (NAA), 30 g/L sucrose and 9 g/L agar through single-node culture. The stem explant (S_Explant) were inoculated onto the callus inducing medium (CIM) for obtaining the callus (designated as S_Primary callus), and the stem-derived callus was then inoculated onto the root inducing www.nature.com/scientificreports www.nature.com/scientificreports/ medium (RIM) in order to regenerate the roots (designated as S_Root). The detailed culture operations were described as previously paper 37 . Meanwhile, the stem-derived callus was also inoculated onto the shoot inducing medium (SIM) in order to obtain the regenerated shoot. The SIM was prepared and autoclaved at 121 °C for www.nature.com/scientificreports www.nature.com/scientificreports/ 20 min after adding 3 mg/L BAP, 0.1 mg/L NAA, 30 g/L sucrose and 9 g/L agar into MS medium. The shoots were regenerated after inoculation of the stem-derived callus on the SIM for around 60 days. Sampling was then carried out and designated as S_Shoot. The obtained samples were immediately frozen in liquid nitrogen and stored at −80 °C for further use. All the tests were conducted triply.
total RnA extraction. Total RNAs for sequencing mRNAs and miRNAs were extracted from the S-Explant, S-Primary callus, S-Root and S-Shoot samples according to the manufacturer's protocols of the RNAprep Pure Plant Kit (TIANGEN Biotech Co., Ltd., Beijing, China) and RNAiso Plus (TAKARA Bio Inc., Shiga, Japan), respectively. The purity and concentration of RNAs were checked by NanoPhotometer spectrophotometer (Implen, CA, USA) and Qubit 2.0 Flurometer (Life Technologies, CA, USA), respectively. Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was also used for assessing the RNA integrity via RNA Nano6000 Assay Kit. transcriptome analysis. The sequencing libraries were constructed by using 3 μg total RNAs for each sample according to the protocol of NEBNext Ultra RNA Library Prep Kit. The libraries were then sequenced on Illumina Hiseq 2500 platform. After removing low quality reads and the reads containing adapter and ploy-N from raw data, clean data were used to assemble transcripts through Trinity software 64 , and then function of assembled transcripts were annotated using BLAST software 65   www.nature.com/scientificreports www.nature.com/scientificreports/ were performed with primers that annealed to the ends of the adapters. Size selection of the PCR products were performed on PAGE gel, and sRNA libraries were recovered by gel-cutting and purified through AMPure XP Kit (Beckman Coulter, Australia). Finally, quality of the cDNAs libraries was ensured by examining the size, purity, and concentration through an Agilent2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The quality-checked libraries were then sequenced using a HiSeq X-Ten System (Illumina, SanDiego, CA, USA) in the Biomarker Technology Co. (Beijing, China, http://www.biomarker.com.cn). Three biological replicates at each culture stage were used for sRNA sequencing.
Small RnA assembly and analysis. The quality of original data was verified after sequencing, and clean reads were obtained by removing the low-quality reads, reads with ploy-N content greater than or equal to 10%, and sequences shorter than 18nt or longer than 30nt and cutting off the 3′ joint sequence. In order to obtain unannotated reads containing miRNAs, through BowTie software 66 , the clean reads were aligned with Silva database, GtRNAdb database, Rfam database and Repbase database respectively to filter the ncRNAs including rRNA, tRNA, snRNA, snoRNA and repetitive sequences. The unannotated reads were then mapped with the reference transcriptome library to obtain positional information by BowTie software 66 . The amounts of common sequences and the type of unique sequences were counted between different samples.
TPM value was used as expression level of miRNA, which was calculated according to equation: TPM = read count *1000000/mapped reads. DEGseq R package was used to measure the miRNA expression difference between two samples. The differentially expressed miRNAs were obtained according to the threshold of q value < 0.005 & |log2(Fold change) > 1| between different samples.
Identification of conserved and novel miRNAs. The mapped sequences were further used to align with the miRNAs from all species in the miRBase to identify the known miRNAs. miRDeep2 34 was used to identify novel miRNAs from each sample, and to carry out the structural prediction of miRNAs and their precursors.
Annotation of the miRnA targeted genes. Based on known and novel miRNAs as well as gene sequence information of corresponding transcriptome library, targeted genes were predicted by using TargetFinder software 67 . The predicted target genes were then compared with Nr, Swiss-Prot, GO, COG, KEGG, KOG and Pfam databases to obtain annotation information through BLAST software.
Verification of the miRNAs and their targeted genes by qPCR. Expression of the mature known and novel miRNAs was verified through poly (A) RT-qPCR according to the method described by Shi and his colleagues 68 . qPCRs were carried out by using TB Green Premix Ex Taq II (Clontech, USA) on an Applied Biosystems StepOnePlus Real-Time PCR System (ABI, Carlsbad, CA, USA) according to the kit protocol. The reaction was performed at 95 °C for 30 s, and 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The threshold cycle (Ct) was determined as the cycle number at which the fluorescence intensity passed a pre-determined threshold. Three biological samples for each treatment were used and all reactions were assayed in triplicate for each biological sample, and 5.8 S rRNA was selected as the reference gene 69,70 after comparison had been conducted in our pre-experiment where expression of 5.8 S rRNA was much steadier and easier to be detected than U6 in each culture stage of tea plant. The primers for the miRNA qRT-PCR were shown in Supplementary Table S3.
To validate the expression of the miRNA targeted genes, qRT-PCR was also performed on Applied Biosystems StepOnePlus Real-Time PCR System by using SYBR Premix Ex Taq II Master Mix (TaKaRa Biotechnology Co., Ltd., Dalian, China) according to the manufacturer's protocol. β-Actin 37 was used as reference gene and relative expression of the targets was calculated using equation 2 −ΔΔct . Detailed information about the primers used in this study was presented in Supplementary Table S4.
The data were statistically analyzed with SAS Version 9.0 software (SAS Institute, Cary, NC, USA) using Duncan's multiple range test at significance level of p < 0.05.