Gene processing control loops suggested by sequencing, splicing, and RNA folding

Background Small RNAs are known to regulate diverse gene expression processes including translation, transcription, and splicing. Among small RNAs, the microRNAs (miRNAs) of 17 to 27 nucleotides (nts) undergo biogeneses including primary transcription, RNA excision and folding, nuclear export, cytoplasmic processing, and then bioactivity as regulatory agents. We propose that analogous hairpins from RNA molecules that function as part of the spliceosome might also be the source of small, regulatory RNAs (somewhat smaller than miRNAs). Results Deep sequencing technology has enabled discovery of a novel 16-nt RNA sequence in total RNA from human brain that we propose is derived from RNU1, an RNA component of spliceosome assembly. Bioinformatic alignments compel inquiring whether the novel 16-nt sequence or its precursor have a regulatory function as well as determining aspects of how processing intersects with the miRNA biogenesis pathway. Specifically, our preliminary in silico investigations reveal the sequence could regulate splicing factor Arg/Ser rich 1 (SFRS1), a gene coding an essential protein component of the spliceosome. All 16-base source sequences in the UCSC Human Genome Browser are within the 14 instances of RNU1 genes listed in wgEncodeGencodeAutoV3. Furthermore, 10 of the 14 instances of the sequence are also within a common 28-nt hairpin-forming subsequence of RNU1. Conclusions An abundant 16-nt RNA sequence is sourced from a spliceosomal RNA, lies in a stem of a predicted RNA hairpin, and includes reverse complements of subsequences of the 3'UTR of a gene coding for a spliceosome protein. Thus RNU1 could function both as a component of spliceosome assembly and as inhibitor of production of the essential, spliceosome protein coded by SFRS1. Beyond this example, a general procedure is needed for systematic discovery of multiple alignments of sequencing, splicing, and RNA folding data.


