Landscape of RNAs in human lumbar disc degeneration

Accumulating evidence indicates noncoding RNAs (ncRNAs) fine-tune gene expression with mysterious machinery. We conducted a combination of mRNA, miRNA, circRNA, LncRNA microarray analyses on 10 adults' lumbar discs. Moreover, we performed additional global exploration on RNA interacting machinery in terms of in silico computational pipeline. Here we show the landscape of RNAs in human lumbar discs. In general, the RNA-abundant landscape comprises 14,635 mRNAs (37.93%), 2,059 miRNAs (5.34%), 18,995 LncRNAs (49.23%) and 2,894 (7.5%) circRNAs. Chromosome 1 contributes for RNA transcription at most (10%). Bi-directional transcription contributes evenly for RNA biogenesis, in terms of 5′ to 3′ and 3′ to 5′. Despite the majority of circRNAs are exonic, antisense (1.49%), intergenic (0.035%), intragenic (1.69%), and intronic (6.29%) circRNAs should not be ignored. A single miRNA could interact with a multitude of circRNAs. Notably, CDR1as or ciRS-7 harbors 66 consecutive binding sites for miR-7-5p (previous miR-7), evidencing our pipeline. The majority of binding sites are perfect-matched (78.95%). Collectively, global landscape of RNAs sheds novel insights on RNA interacting mechanisms in human intervertebral disc degeneration.


IntroductIon
Intervertebral disc degeneration (IDD) is the chief contributing factor to low back pain with deleterious medical and social impact [1].Increasing evidence reveals that IDD is a multifaceted spinal disease.Both genetic and environmental factors contribute to IDD, the chief of which are genetics [2,3].The genetics machinery ranges from single-nucleotide variants [4] and coding genes [5][6][7][8], to newly defined noncoding RNAs (ncRNAs).
Recent whole-genome sequencing (WGS) studies have greatly expanded our scopes upon human genome [9,10].As a multitude set of transcript products of the human genome, ncRNAs account for 98% of the human genome devoid of protein-coding function.

Research Paper
exist in living things.Moreover, ncRNAs play critical roles in a variety of biological processes pertaining to gene expression [11,12].
The ncRNA superfamily comprises small ncRNAs (miRNAs), long ncRNAs (LncRNAs) and the newly defined circular RNAs (circRNAs) in terms of length and structure.Indeed, RNAs play fundamental roles in organisms.Not only different lengths, but sequences or structures (e.g., riboswitches) of RNAs, are critical [13,14].In structure, miRNAs are ~21-nucleotidelong noncoding RNAs.LncRNAs contain nucleotides over 200, even more than 100,000 [15].CircRNAs are a class of ncRNAs with covalent linked ends.Emerging as the first representative of ncRNA family, miRNAs direct an effector protein Argonaute (AGO) to suppress their targets' expression [16].The mysterious class of LncRNAs has been addressed recently with increasing uncovered numbers.Furthermore, a growing body of evidence indicates that LncRNAs play important roles in each issue of human biology [17].LncRNAs can be divided into different subgroups in terms of their genomic contexts: stand-alone LncRNAs or lincRNAs [18]; Natural antisense transcripts [19]; Long intronic ncRNAs [20]; Divergent transcripts, promoterassociated transcripts and enhancer RNAs [21].In mechanism, LncRNAs are key regulators in epigenetics as recruiters, tethers and scaffolds [22]; in transcription as decoys, coregulators and Pol II inhibitors [23]; in posttranscriptional regulation involving in mRNA processing, stability and translation [23].
In contrast, studies on circRNAs are at their early stage.It has been noted that circRNAs act as posttranscriptional regulators.CircRNA expression is tissue/ developmental stage specific.They can interact with miRNAs via miRNA sponges and competing endogenous RNA (ceRNA) in the cytoplasm [24][25][26].In details, miRNA sponges indicate that a circRNA with miRNA binding sites could absorb the miRNA and eliminate the original repression on the miRNA-targeted gene.
Regarding ncRNAs in human IDD, we have addressed the expression profiles of miRNAs previously using scoliotic nucleus pulposus (NP) as control [27].Moreover, we addressed the expression profiling of LncRNAs and mRNAs in IDD using normal NP control [28].It has been noted that there might be a regulatory link between LncRNAs, circRNAs and miRNAs.LncRNAs and circRNAs could repress miRNAs via various mechanisms, e.g.miRNA sequestration.Furthermore, the interaction of circRNAs and LncRNAs remains mysterious.
On the other hand, it should be stressed that discs from scoliotic patient are abnormal [29].Accordingly, we conducted miRNAs and circRNAs microarray analyses using the same total RNA samples from lumbar NP tissues as LncRNA-mRNA microarray.Moreover, advanced bioinformatics data analyses were employed.

