Root phenotype under conventional and root restricted cultivation
To observe the root phenotype differences under conventional cultivation (referred to NR) and root restricted cultivation (or root restriction, referred to RR) (Figure 1 and Figure 2). One-year-old self-rooted grapevine cv. ‘Muscat Hamburg’ was planted by two different ways, respectively. The entire root system was sampled and photographed to record root phenotype during different developmental stages. A total of twelve sampling time points were conducted (at one development stage, root systems were simultaneously collected from two different cultivation methods). The phenotype of root samples from twelve different times of NR and RR were marked as NR1-12 (Figure 1) and RR1-12 (Figure 2), respectively. The two cultivations results showed that similar root formation orders including new adventitious roots, absorbing roots and secondary lateral roots. Finally, the old roots degenerated, and young roots developed into the main root system. However, from the the seventh sample (NR7 and RR7), the root morphology was significantly changed between NR and RR cultivation during the developmental process. Before the seventh sampling, new adventitious roots were generated and elongated continuously under NR and RR. After the seventh sampling, the old roots gradually aged and disappeared from bottom to top. The adventitious roots of NR became coarse and gradually lignified to brown, while the adventitious roots under RR further elongated to brown, differentiated more secondary lateral roots and produced a great deal of adventitious roots continuously.
In order to further detail the root architecture variation of cv. ‘Muscat Hamburg’ in response to root restriction. The weight of young roots, old roots and total roots were further counted (Figure 3). The weight of all three kinds of roots increased during developmental process. In the later stage of development, the total weight of roots under RR8-12 was significantly higher than that under NR8-12. However, there was no obvious difference in the weight of old roots during the whole development stage. Consistent with the total weight of root system, the weight of young roots, largely derived from the lateral root morphogenesis, significantly increased in the later stage of development under RR8-12. It shows that the difference in total roots weight were caused by enhancing the induction of young lateral roots, which were in keeping with phenotype of root system (Figure 1 and Figure 2). Combining the above results, The results shows that the root weight increasing and architecture remodeling of grapevine, to a great extent, contributes to the lateral root morphogenesis on the condition of RR.
Sequencing statistics in different roots samples
To systematically investigate the potential regulatory roles of genes and ncRNAs involved in root development in cv. ‘Muscat Hamburg’, we first collected the roots derived from two different developmental atages (NR7 and NR12) in NR and roots derived from the stage with abundant lateral root in RR (RR12) and performed integrative transcriptomic analyses. As a result, 59.62 G CleanReads from a total of 9 samples (NR7_1, NR7_2, NR7_3, NR12_1, NR12_2, NR12_3, RR12_1, RR12_2 and RR12_3) were obtained. The number of CleanBases ranged from 5.35 to 8.69 G in each samples, the Q30 base distribution ranged from 93.56% to 94.28%, and the average GC content was 47.07% (Table 1). The statistical assessment confirmed that the datasets were highly reliable.
Based on the expression pattern of mRNAs, lncRNAs and circRNAs, correlation and principal component analysis (PCA) of the 9 samples were conducted. The results indicated that the three independent biological triplicates were reliably repeatability (Figure 3). Meanwhile, The 9 samples were divided into three clusters (NR7, NR12 and RR12), indicating that the closer the clustering distance or PCA distance of the samples, the more similar the samples were (Figure 4).
Identification and characterization of lncRNA and circRNAs
There were different distributions of lncRNAs and circRNAs on each chromosome, circRNA distributed on chromosome 18 was the most abundant and lncRNA located on chromosome 8 was the most (Figure 5a). The length of lncRNAs and circRNAs in grapevine were different, and the distribution of most lncRNAs and circRNAs were regular (Figure 5b). The length of 612 circRNAs (23.40 %) and 149 lncRNAs (7.56 %) were greater than 2000 nt, and 350 circRNAs were less than or equal to 200 nt. And the number of lncRNAs and circRNAs in the range of 201-2000 nt decreased with their sequence length (Figure 5b). Based on the position relationship between lncRNA and known protein coding transcripts, FEELnc software was used to count lncRNA types according to three levels, including direction, type and location. There are 657 antisense lncRNAs and 1308 sense lncRNAs were identified (Figure 5c). Bedtools software was used to classify the position relationship between circRNA and the coding genes, and five types of circRNAs were found. The largest number of circRNA of the sense-overlapping type containing 1878 circRNAs (72%). Intronic circRNAs was the lowest proportion containing 35 circRNAs (1.34%). Moreover, the percentage of exonic-derived circRNAs, antisense-derived circRNAs and intergenic-derived circRNAs accounted for 15.64% (409 circRNAs), 7.88% (206 circRNAs), and 3.33% (87 circRNAs), respectively (Figure 5d).
Differential expression analyses of mRNA, lncRNAs and circRNAs
There were three sample groups (NR12_vs_NR7, RR12_vs_NR7, RR12_vs_NR12) were used to analyze the differential expression of mRNA, lncRNAs and circRNAs. A total of 26588 mRNAs were detected in all groups. The numbers of differentially expressed genes (DEGs) were 2320, 1864 and 2440, respectively. For ncRNAs prediction, 2615 circRNAs were also characterized. However, a few circRNAs were differentially accumulated in all sample groups, which were 16, 17 and 9, respectively. Additionally, 1971 lncRNAs were identified in all sample pools, and 176, 173 and 137 were differential expression, respectively (Table 2).
We found 170 mRNAs, 5 lncRNAs, and 0 circRNA that were shared by three various groups (Figure 6a-c). There were 2320 mRNAs (1329 up-regulated and 991 down-regulated), 176 lncRNAs (110 up-regulated and 66 down-regulated), and 16 circRNAs (9 up-regulated and 7 down-regulated) that were differentially expressed in NR12_vs_NR7. 2440 mRNAs (1237 up-regulated and 1203 down-regulated), 173 lncRNAs (84 up- and 89 down-regulated), and 17 circRNAs (9 up-regulated and 8 down-regulated) that were differentially expressed in RR12_vs_NR12. Moreover, there were 1864 mRNAs (1222 up-regulated and 642 down-regulated), 137 lncRNAs (79 up-regulated and 58 down-regulated), and 9 circRNAs (including 7 up-regulated and 2 down-regulated) that were differentially expressed in RR12_vs_NR7 (Figure 6d-f) .
Functional enrichment analyses
To explore the functions of the mRNAs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed. The GO terms for DEGs mainly related to biological regulation, metabolic process, cellular process, cell, cell part, binding, catalytic activity and transporter activity in three differential groups (NR12_vs_NR7, RR12_vs_NR7, RR12_vs_NR12). KEGG pathway analysis was also performed to further explore the functions of mRNAs. The results showed that biological functions of mRNAs were significantly enriched in three differential groups, including lipid metabolism, carbohydrate metabolism, biosynthesis of other secondary metabolites and signal transduction (Figure 7).
To investigate the potential function implication of lncRNAs, GO and KEGG enrichment analyses of lncRNAs and their total differentially expressed lncRNAs were also performed. The results showed that these lncRNAs were enriched in the cellular process, metabolic process, cell, cell part, binding and catalytic activity in three differential groups. The KEGG enrichment analysis showed that these lncRNAs were clustered into nucleotide metabolism, carbohydrate metabolism and signal transduction in three differential groups (Figure 8).
To further understand the potential functions of circRNAs, GO and KEGG analyses of the host genes of circRNAs were performed. Base on GO annotation, most the host genes of circRNAs in the root were annotated to the cellular process, metabolic process, cell, cell part, binding and catalytic activity in three differential groups. The parent genes of circRNAs were involved in 5 KEGG pathways and significantly enriched in metabolism of other amino acids, lipid metabolism, global and overview maps, energy metabolism and amino acid metabolism (Figure 9).
Experimental validation of the circRNA candidates
In the present study, we experimentally tested and verified the circRNA predictions in the grapevine. Since the majority of circRNA were derived from back-splicing of mRNA transcripts, convergent primers and divergent primers were designed and used to amplify the linear transcript and backsplice sites using cDNA and genomic DNA (gDNA) as template, respectively. Unlike convergent primers, divergent primers could amplified the back-splicing junctions in cDNAs synthesised by the random primers but did not in gDNAs. Additionally, We successfully validated the circRNA - circrRNA_2377 associated with root development. The PCR amplification products were further analyzed by agarose gel electrophoresis and Sanger sequencing to confirm the occurrence of back-splicing (Figure 10).
Quantitative real-time PCR (qRT-PCR)
To verify the reliability of the sequencing data profiles, qRT-PCR was used to analyze the candidate genes related to root development (Figure 11). The expression pattern were generally consistent with the results of high-throughput sequencing, indicating that the research method used in this experiment is reliable. However, circRNA_2377 in NR12_vs_NR7 was not consistent with the high-throughput sequencing, which might be due to the expression level of circRNA itself is low (Figure 11).
CeRNA network analyses
It was supposed that lncRNA and circRNA generally contain one or multiple miRNA binding sites, and the method of miRNA target gene prediction can be used to identify lncRNA and circRNA which bound to miRNA. The functions of lncRNA and circRNA were clarified according to the functional annotations of miRNA target genes. In the present study, all the candidate ceRNA pairs (circRNA sequesters miRNA and miRNA targets mRNA) were first identified. Top 5 miRNA-lncRNA-circRNA interaction combinations were extracted based on minimum free energy according to website prediction. Furthermore, 3 miRNA-mRNA network with a high degree of confidence were also predicted and extracted. Finally, the lncRNA/circRNA-miRNA-mRNA target interaction network diagram was obtained and plotted by Cytoscape software, and constructed a hub prognostic ceRNA network consisting of 21 mRNAs, 14 miRNAs, 2 lncRNAs, and 1 circRNA (Figure 12).
As showed in Figure 12, a large proportion of mRNAs communicated with individual lncRNA and circRNA. LncRNA and circRNA acted as ceRNAs to communicate with multiple mRNAs by competing for specific miRNAs. These results suggested that the expression of lncRNA and circRNA potentially regulated gene expression through miRNA-mediated lncRNA/circRNA-mRNA ceRNA interactions, implying that ceRNA might important for lncRNA and circRNA function during root development.