Background
The numerous, very dissimilar types of bioinformatic data conspire to make integration a central problem for efficient and effective application of biological findings. Integration of data of three particular types is the goal of this paper. Gene splicing is the focus, held up as an example of how sequencing, splicing, and RNA folding data types might be used to guide research that could illuminate major mechanisms of cell biology such control of levels of ribonucleoprotein species.
Function and dysfunction of gene splicing impact embryogenesis, cell motility and viability, cell cycle arrest, and many other mechanisms of metazoan cell biology [1]. This paper stems from three remarkable observations involving splicing. The spliceosome is a large complex of protein subunits and five ribonucleoprotein subunits, the latter incorporating snRNAs. One of the snRNAs is the 164-nt RNU1. Predicted 2D molecular shapes of RNU1 include four "hairpins," conformations in which pairs of nucleic acids form a double-stranded stem while single-stranded nucleic acids form a loop. The first two of the RNU1 hairpins are already known to be bioactive through functional assays of regulation of the gene cyclin H (CCNH) [2]. The fourth hairpin, denoted herein as H, has a loop of four nts and a stem of 12 pairs of nts including eight C-G bonds (hence is very stable).
Our deep sequencing to detect small RNAs in three samples of post-mortem human prefrontal cortex produced abundant reads corresponding to a 16-nt sequence from the 3' side of the stem of H. We denote herein the 16-nt sequence as S.
Regarding small RNA context [3], Kawaji et al. engaged in unbiased exploration of 19-to 40-base sequences from small RNAs. Their pioneering report provided evidence of abundant small RNAs originating from familiar noncoding RNAs (ncRNAs) including tRNAs, snoRNAs, snRNAs, and rRNAs. Regarding tRNAs, 3' ends fragments are transported from the nucleus to accumulate in the cytoplasm, as reported by Liao et al. [4]. Bidirectional promoters suggested that small RNAs can be derived from double stranded RNAs (dsRNAs) with subsequent cleavage. Shi et al. [5] found abundant transcriptional representation of sequences immediately adjacent to-that is, offset from-predicted pre-miRNAs in the simple tunicate Ciona intestinalis (sea squirt). Langenberger et al. [6] also found transcripts offset from miRNAs in human samples, albeit at low levels unrelated to levels of the adjacent miRNAs. Taft et al. [7] first reported~18 nt RNAs in FANTOM4 data that map within -60 to +120 nt of transcription start sites of genes of humans and other metazoans. Taft et al. [8] then found miRNA-like small RNAs derived from the ends of snoRNAs in humans and other eukaryotes. Moreover, Taft et al. [9] reported 17-or 18nt RNAs with 3' ends that map precisely to the splice donor site of internal exons of mice and other metazoans. Regarding snoRNAs, Ender et al. [10] assayed human cancer cell RNAs and reported a number of human snoRNAs with miRNA-like processing signatures, evidently targeting an mRNA. Likewise, Saraiya et al. [11] used sequencing to find a 26-nt RNA from the flagellated protozoan Giardia lambia, again with miRNA-like processing and apparent RNAi activity. Other non-miRNAs of about 16 nts that are subsequences of known miRNAs have been shown by Li et al. to participate in gene regulation, targeting the 3'UTRs of target genes as efficiently as sequentially enclosing miR-NAs [12]. Importantly, Li et al. documented a long list of small RNAs, some with known sources and some not. In a generalising study, Langenberger et al. [13] discovered from sequencing data that certain small RNA subsequences of a variety of human ncRNAs are highly overrepresented in the transcriptome, extending all the above reports. They analysed low molecular weight RNAs isolated from frozen prefrontal cortex, as did we in preparation of the present report. A rapidly developing line of research on small RNAs derived from tRNAs is represented by work of Haussecker et al. [14].
Additional sources of small ncRNA are the vault RNAs,~100-nt Pol III transcripts in the enigmatic vault organelles of eukaryotic cells. There are three described human vault RNAs from a cluster on chromosome 5 [15]. Stadler et al. [16] reported differential vault RNA expression in five human cancer cell lines and consensus patterns of small RNAs from vault RNAs across species. Vault particles are associated with multidrug resistance and intracellular transport. Persson et al. [17] discovered that human vault RNAs produce several small RNAs via mechanisms different from the canonical miRNA pathway, but at least one such small RNA associates with Argonaute proteins and guides sequence-specific cleavage of mRNAs to regulate gene expression. In particular Persson et al. discovered regulation of CYP3A4 (one of 57 human cytochrome P450 proteins) in MCF7 cells by a small byproduct of vault RNA transcription. The CYP3A4 enzyme is important in the initial metabolism of many marketed drugs [18]. Importantly, the experiments of Persson et al. might explain the association of abundance of vault particles with drug resistance.
It seems quite likely that nature must put such abundant, selected subsequences of the above types to some purpose, implying unrevealed pathways that are presently without definitive annotations or even realisation [3]. For example, nuclear-localized small RNAs might be epigenetic regulators of gene expression [9]. Thus block patterns of small RNA transcription sources might greatly improve and simplify ncRNA annotation [13].
Regarding neurological bioactivity, Smalheiser et al. [19] discovered in adult mouse hippocampus that certain species of 25-to 30-nt small RNAs derived from specific sites within well known noncoding RNAs were dramatically increased as a consequence of odorant discrimination training. This work reveals the potential importance of byproducts of ncRNA synthesis in neuroscience, possibly a universe of gene regulation parallel to that of the miRNAs. Consistent with the above prior work, we found that reads representing the 16-nt sequence S appear in every sample more than ten times as frequently as reads from the other three RNU1 hairpins and at frequencies comparable to those of abundant brain miRNAs. Further compounding interest in the 16-nt sequence S from hairpin H are, in the manner of miRNA target predictions, two putative target regions (lengths 9 and 11 nts) in the 3'UTR of splicing regulator gene SFRS1. Thus the 16-nt byproduct of RNU1 synthesis (from promotion of splicing) might also inhibit expression of SFRS1 (inhibition of splicing or at least inhibition of formation of spliceosome components). This might be a form of auto-regulation essential to homeostasis of splicing. Our neuroscience interests provide focus on SFRS1 protein product because it modulates several forms of synaptic plasticity considered to be involved in the very essence of memory [20].
Thus there is a triple intersection of bioinformatics: annotated function of an ncRNA, abundance in brain of a small RNA evidently processed from the same ncRNA source, and sequence alignment of the complement of the same small RNA with the 3'UTR of a major gene having the same function. These in silico coincidences demand investigation of potential miRNA-like mechanisms involving the RNU1 hairpin H, especially with regard to SFRS1. Needed are functional validations of nuclear RNU1 targets. Considering the huge impact of splicing function in nature and dysfunction in disease, elucidation of splicing homeostasis would carry a significant potential for progress toward novel diagnostic tools and drug platforms.
Regarding RNU1 context, hairpins studied by O'Gormann et al. [2] (which do not include S) were found to be bioactive, as mentioned above. Additionally, it has long been known that pre-mRNA splicing can be regulated both positively and negatively by reversible phosphorylation of spliceosomal SR proteins [21,22]. Thus it would be no surprise that additional layers of complexity might exist to regulate bioactivity of SFRS1 protein.
Moreover, Kohtz et al. [23] showed at an early date that SFRS1 protein cooperates with U1 small nuclear ribonucleoprotein particle (snRNP) in binding pre-mRNA, so there is already a direct, mechanistic link of RNU1 in U1 with SFRS1 protein. However, demonstrating that a small RNA byproduct of RNU1 transcription goes on to bind to SFRS1 mRNA and inhibit expression of that gene would be, to our knowledge, a novel splicing feedback loop discovered by virtue of modern, unbiased sequencing.
In summary, alignments of abundant reads, hairpin structures, and logical targets are known to be important in some cases and as yet unrecognised alignments are likely to be important in others-provided such colligations can be efficiently discovered.