results superseries of human rnAs in lumbar disc samples
Based on tetrad array platforms using the same 10 lumbar disc samples, we successfully established a SuperSeries in GEO (GSE67567) [28,30] (Figure 1A).The BioProject accession of this SuperSeries is PRJNA280271, which encompasses and links the expression profiles of miRNAs (GSE63492), LncRNAs and mRNAs (GSE56081) and circRNAs (GSE67566).We addressed RNAs expression profiles by comparing RNAs in 5 degenerate disc samples with those in healthy disc samples.

Global analysis reveals rnA transcription and binding truth
Functional annotation generated 3,146 potentially interactive networks of differentially expressed circRNAs and miRNAs in terms of computational pipeline (Supplementary PDF S1).
Chromosome 1 transcribes RNAs at the highest level (10%); whereas Chromosome Y contributes the least.For circRNA transcription, Chromosome 17 ranks as the second candidate, following by Chromosome 2. For mRNA transcription, Chromosome 19 ranks as the second one, following by Chromosome 11.For LncRNAs, the ranking order is Chromosome 2 and 3 (Figure 2A, 2B, normalized value = 1).Bi-direction transcription contributes nearly equally for RNA biogenesis in terms of 5ʹ to 3ʹ and 3ʹ to 5ʹ (Figure 2C).

Interacting hallmarks from the point of circrnAs
We noted that circRNAs commonly combined with the 3p or 5p of miRNA (Supplementary Table S9).It is well established that miR-3p or miR-5p derives from the arms of stem loops of pre-miRNAs.We termed the phenomena as Stem Loop Arm Sponge (SLAS, Figure 4A).Moreover, we found that a single circRNA commonly interact with at most 5 miRNAs (99.31%) as sponges following intensive exploration.In total, there were 20 circRNAs interacting with less than 5 miRNAs (range 1-4, Supplementary Table S10).Considering the available evidence [26] and the principle of complementary base pairing, we putatively put forward the phenomena as sponge saturation (Figure 4B, Supplementary Figure S1).

binding principle for ncrnAs and mrnAs
It is well established that RNAs have polarity with one 3ʹ ends and one 5ʹ ends.At the 3ʹ ends, most eukaryotic mRNAs have a sequence of polyadenylic acid, referred to as the poly (A) tail.Appropriate poly (A) tail length is crucial for efficient translation, being ~70 residues [35,36].As well, the methyl guanosine 'cap' at the 5ʹ ends of all eukaryotic mRNAs plays critical roles during mRNA processing and metabolism [37].
Pre-RNAs undergo several steps of RNA splicing during which the phosphodiester bonds at exon-intron boundaries are cleaved and the introns are excised.The process is transcription during which RNAs mature.circRNAs are chiefly the products of exons splicing cotranscriptionally mediated by the spliceosome and flanking introns [24].
Furthermore, one RNA binds with another RNA in accordance with the polarity principle.In details, when a single miRNA meets its targeted mRNA, the 3ʹ end of the miRNA binds with the 5ʹ end of the mRNA; whereas the 5ʹ end of the miRNA binds to the 3ʹ end of the mRNA.
Traditionally, miRNAs are reported to repress gene expression by binding with the seed regions of the 3ʹUTRs of their targeted genes.Accumulating evidence suggests that miRNAs can suppress gene expression by complementary interactions with the coding sequence (CDS) and 5ʹUTRs of mRNAs [38][39][40].
Based on previous lines of evidence and our findings, back-spliced exons and introns produce circRNAs co-transcriptionally. Subsequently, circRNAs accumulate in the most crucial regions for mRNAs processing and metabolism, i.e., cap and poly (A) tail regions, acting as miRNA sponges via the SLAS machinery by interacting with miR-3p or miR-5p.Consequently, circRNAs frequently get saturation by interacting with at most 5 miRNAs.These fine tuning processes reflect the regulatory networks between circRNAs, miRNAs and mRNAs (Figure 4C).