Results and discussion
We note that the topic of processing small RNA sequencing results is a very active area of research with several important, powerful search and alignment engines developed along lines of analysis somewhat different from ours. These are represented by development of software for efficient and fast selection of abundant core sequences within numerous short reads by Hoffmann et al. [24], the description of novel miRNA discovery methods with mirTools from Friedländer et al. [25], and comprehensive statistical and annotative methods in mir-Tools by Zhu et al. [26] and Gunaratne et al. [27].
The conformations of the full RNU1 molecule are predicted by mfold [28] to include four hairpins, of which the 3' end hairpin H is of interest due to sequencing abundance and high predicted hairpin stability. The nascent RNU1 transcript is presumably chaperoned by proteins, but this hairpin might be so stable that it immediately forms and remains folded in most or all RNU1 conformations. Regardless, the strong signal for S in RNA from human brain and from our custom assays of human neural stem cells suggests that some mechanism isolates S from the hairpin H and protects it as a 16-nt single-stranded RNA from digestion.
In summary, we advocate development a rigorous methodology leading to the general discovery of multiple alignments of the S type as depicted in Figure 1.

Conclusions
As expressed by Kawaji et al. [3], nature seems to shun dogmatic classification of small biological RNA molecules. They point out that the likelihood of unrevealed pathways, implied by discovery of abundant small "core" small RNA sequences derived from deep sequencing 2D hairpin conformations of noncoding RNAs proteins complexed with ncRNAs in ribonucleoprotein particles and displaying targetable 3'UTRs in mRNAs RNA "cores" appearing in RNA hairpin stems small RNA sequences targeting 3'UTRs ncRNAs and proteins forming ribonucleoprotein particles Figure 1 Three types of bioinformatic data linked pairwise in a triangle. "Core" RNAs are short sequences appearing abundantly in intersections of sets of deep sequencing reads. Noncoding RNAs (ncRNAs) generally have predicted 2D structures and certain well-formed hairpins of some 20 to 100 nt, possibly with cores embedded in stems. The same ncRNAs can also be associated with proteins in ribonucleoprotein particles (RNPs). Lastly, mRNAs coding the same RNP proteins have 3'UTR regions that might be targeted in an miRNA-like manner by short RNAs, including some core RNAs.
subsequences from particular regions of ncRNAs, requires us to avoid premature definitions, specifications, and annotations of the transcriptome and its regulation. The present paper reports a remarkable colligation, a suggestive triple alignment of disparate bioinformatic observations, and it points to a line of investigation regarding the stable abundance of an organelle that is a hallmark of metazoan diversity, the spliceosome. What other colligations that include ncRNA fragments remain to be discovered? How can they be efficiently and systematically discovered in the contexts of other organelles and cell functions? These questions, partly answered by the research of others and the present work, demand attention and resources appropriate to elucidation of the foundations of cell biology.

Methods
We have sequenced Dicer-processed small RNAs from dorsolateral prefrontal cortex (Brodmann area 9) of three persons who had no mental illness at time of death. Samples were generously made available from the Stanley Medical Research Institute. The spliceosome is involved in synaptic plasticity, critical to normal brain function. It is of interest that altered expression of SFRS1 has been reported in post-mortem schizophrenia studies involving variable isoforms of DISC1 [29] or NCAM1 [30]. We discovered a 16-nt sequence-herein denoted as S = TTCGCGCTTTCCCCTG or UUCGCGCUUUCCC-CUG-in thousands of reads, and we confirmed by PCR the presence of S in commercially available neural stem cells derived from human embryonic stem cells. S includes 9 and 11 nucleotide subsequences that appear exactly as the reverse complements of two subsequences of the 3'UTR of the important splicing gene SFRS1. Furthermore, S is included as part of the loop and the 3' side of a well-formed 28-nt hairpin H from the snRNA RNU1. Parallel to other studies cited above, this level of complementary suggests that S could bind to the 3'UTR of SFRS1 and inhibit processing or accelerate sequestration or degradation of SFRS1 mRNA, possibly in an miRNA-like manner.

Methods of sequence analysis
As yet undiscovered regulatory small RNAs might function in normal human brain or might provide signatures of human brain diseases. Total RNA from dorsolateral prefrontal cortex was isolated using Trizol under the auspices of Stanley Medical Research Institute. We derived cDNA libraries using the Illumina Small RNA Sample Preparation Kit (San Diego CA), following the manufacturer's version 1.5 protocol. The cDNA libraries were sequenced with an Illumina Genome Analyzer at our university core facility. Resulting 35-base reads were first filtered to remove concatenated adaptors at the 3' ends. Specifically, we sought the first six or more bases of the Illumina 3' adapter TCTCGTATGCCGTC TTCTGCTTGAAA... in valid reads. Next a variable number of consecutive As were trimmed from 3' ends. Those remaining sequences with 16-29 bases were aligned with ClustalW2 [31] and then further filtered by requiring at least three exact copies. Frequently appearing subsequences of length ≥16 bases were designated "core sequences," including S itself. Our algorithm was a modification of published methods [12,32]. Sources of cores from all three subjects were about 64% from mature miRNAs, 27% from unknown transcripts, 4% from ends of tRNAs, 3% from snoRNAs, 1% from snRNAs (including S), and 1% from mitochondrial tRNAs.
We generated 1000 distinct, random permutations of the 16 bases corresponding to S and sought them among all raw reads of all three samples. Not one was found. The probability of finding not even one of 1000 16-base random sequences within a million, distinct, random, 30-base sequences is 0.0304 (but of course read sequences are not random). Thus while S itself is plentiful among our reads, distinct permutations of the bases corresponding to S yield random sequences that are consistently absent.
We used BLAT in UCSC Human Genome Browser and found full 16-base sources of S to occur only when preceded by TGCA or TGCG and only in 14 genomic locations, all 14 listed within the RNU1 genes in wgEncodeGencodeAutoV3.
Given our interest in neuroscience, we sought to determine if human embryo-derived neural stem cells also generated copies of S. Cells used in these experiments were commercially available neural stem cells called hNP1 cells from ArunA Biomedical (Athens GA) derived from human embryonic WA09 cell lines. hNP1 cells grow as an adherent monolayer and are homogeneous, uniformly expressing various neural stem cell proteins (e.g. NES, SOX2) and low levels of embryonic stem cell proteins (e.g. POU5F1 (alias OCT4)) [33,34]. hNP1 neural stem cells are subject to rigorous quality control including DNA fingerprinting, viral testing, and maintenance of a stable karyotype.
We cultured hNP1 cells according to the manufacturer's protocol. We harvested about 1 million passage 4 cells. Cytoplasm and nuclear lysate were separated using the Norgen Cytoplasmic and Nuclear RNA purification kit (Thorold ON). We verified partitioning with an Agilent Bioanalyzer (Santa Clara CA) by observing with a DNA chip a DNA:DNA concentration ratio of >100X of nucleus over cytoplasm and for the same samples processed with a Bioanalyzer RNA chip an RNA:RNA ratio of >10X for 18S peak in cytoplasm over nucleus. These two ratios for the same samples imply good nuclear and cytoplasmic partitioning.
We used a Qiagen (Alameda CA) miScript PCR system with custom primer to determine the presence of S in the samples. For negative control, we used a miScript primer for a 16-nt subsequence Arabidopsis thaliana miRNA-159a (mature previously shown by us to be absent in nuclear and cytoplasmic extracts from the same cell type, data not shown) and a 16-nt subsequence of human miRNA-128 (mature previously shown by us to be present in both nuclear and cytoplasmic fractions, data not shown). Input RNA templates were 0.76 pg/uL for nuclear extract and 1.01 pg/uL for cytoplasmic extract. All 12 samples were run simultaneously using an Applied Biosystems (Foster City CA) 7900HT Fast Real-Time PCR System. As shown in Figure 2, we readily and consistently detected S in both nuclear and cytoplasmic fractions from hNP1 cells. All negative controls were not detected or very weakly detected (cycle count >35).