landscape of rnA length
We analyzed the length of each type of RNAs.Consequently, we found that the length for mRNA, circRNA and LncRNA varies greatly.For mRNAs, the median length is 2,446 nt, ranging from 80 to 43,816 nt.For circRNAs, the median length is 3,648 nt, ranging from 74 to 433,729 nt.For LncRNAs, the median length is 833 nt, ranging from 61 to 106,351 nt (Figure 4D, Supplementary Table S11, S12 and S13).

mirnA microarray data
Based on miRNA microarray (GSE63492), we found there were 50 deregulated miRNAs passing student t test with a P value < 0.05.There were 31 up-regulated and 19 down-regulated miRNAs.The top 10 up-regulated and top 10 down-regulated miRNAs were listed in Supplementary Text S1.
We addressed the expression profiles of miRNAs in IDD using scoliotic NP as control [27] (GSE 19943), amongst which 10 miRNAs were upregulated and 67 miRNAs were downregulated.In addition to the difference of microarray platform (miRCURY™ LNA Array v.11.0 vs. miRCURY™ LNA Array v.18.0), we found that the expression profiles of miRNAs were largely different in NP samples in terms of normal, scoliotic and IDD.Collectively, these findings indicate that miRNAs play a role not only in the etiology of IDD, but in the underlying mechanisms of scoliosis.

CircRNA expression profiles in IDD
Using the Arraystar Human Array analysis, we identified the expression profiles of 2,894 circRNAs (GSE67566).Amongst the 636 differentially expressed circRNAs, there were 354 up-regulated and 282 downregulated.The top 10 up-regulated and top 10 downregulated circRNAs were listed in Supplementary Text S2, which includes annotations as pertaining genes.

core interacting mirnAs and circrnAs in human Idd
Given that one miRNA could interact with a number of circRNAs and vice versa, we tried to narrow the scope and find valuable clues.We screened top 20 deregulated miRNAs in Supplementary Text S1 and top 20 deregulated circRNAs in Supplementary Text S2.Consequently, we found top 7 deregulated miRNAs may interact with 32 circRNAs via 84 potential binding sites; whereas 20 deregulated circRNAs may interact with 99 miRNAs via 143 potential binding sites.Notably, miR-328-5p could interact with 12 circRNAs via 49 binding sites; whereas circRNAs generally interact with 5 miRNAs (Supplementary PDF S2).

dIscussIon
The study is the first addressing mRNAs and ncRNAs in a type of human tissues rather than cultured cells.Furthermore, the expression profiling of RNAs sheds novel light on the understandings of not only the broadly influencing disease as lumbar IDD, but RNAs per se.These RNA hallmarks might reflect the true scenarios in human lumbar disc diseases and low back pain.Moreover, the novel vision might open a new page for RNA studies, based on RNA derivation from different Chromosomes.Strikingly, the length of mRNAs, circRNAs and ncRNAs varies to a great extent, the shortest of which is 2 lncRNAs (mascRNA and lincRNA-NFIA-2) with 61 nt.In contrast, the longest RNA is a circRNA (circRNA-100723) with 433,729 nt.These extreme RNAs might be the key to the treasure door of the truth on RNAs.
Traditionally, it has been reported that circRNAs derive from back-spliced exons [24,25,31,32].Surprisingly, we found nearly 10% circRNAs belong to intronic, antisense, intergenic rather than exonic types.Differentially expressed circRNAs remain the trend as shown in Figure 5B.The finding greatly expands our vision on circRNAs and splicing.
Accumulating evidence indicates that circRNAs expression has spatio-temporal features.Memczak et al. [32] identified 1,903 circRNAs in mouse tissues (brains, fetal head, differentiation-induced embryonic stem cells); 81 of which mapped to human circRNAs.Indeed, the brain is circRNAs abundant in human, mouse [31] and porcine [43].The abundance of circRNAs is spatial and temporal depending, with the maximum in cortex at day 60 of gestation with 4,634 circRNAs [43].As for H9 human embryonic stem cells, the abundance of circRNAs is 1,662 [26].
Bachmayr-Heyda et al. [44] compared the expression abundance of circRNAs in 13 human tissues.The most abundant tissue expressing circRNAs is brain and the lowest is muscle.In testis, Sry acts as miR-138 sponge [32].We identified 2,894 circRNAs in human lumbar discs as a novel circRNA-abundant tissue.
Given that human brain, testis, stem cells and discs belong to immune privileged sites with FasL/Fas expression [45,46], we putatively propose that circRNA expression might be tightly linked to immune privileged organs and tissues.