Methods of RNA targeting analysis
Using TargetScan 5.1 [35] we found the in silico predicted targets for all ten of the seven-nt subsequences of S. Repeatedly appearing was SFRS1. This was because 11 consecutive nts from S were the exact reverse complement of a 3'UTR sequence of SFRS1; a second target was generated from 9 nts. How S might bind to the 3'UTR is shown in Figure 3; high C-G bond content suggests strong affinity.
We note that with one exception (miR-4315 with seed CGCUUUC) no seven-nt subsequence of S is also the seed of any human miRNA (nts 2-8 of all 1,100 miRBase v15 mature miRNAs [36]); thus targeting by S would likely not be redundant to conventional miRNA actions. Also to be noted is the fact that overall alignments are well within the parameters of potential alignments as reported by Thomas et al. [37]. SFRS1 protein is employed in protein-protein interactions and other processes and in particular recruits the U1 snRNP to the 5' splice site [38,39]. The upshot is the suggestion that a byproduct of RNU1 transcription is auto-regulation of spliceosome assembly and function.
Importantly, Ohrt et al. [40] have reported evidence of a nuclear RISC imported from the cytoplasm and consisting of Ago2 and a mature miRNA, thus some 20X smaller than the conventional cytoplasmic RISC. It therefore is possible that the small sequence S is also mounted in a nuclear RISC that includes Ago2.

Methods of RNA folding analysis
RNA hairpin structures are key features of RNA functions and processing generally, and miRNA processing in particular. The mfold engine uses dynamic programming and large tables of empirically derived binding affinities for small dsRNAs [28]. mfold predicts four hairpins formed from the RNU1 consensus (shown in Figure 4). mfold also predicts five 2D conformations for the full RNU1 sequence, four of which include the hairpin H with S as loop and 3' side; one full conformation is shown in Figure 5.
Regarding stability of the various hairpins, the stability ratio defined by -Gibbs free energy divided by number of nts (length) of a hairpin is -dG/L. For the fourth hairpin H in isolation, the number is much stronger (.59) than that of any other RNU1 hairpin (average~.35).