MAterIAls And Methods human nP sample collection
We collected human normal (cadaveric donors) and degenerative lumbar nucleus pulposus (NP) with patient demographics and IDD grading as we previously described [47].Human disc bank was established with normal and degenerative samples and demographic records (Supplementary Figure S2).

rnA isolation and quality control assay
Total RNA isolation and quality control assay were conducted as we previously describedas TRIspin method [47,48].

circrnA array analysis
The sample preparation and microarray hybridization were performed based on the Arraystar's standard protocols.Labeled cRNAs were hybridized onto the Arraystar Human circRNA Array (8 × 15 K, Arraystar).Subsequently, Agilent Feature Extraction software (version 11.0.1.1)was used to analyze acquired array images.Quantile normalization and subsequent data processing were performed using the R software package.Hierarchical Clustering was performed to show the distinguishable circRNAs expression pattern among samples.

lncrnA and mrnA microarray analysis
LncRNA and mRNA microarray analysis was conducted as we previously reported [47].

Annotation for circrnA/mirnA interaction
Accumulating evidence indicates circRNAs play a crucial role in the delicate tuning of miRNAmediated regulation of gene expression by sequestering corresponding miRNAs.The interacting machinery between circRNAs and miRNAs associated with diseases suggests that circular RNAs are important for disease regulation [49].The circular RNA ciRS-7 contains various, tandem miRNA-7 binding sites, thereby acting as an endogenous miRNA "sponge" to adsorb, and hence quench, normal miRNA-7 functions [50].ciRS-7 might serve as a crucial factor engaged in the functioning of neurons as well as a candidate in neurological disorders and tumor development on basis of several lines of evidence as follows: the widespread involvement of miR-7 as a key regulator of various cancer pathways; the suggested implications of miR-7 in Parkinson's disease by direct targeting of a-synuclein protein expression; direct targeting of the ubiquitin protein ligaseA (UBE2A) in Alzheimer's disease.

Quantitative real-time Pcr (qrt-Pcr) assay
Candidate target genes were validated with highly reliable biotechniques such as qRT-PCR [53].The primers were listed as Supplementary Text S3.

statistical analysis
Statistical analyses were performed using SPSS 19.0 software package.To analyze the expression difference of the particular lncRNAs or mRNAs in microarray and PCR analysis, student's t-test was applied to the study.Furthermore, the Benjamini-Hochberg FDR (the FDR cutoff was 0.05) was used for multiple-testing correction.The threshold of significance was set as P-value < 0.05.

Figure 1 :
Figure 1: superseries of ncrnAs in human lumbar discs.(A) Schematic diagram of human lumbar disc RNA SuperSeries based on tetrad platforms.(b-e) delineate the Hierarchical clustering of mRNAs, miRNAs, circRNAs and lncRNAs; whereas (F-I) represent the boxplots of mRNAs, miRNAs, circRNAs and lncRNAs.(J) indicates the ratio of RNA subgroups in human lumbar discs.(K) indicates the ratio of differentially expressed RNA subgroups in human IDD.

Figure 2 :
Figure 2: the landscape of rnAs and their binding sites.(A) Represents the distribution diagram of mRNAs, circRNAs and LncRNAs in each Chromosome in terms of each type of RNA; whereas (b) represents the constituent ratio of mRNAs, circRNAs and LncRNAs in each Chromosome in terms of Chromosome (normalized value = 1.0).(c) delineates the constituent ratio of transcription direction of RNAs.(d) indicates the constituent ratio of miRNAs interacting with circRNAs.(e) represents the hallmarks of binding sites in terms of miRNA types; whereas (F) represents the features of binding sites in terms of global, top 10 differentially expressed and exceptional views.

Figure 3 :
Figure 3: circrnA families interacting with mirnAs.(A-d) delineate circRNA families binding with corresponding miRNAs.(e and F) indicate exceptional cases of multiple binding sites between circRNAs and miRNAs.

Figure 4 :
Figure 4: the mirnA-circrnA binding machinery and landscape of rnA length.(A-b) indicate miRNA-circRNA interacting machinery.(c) is the diagram of RNA transcription and splicing.(d) delineates the length scope of each type of RNAs.

Figure 5 :
Figure 5: circrnAs and mirnAs in human intervertebral disc degeneration. (A-b) Represent the constituent ratio of circRNAs and differentially expressed circRNAs.(c) indicates the core interacting miRNA, circRNA, and mRNA in human lumbar discs.