Methods of bioinformatic control mechanism analysis
Again, integration of biochemical and bioinformatic data to yield information about mechanisms of control of gene expression is the goal of this paper. To that end, the above example provides a paradigm for seeking causal relationships among biochemical concepts and bioinformatic concepts. Figure 6 will be described below to substantiate the paradigm.
Deep sequencing applied to cDNA libraries constructed from small RNA molecules provided us with millions of reads of observed sequences. From tables of Illumina 35-base reads, we trimmed from 3' regions the 5' ends of adaptors (at least six bases) and zero to four As, and then retained only the results that appeared at least three times. In our procedure, this yielded thousands of distinct sequences of 16 to 29 bases, each with an instance count ranging from three to more than 100,000. The sequences were aligned with ClustalW2, and then commonly occurring subsequences of at least 16 bases were deduced, each with a total instance count; we called each such subsequence a "core." This procedure includes multiple tuned parameters and an inevitable degree of arbitrariness. Conclusions reached with an algorithm containing tunable parameters must always be tested for invariance with respect to reasonable retuning.
In a separate line of investigation, tables of ncRNAs can be processed by submitting sliding windows of 50-nt subsequences to RNA folding engines such as mfold. Sought are simple RNA hairpins incorporating at least 25 nts and meeting certain stability criteria. For example, the above ratio of -Gibbs free energy (kcal/mol) to length (number of nts) in a hairpin should be at least~0. 30   >RNU1 consensus (from 10 very similar RNU1 sequences from 14 that include S): AUACUUACCUGGCAGGGGAGAUACCAUGAUCACGAAGGUGGUUUUCCCAGGGCGAGGCUUAUCCAUUGCACUCCGGAUG UGCUGACCCCUGCGAUUUCCCCAAAUGUGGGAAACUCGACUGCAUAAUUUGUGGUAGUGGGGGACUGCGUUCGCGCUUU CCCCUG Figure 4 Layout of predicted hairpins and observed cores in an RNU1 consensus sequence. Underlined sequences correspond to predicted 2D hairpins, as declared by mfold and approximately as predicted and described by O'Gormann et al. [2]. Nucleotides highlighted in gray correspond to our declared "cores" of Illumina sequences (abundant shared subsequences in Illumina reads). In Sample 1, all cores but S had 45 to 154 instances, but S itself had 4736 instances. In the other two samples the counts were: 10-59 vs. 1317; 60-247 vs. 2454. In short, S was much more popular than any other sequenced cores apparently derived from RNU1.

Figure 5
One of the 2D structures of RNU1 predicted by mfold [28].  Figure 6 Network concepts embedded in a flowchart. The goal is efficient integration of selected types of information from sequencing, RNA folding, and gene function data. Sequencing reveals abundant small RNA sequences (blue box). Inputs for bioinformatic assays are tables of ncRNAs, proteins with which ncRNAs are complexed to yield gene processing moieties, and 3'UTR sequences of the mRNAs for the same proteins (gray boxes). Then biochemical assays and bioinformatic assays are compared to suggest agents of regulation (yellow box) and control mechanisms (green box). similar to mature miRNAs. A survey by Zhang et al. of 513 pre-miRNAs found average MFEI was 0.97, significantly higher than sampled mRNAs (0.65), tRNAs (0.64), or rRNAs (0.59). The distinguished hairpin in Figure 5 has MFEI = 0.85; all other predicted hairpins from RNU1 have MFEI ≤ 0.71. Thus sequencing and hairpin searches can lead to comparisons of cores with sides of predicted hairpins to deduce potential agents of gene regulation, possibly processed in an miRNA-like manner.
In a third line of investigation, ncRNAs can be considered for membership in well annotated ribonucleoprotein particles of many types with wide-ranging functions [42]. The proteins in the same particles can be examined for subsequences in the 3'UTRs of their mRNAs that align with the putative miRNA-like agents formed from processed cores.
The complete flowchart of the proposed triple colligation is shown in Figure 6.