Structural analysis of the lncRNA SChLAP1 reveals protein binding interfaces and a conformationally heterogenous retroviral insertion

Evolutionary analysis reveals conserved and potentially functional exons of SChLAP1

Conserved RNA motifs within prokaryotic and eukaryotic transcripts are often indicative of critical function. Thus, we performed a BLAST search of the NCBI RefSeq database to identify potential homologous SChLAP1 sequences in nonhuman genomes (Zhang et al. 2000; O'Leary et al. 2016). This search revealed several nonhuman primate (NHP) sequences of SChLAP1, all of which were predicted to be noncoding but have no functional annotation to date (Fig. 1A; Supplemental Table S1). SChLAP1 sequences were not called in any nonprimate species, suggesting that SChLAP1 emerged during primate evolution, consistent with general evidence of greater lncRNA content with increasing organismal complexity and evidence that many lncRNAs are primate specific (Derrien et al. 2012; Liu et al. 2013; Mattick et al. 2023). To further evaluate SChLAP1 conservation, we examined phyloP scores for SChLAP1 in the UCSC Genome Browser (Kent et al. 2002; Pollard et al. 2010). Beginning with a 447-mammal set, we observed generally near-zero phyloP scores throughout the transcript, indicative of poor conservation (Fig. 1B; Supplemental Appendix). In contrast, phyloP scores were generally more positive in a more restricted 30 mammal (27-primate) set, indicative of conservation (Fig. 1B; Supplemental Appendix). These results suggested that SChLAP1 shows stronger conservation when comparing within the primate lineage and further supports SChLAP1 as a primate-conserved lncRNA. However, in both phyloP analyses, we observed generally negative phyloP scores for exon 1 of SChLAP1, indicative of accelerating evolution in this region of the transcript.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 1.

Phylogenetic analysis of SChLAP1. (A) BLAST search of isoform 4 of SChLAP1 (containing all possible human exons) against the RefSeq database. (Top) Isoform 4 of SChLAP1 with each exon shown for reference. The 3′-positions of each exon in the human SChLAP1 are labeled. (Bottom) Sequence alignments of all identified SChLAP1 homologs against human SChLAP1 isoform 4. The 5′- and 3′-positions of alignment against human SChLAP1 are noted. All BLAST alignment scores were >200. Asterisk denotes Pan troglodytes aligning at nucleotides 1006 rather than 1007. (B) PhyloP analysis of SChLAP1 isoform 4 using phyloP 30-way (top) or 447-way (bottom). Figure generated in Integrative Genomics Viewer (Robinson et al. 2011). (C) RepeatMasker analysis of SChLAP1 Iso. 1 (see also Supplemental Fig. S1).

While our BLAST search uncovered several complete putative SChLAP1 homologs in NHPs, we noted that multiple other genes (both human and NHP, including protein-coding and noncoding genes) aligned with exon 1 of SChLAP1 (not shown). This result suggested that repeat insertions, occurring not only within SChLAP1 but also in other genomic regions, may have resulted in spurious alignment in our BLAST search. These insertions (i.e., transposable elements [TEs]) are highly promiscuous DNA elements that can replicate independently from the host cell. TEs can become incorporated into coding and noncoding genome sequences and result in the formation of novel functions, such as protein recognition and/or RNA structures in lncRNAs (Johnson and Guigó 2014). To search for potential TE insertions, we analyzed the SChLAP1 Iso. 1 sequence using RepeatMasker (http://www.repeatmasker.org), focusing on SChLAP1 Iso. 1 hereafter as it is the most abundant isoform of SChLAP1 (Prensner et al. 2013). We identified two long terminal repeats (LTRs) within SChLAP1: LTR12C (DFAM accession DF0000402) within exon 1 (nucleotides 1–322; 95.3% of exon 1) and THE1B (DFAM accession DF0000818) within exon 7 (nucleotides 1123–1436; 46.9% of exon 7; Fig. 1C; Supplemental Fig. S1; Storer et al. 2021). LTRs occur at the 5′- and 3′-ends of retroviral genomes and are involved in viral processes such as replication. After integration into the host genome, the central components of the retroviral genome may be removed across evolution by recombination, such that only the LTR sequences remain (so-called solo-LTRs) (Babaian and Mager 2016; Johnson 2019). LTRs constitute over 8% of the human genome, and these sequences may be epigenetically silenced or, in contrast, “exapted” (i.e., co-opted) by the host genome for various functionalities including as promoters/enhancers (Babaian and Mager 2016; Johnson 2019). The occurrence of these sequences in SChLAP1 suggested that some of SChLAP1 function may be the result of such exaptation from these insertions and resultant RNA sequence/structure.

Our BLAST search results are consistent with evidence of retroviral insertion in exons 1 and 7. LTR12C, derived from human endogenous retrovirus 9 (HERV-9), is a Hominoidea (i.e., ape)-specific insertion. Our BLAST search only observed this sequence in Homo sapiens, Gorilla gorilla gorilla, Pan troglodytes, P. paniscus, Nomascus leucogenys, and Hylobates moloch genomes, all of which are of Hominoidea phylogeny (Schoch et al. 2020; Storer et al. 2021). In comparison, the BLAST results did not identify LTR12C insertions nor the SChLAP1 exon 1 sequence within the Colobus angolensis palliates, Cercocebus atys, or Mandrillus leucophaeus genomes (Fig. 1A); these results are consistent with known phylogeny as these species are Old World monkeys and thus of non-Hominoidea phylogeny (Schoch et al. 2020). The other retroviral insertion in SChLAP1 (THE1B) belongs to the MaLR (mammalian apparent LTR-retrotransposons) family of repeat elements. In contrast to LTR12C, THE1B is specific to Simiiformes (i.e., simians), which contains both Hominoidea and Old-World monkeys (Schoch et al. 2020; Storer et al. 2021). In turn, a complete exon 7 sequence was observed in both lineages in our BLAST search. Our results indicated that THE1B resulted in the formation of approximately half of exon 7. Given that both LTR12C and THE1B insertions are selectively found in Hominoidea and Simiiformes, respectively, we believe these data may explain why SChLAP1 did not show evidence of conservation outside of the primate lineage. TE insertions also explain the low/negative phyloP scores (i.e., accelerating evolution) in these regions of the transcript.

Using the USCS Genome Browser (Kent et al. 2002), we found that the LTR12C sequence alignment to the 5′-end of human SChLAP1 not only overlaps with the majority of SChLAP1 exon 1 but also includes the promoter region of the SChLAP1 gene (Supplemental Fig. S2). Potential retroviral incorporation within the SChLAP1 promoter has been mentioned in previous work (Prensner et al. 2013; Babaian and Mager 2016), and this analysis suggested the insertion may have also introduced novel sequence/structures into the 5′-end of the SChLAP1 transcript.

Sequence and bioinformatic analyses identify potentially functional SChLAP1 substructures

To complement our phylogenetic analysis, we used additional sequence and bioinformatic analyses to identify regions of SChLAP1 that may be crucial to its function. As above, we focused our structural analysis on SChLAP1 Iso. 1. To begin, we used ScanFold 2 (Andrews et al. 2022) to predict thermodynamically favorable and likely functional RNA secondary structures within SChLAP1. Specifically, ScanFold 2 folds a given transcript in sliding folding windows, shuffles the sequence of each window (100 times each in our case), and thereafter calculates the folding energy of each shuffled sequence. The folding energies for the native versus shuffled sequences are compared via Z-score (defined as the average minimum free energy [MFE] of the shuffled sequences subtracted from the MFE and divided by the standard deviation), where Z-scores of −1 and −2 indicated a folding window 1 and 2 standard deviations more stable than the shuffled folds, respectively. Previous work has found that regions with significantly lower energies of folding in the native versus shuffled sequence are often enriched in functional roles (Andrews et al. 2018, 2022). From this work, we identified multiple SChLAP1 structures with significantly low Z-scores (Fig. 2, top), including (1) predicted stem–loops in 5′-end of SChLAP1 (i.e., containing the LTR12C insertion); (2) the junction between exon 2 and exon 5 of SChLAP1 (hereafter referred to as the E2–E5 junction); (3) a predicted helix from nucleotides 989–1056 (within exon 7); and (4) an adjacent predicted three-way junction (3WJ, nucleotides 1141–1213) and stem–loop (nucleotides 1224–1242), both located within exon 7 and the THE1B insertion. We obtained concordant structure predictions with three independent runs of the ScanFold 1 algorithm on the SChLAP1 Iso. 1 sequence (Supplemental Fig. S3; Andrews et al. 2018). These analyses support the presence of potentially functional RNA structures within SChLAP1 on the basis of their favorable sequence-based folding properties. The above RepeatMasker analysis also suggested that some of these structures may have been derived through retroviral insertion.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 2.

ScanFold 2 (Andrews et al. 2022) prediction of SChLAP1 Isoform 1. (Top) Significant base pairs identified using ScanFold 2. Second and third row depict average base-pair Z-score and minimum free energy of folded windows calculated by ScanFold. Bottom row shows GC content calculated over 51 nt windows. Exon coordinates of SChLAP1 Iso. 1 are depicted along with coordinates of LTR12C and THE1B insertions. Figure generated in Integrative Genomics Viewer (Robinson et al. 2011).

Additionally, we calculated the GC content of SChLAP1 in 51 nt sliding windows to identify regions of structure correlated with the more stable hydrogen bonding provided by canonical GC base pairs versus AU base pairs (Fig. 2, bottom). Overall, SChLAP1 comprises 44.5% GC but contains windows as high as ∼70% GC. Several ScanFold predicted structures occur in such high-GC regions; however, we observed that one highly predicted structure (nucleotides 989–1056) occurs within a remarkably GC-depleted region of the transcript. We refer to this predicted structure as the AU-rich helix where it appears in our analysis below.

Development of a secondary structure model for SChLAP1 in vitro using SHAPE-MaP

Our work so far has indicated that there are thermodynamically predicted functional structures within the SChLAP1 RNA transcript, some of which were derived from retroelements inserted into the primate lineage. We next set out to generate an experimentally informed secondary structure model of SChLAP1 using SHAPE-MaP. We first began by generating an in vitro structure SChLAP1 model as a preliminary reference to evaluate the propensity of these structures to form both in vitro and in cells (analyzed below). To this end, we performed SHAPE-MaP and DMS-MaP on in vitro transcribed SChLAP1 Iso. 1 isolated using a semi-native purification protocol, where both heat denaturing and harsh buffer exchanges are avoided, as previous work found this protocol superior in maintaining a homogenously folded RNA as compared to denaturing protocols (Somarowthu et al. 2015; Adams et al. 2019). We chose the SHAPE reagent 5-nitroisatoic anhydride (5NIA) for its enhanced cell permeability compared to other SHAPE reagents and precedence for use of 5NIA in prostate cancer cells (Busan et al. 2019), which enabled comparable use of 5NIA in both in vitro and in cellulo experiments.

The in vitro secondary structure model of SChLAP1 was generated with SHAPE and DMS probing, which were performed with independent batches (i.e., transcribed, purified, and probed on different days) of SChLAP1 RNA. As SChLAP1 is >1400 nt, which hinders sequencing analysis on Illumina-based platforms, we generated four overlapping amplicons (positions 1–500, 400–903, 800–1300, and 1200–1436) from the probed RNA to facilitate sequencing. These four amplicons were independently reverse transcribed (SuperScript II and Mn²⁺ for SHAPE probing, TGIRT for DMS) and amplified for sequencing. The SHAPE sequencing data were processed using the ShapeMapper (Busan and Weeks 2018) and SuperFold (Smola et al. 2015b) pipelines, whereas the DMS-MaPseq program was used to process the DMS data (see Materials and Methods; Zubradt et al. 2017). In both SHAPE and DMS experiments, we observed a lack of chemical reactivity in a poly(A) stretch of SChLAP1 (position 1088–1104), which is consistent with work from Kladwang et al. showing that chemical modifications in poly(A) regions are bypassed by reverse transcriptases and result in incongruous chemical modification frequency and mutational profiling results (Kladwang et al. 2020). To avoid biasing our resultant model as predicted by SuperFold/RNAstructure, we manually set this poly(A) stretch to “undefined” in the input .shape and .map files for all SHAPE-informed structure predictions of SChLAP1; setting these nucleotides to undefined removes all chemical reactivity data from these nucleotides, preventing these data from biasing the predicted secondary structure.

The MFE model (Supplemental Fig. S4A) and predicted base pairs (Supplemental Fig. S5A) of in vitro SChLAP1 reveal a complex structural landscape, including significantly base-paired regions of varying secondary structure types as well as large single-stranded regions. Pseudoknot prediction was performed as published (Smola et al. 2016; Wan et al. 2022), but predictions were not incorporated into the structure models reported herein due to lack of evidence for pseudoknot formation in SChLAP1 in our downstream in cellulo analysis (see Supplemental Text 1; Supplemental Fig. S6). The DMS data showed good agreement with the in vitro SHAPE data (Supplemental Fig. S4A) and resulted in greater modification of adenosines and cytidines over guanosines and uridines as expected (Supplemental Fig. S4B). We also observed agreement for several local structures between our SHAPE-informed partition functions and ScanFold models (Supplemental Fig. S5A). These include the aforementioned AU-rich helix, the 3WJ and adjacent stem–loop, and the E2–E5 junction. The agreement between these models gives independent support for the formation of these various substructures in in vitro SChLAP1 as predicted by ScanFold 2.

Ex cellulo SHAPE probing of SChLAP1 reveals cell-derived RNA structures

As our in vitro chemical probing data supported the presence of functionally important structures within the SChLAP1 RNA, we next set out to generate a secondary structure model for cell-derived SChLAP1 RNA to identify whether these structures or others fold in a physiologically relevant context. We performed SHAPE probing of endogenous SChLAP1 isolated from LNCaP prostate cancer cells, a common model cell line with high SChLAP1 expression as compared to normal prostate cells (Prensner et al. 2013) and with a known metastatic phenotype. Due to the retroviral insertion in exon 7, which overlapped completely with amplicon 4 used in our in vitro probing, we extended the coordinates of amplicon 4 from what was used in our in vitro experiments to avoid amplification of other THE1B-containing transcripts (see Materials and Methods).

We began our in-cell structure modeling with RNA that was gently isolated from LNCaP cells and Proteinase K treated (i.e., ex cellulo RNA) using previously established protocols to extract cellular RNA without heat-denaturing (Smola et al. 2015a, 2016). Two independent replicates of ex cellulo probing of SChLAP1 yielded reactivity profiles with good agreement (Pearson's R = 0.72). The resultant secondary structure arc diagrams (Fig. 3A; Supplemental Fig. S5B,C) and MFE model (Fig. 4, see below) revealed a complex structure for cell-derived SChLAP1, with an array of predicted structures across the length of the transcript. We also observed agreement among several structures identified in our ex cellulo SHAPE data, in vitro SHAPE data, and ScanFold data (Supplemental Fig. S5), including the E2–E5 junction, AU-rich helix, and 3WJ, giving independent support for the formation of these structures and their ability to form in cells.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 3.

SChLAP1 secondary structure generated through ex cellulo SHAPE-MaP from LNCaP cells (representative of two independent replicates). (A) Arc diagram showing the SHAPE-informed base-pair predictions across ex cellulo SChLAP1. Arc diagram of second replicate shown in Supplemental Figure S5C. (B) Median SHAPE reactivity of 51 nt windows across the length of SChLAP1 subtracted from the global median (SHAPE = 0.31). (C) Shannon entropy plot across SChLAP1. Blue boxes denote low SHAPE/Shannon regions that were reproducible across both biological replicates. Figure generated in Integrative Genomics Viewer (Robinson et al. 2011).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 4.

ΔSHAPE analysis of SChLAP1. Significant in-cell protected or enhanced regions identified in both replicates of ΔSHAPE (see Materials and Methods) are mapped over a representative MFE structure of ex cellulo SChLAP1. SChLAP1 structure visualized in VARNA (Darty et al. 2009).

We further analyzed our ex cellulo SHAPE data to identify putative functional structures within cell-derived SChLAP1. Work by Siegfried et al. characterizing HIV-1 genomic RNA found that RNA structural elements classified as (1) highly structured (i.e., low SHAPE reactivity) and (2) well-determined (i.e., low Shannon entropy), which are referred to as low SHAPE/Shannon (lowSS) regions, were significantly enriched in previously unknown functional roles (Siegfried et al. 2014). This metric has since been used to identify novel functional regions, including in the Dengue virus RNA genome (Dethoff et al. 2018), XIST lncRNA (Smola et al. 2016), and HCV genome (Wan et al. 2022). To identify if SChLAP1 contained lowSS regions, we calculated the SHAPE reactivity (51 nt local median compared to global median, Fig. 3B) and Shannon entropy (51 nt windows, Fig. 3C) to identify lowSS sites. In our SChLAP1 model, we identified three lowSS regions: nucleotides 1–165 (within LTR12C insertion), nucleotides 398–478 (E2–E5 junction), and nucleotides 989–1056 (AU-rich helix) (Fig. 3, blue bars). These lowSS regions agree well with ScanFold-predicted structures (Fig. 2).

Additionally, we observed that the nucleotides between the above-mentioned lowSS regions (i.e., nucleotides 166–397, which contains portions of exon 1 and 2, and nucleotides 479–988, which contains portions of exons 5 and 7 and the entirety of exon 6) display both elevated SHAPE reactivities and Shannon entropies relative to the median (Fig. 3B,C). This result indicated that these regions are generally single-stranded and/or conformationally heterogenous. Consistent with these results, the arc diagrams of the two ex cellulo probing replicates vary considerably in these regions (Supplemental Fig. S5C). In contrast, downstream from the lowSS site located at nucleotides 989–1056 (i.e., nucleotides 1057–1436), we observe low SHAPE reactivity alongside high Shannon entropy (Fig. 3B,C), indicating a highly structured but poorly predicted and/or conformationally heterogenous structure. This highly structured region correlates well to the THE1B insertion within SChLAP1 (nucleotides 1123–1436).

We next compared the consistency of our in vitro and ex cellulo SHAPE profiles using Pearson correlation coefficients. We calculated the correlations using 51 nt sliding windows for both replicates of our ex cellulo SChLAP1 against our in vitro SHAPE data. We observed regions of high correlation between our ex cellulo and in vitro data alongside regions of low/modest correlation (Supplemental Fig. S7). In particular, SChLAP1 nucleotides 1–204, 318–477, and 960–1105 were the largest regions with high correlation between our in vitro and ex cellulo data (defined as having Pearson R ≥ 0.5 in both replicates, a threshold chosen heuristically based on previous comparisons of in vitro and in/ex vivo chemical probing studies [Sherpa et al. 2018; Frank et al. 2020; Manfredonia et al. 2020]). Thus, these regions are more likely to recapitulate a similar structure in vitro and ex cellulo. These regions also align with the aforementioned lowSS regions in our ex cellulo data (Fig. 3). We also identified smaller regions that show correlation between the in vitro and ex cellulo models, particularly within nucleotides 803–857, which corresponds to a predicted stem and 5′-end portion of a large, single-stranded loop in our in vitro (Supplemental Fig. S4A) and ex cellulo (Fig. 4) MFE models.

Identification of protein binding regions via in cellulo SHAPE

SChLAP1 has been documented to interact with several proteins (see Introduction). As our in vitro and ex cellulo work identified potentially functional secondary structures within SChLAP1 (in particular, ScanFold and lowSS regions) as well as single-stranded regions that could additionally function in protein binding, we next performed in cellulo SHAPE probing to identify whether these structures persist in cells and/or function as putative protein binding regions that may be involved in any of the aforementioned interactions. For consistency with the ex cellulo experiments, in cellulo SHAPE probing was performed in the LNCaP cell line. The resultant SHAPE reactivity profiles between two independent replicates showed good agreement (Pearson's R = 0.59). For reference, replicate in cellulo SHAPE probing experiments of endogenously expressed Xist lncRNA yielded a Spearman correlation coefficient of 0.50 (Smola et al. 2016). We then used the ΔSHAPE algorithm to identify protein binding regions/structures of SChLAP1 by comparing our in cellulo reactivity data and ex cellulo data, as the ΔSHAPE algorithm has previously identified protein binding sites in other lncRNAs and viral RNAs (Smola et al. 2016; Frank et al. 2020; Jones et al. 2020; Schmidt et al. 2020; Przanowska et al. 2022; Wan et al. 2022).

The results of two replicates of ΔSHAPE analysis are shown on a representative MFE structure from our ex cellulo probing (Fig. 4); each replicate (see Materials and Methods) is individually shown in Supplemental Figure S8. Across both replicates of ΔSHAPE analysis, we observed indications for protein binding (in-cell protections) in both structured regions and relatively unstructured regions as previously identified in our ex cellulo analysis. For example, the three lowSS regions in our ex cellulo model, that is, the beginning of exon 1/LTR12C insertion, the E2–E5 junction, and the AU-rich helix, each contain in cellulo protected nucleotides, indicating that these structures likely participate in intermolecular interactions and suggesting a role for these highly predicted RNA structures in protein binding. Additionally, we observed significant in-cell protections within the remainder of exon 5, which was found to be relatively high in SHAPE reactivity and Shannon entropy and suggestive of an unstructured, single-stranded region (Fig. 3). These data suggested that this region functions as a single-stranded “landing pad” for protein binding, with similar regions observed in other RNAs such as Xist (Smola et al. 2016; Weeks 2021). Lastly, we observed significant in-cell protections in the low SHAPE, high Shannon 3′-end of SChLAP1 (within THE1B insertion), indicating that this region plays important roles for protein binding. In particular, these protections were observed in the aforementioned 3WJ (supported by ScanFold [Fig. 2] as well as our in vitro/ex cellulo predicted base pairs [Fig. 3; Supplemental Fig. S5]) and a CU-rich single-stranded region (nucleotides 1313–1335, observed in our in vitro/ex cellulo SHAPE-informed MFE models [Fig. 4; Supplemental Fig. S4A]). As this region corresponds to the THE1B LTR insertion (Fig. 1C; Supplemental Fig. S1), these results suggested that the formation of this SChLAP1:protein recognition interface occurred de novo through a retroviral insertion within primate evolution.

Several putative protein binding regions identified herein are consistent with previous work from other groups, supporting the validity of our ΔSHAPE data. For example, in-cell protections observed in exon 2 are consistent with recent work from Ji et al., where binding between SChLAP1 exon 2 and HNRNPL facilitated activation of the NF-κB pathway in glioblastoma (Ji et al. 2019). In separate work, HNRNPL was found to promote prostate cancer cell growth, and knockdown of HNRNPL significantly diminished LNCaP cell growth but not RWPE-1 cell growth, indicating cancer-specific functions for HNRNPL (Fei et al. 2017). While exon 2 lacks the full poly(CA) consensus sequence found for HNRNPL via RIP-seq (Fei et al. 2017), a CA-rich region (nucleotides 341–349) overlaps with an in-cell protected region in our ΔSHAPE data (nucleotides 347–349) and may support weak or indirect interaction with HNRNPL. In addition, protein binding with exon 7 is in line with previous work from the Chinnaiyan group. Preliminary work from Sahu revealed that a deletion of 250 nt (from position 1001–1250, termed the Deletion 5 region) of SChLAP1 Iso. 1 inhibited SChLAP1-driven invasion and binding to the SWI/SNF complex (Sahu 2015). These coordinates overlap with two regions of in-cell protection, the former within the AU-rich helix and the latter within the predicted 3WJ in the THE1B insertion. Both sites are also contained within favorable ScanFold-predicted structures (Fig. 2). While interaction with the SWI/SNF complex is contested (Raab et al. 2019), the role of this particular region in protein recognition is supported by our work.

Analysis of unfolded regions in cellular SChLAP1 compared to ex cellulo

Alongside the above in-cell protections, we also observed regions of in-cell enhancement of reactivity (i.e., nucleotides more reactive in cellulo as compared to ex cellullo data) by ΔSHAPE, indicating the unfolding of some structures in cells as compared to the ex cellulo model. In particular, in-cell enhancements appeared to be enriched at the 5′-end of the transcript and more limited in the remaining sequence (Fig. 4; Supplemental Fig. S8B). Within the first 500 nt of SChLAP1 (i.e., the first 34.8% of sequence), 76 of 130 (58.5%) and 82 of 120 (68.3%) instances of in-cell enhancement are observed in replicates 1 and 2 of ΔSHAPE analysis, respectively (Supplemental Fig. S8B). While the ΔSHAPE algorithm emphasizes robust local changes in SHAPE reactivity (Smola et al. 2015a, 2016), we hypothesized that the abundance of in cellulo enhancements called by ΔSHAPE at the 5′-end of SChLAP1 may be reflective of larger scale structural changes. To analyze broader conformational changes between ex cellulo and in cellulo SChLAP1, we calculated the overall SHAPE reactivity change in 51 nt windows between our ex cellulo and in cellulo samples as was used with the Xist lncRNA (Smola et al. 2016). Several regions of significant SHAPE reactivity change were identified (Supplemental Fig. S8C). Of particular note, we observed significant enhancement of SHAPE reactivity of exon 1 in cells, which was consistent with the multiple ΔSHAPE-identified in cellulo reactivity enhancements we observed (Fig. 4; Supplemental Fig. S8B). These results suggested a broad structural unfolding of this region in cells, in contrast to the significant indications for structuredness of this region in vitro and in silico, including the ScanFold predictions (Fig. 2, top; Supplemental Fig. S3), high GC content (Fig. 2, bottom), in vitro probing data (Supplemental Figs. S4A, S5A), and ex cellulo probing data (Figs. 3, 4), where significant structure was indicated. Additionally, we identified several in-cell protections by ΔSHAPE in this region (Fig. 4; Supplemental Fig. S8B), indicating that local sites of in cellulo protection occur and suggesting participation in intermolecular interactions.

Similar to exon 1, we observed that the E2–E5 junction showed enhancement of SHAPE reactivity in cellulo by ΔSHAPE (Fig. 4; Supplemental Fig. S8B). However, our ΔSHAPE data indicated that the majority of in-cell enhancements in this region are located on the 5′-end of the E2–E5 junction. Despite these in-cell enhancements, the structure of the E2–E5 junction is still strongly predicted in the arc diagram for our in-cell data (Supplemental Figs. S8A, S9A). We examined the SHAPE reactivities of the largest strands in this structure (nucleotides 398–412 and 464–478, the former containing multiple in-cell enhancements) and observed statistically significant enhancement of SHAPE reactivity in cellulo compared to ex cellulo for the 5′-strand, but not the 3′-strand in this region (Supplemental Fig. S9B). In addition, we observed ΔSHAPE in cellulo protections in a small stem–loop in the E2–E5 junction (nucleotides 449–461, Fig. 4), suggesting that protein binding occurs within this structure despite the presence of in cellulo enhancements elsewhere. From these data, we hypothesize that protein binding within the E2–E5 junction occurs primarily on the 3′-portion of the structure but not the 5′-portion. As a result, the 5′-stem of the structure, while putatively base-paired to the 3′-stem in vitro and ex cellulo, is more accessible to chemical modification in cellulo. We note that this proposed disruption would be distinct from the broader unfolding observed for exon 1/LTR12C insertion.

While the 5′-end of SChLAP1 appears to be generally unfolded in cells, we observed a few ΔSHAPE in cellulo enhancements at the 3′-end of SChLAP1 (Fig. 4; Supplemental Fig. S8B). This result is supported by the presence of highly predicted base pairs at the 3′-end rather than the 5′-end of SChLAP1 in cellulo (Supplemental Fig. S8A). Thus, we believe that most of the structures at the 3′-end of SChLAP1, such as the AU-rich helix and 3WJ, remain folded in cells, allowing these RNA structures to function within cells as well-folded protein binding substructures.

In vitro folding analysis of SChLAP1 3′-end reveals conformational heterogeneity in the THE1B LTR insertion

As our above work characterized SChLAP1 secondary structures and identified protein binding interfaces, we next set out to evaluate if any of these recognition elements participated in tertiary/higher order structural complexity. We specifically focused on the 3′-end of SChLAP1 for several reasons. First, this region showed generally low SHAPE reactivity yet high Shannon entropy. While the low SHAPE reactivity indicated that this region is structured, the high Shannon entropy suggested the presence of several possible conformations with unique SHAPE reactivities, which may or may not dynamically interconvert. The correlation coefficients comparing our ex cellulo and in vitro SHAPE profiles varied more between replicates at the 3′-end of SChLAP1 than at the 5′-end (Supplemental Fig. S7), suggesting structural variability between experiments in this region. Second, evidence of retroviral insertion (e.g., of THE1B) in this region indicated that the local RNA structures may have emerged from an external source; indeed, recent evidence points to the role of TEs in forming de novo RNA structural domains. In particular, Feschotte and coworkers compared the predicted structures of 100 lncRNAs with high TE content (over 96% of sequence) to 100 low-TE content lncRNAs (<3% of sequence) and observed a greater propensity for statistically significant secondary structures in high-TE lncRNA as evaluated by a sequence randomization approach similar to ScanFold (Kapusta et al. 2013). While THE1B insertions have been previously analyzed for their roles as gene promoters and enhancers for both coding and noncoding genes (Lamprecht et al. 2010; Young et al. 2013; Dunn-Fletcher et al. 2018), to our knowledge they have not been analyzed in the context of modifying RNA structure. Thus, this evidence for a structured RNA domain corresponding to the THE1B insertion, along with evidence that the THE1B insertion remains folded in cells while the LTR12C insertion is unfolded in cells, led us to hypothesize that this insertion could have resulted in the formation of a complex structural landscape leading to functional protein binding capability.

To assess the folding landscape of the 3′-end of SChLAP1, specifically within exon 7 and the THE1B insertion, we used native gel electrophoresis. This method was chosen due to its flexibility to accommodate RNA of various sizes and ability to assess varying conditions (e.g., salt identity and concentration) in an efficient manner (Woodson and Koculi 2009). Using our ex cellulo model along with additional structure predictions in RNAstructure (Bellaousov et al. 2013) as guidance for substructures to study, two strongly predicted structures from the 3′-end of SChLAP1 were chosen for analysis by native gel electrophoresis: Arm A (nucleotides 949–1116), which contains the above-mentioned lowSS AU-rich helix, and Arm B (nucleotides 1117–1428), which contains the 3WJ and corresponds to the THE1B insertion (Supplemental Fig. S10A).

We began by looking at two variables, magnesium (Mg²⁺) concentration and annealing method, and their impact on Arm A/B migration by gel electrophoresis. Mg²⁺ is a divalent cation known to be crucial for the stabilization of RNA tertiary structures (Butcher and Pyle 2011), whereas the method of RNA annealing, specifically whether the RNA is snap-cooled or slow-cooled, may modify the resultant conformational landscape of the RNA population (Somarowthu et al. 2015). In our initial examination of these variables, wherein Arm A/B were snap- or slow-cooled in filtration buffer and thereafter incubated with varying Mg²⁺ concentrations (0, 2.5, 5, 7.5, and 10 mM Mg²⁺; protocol designed for consistency with our chemical probing above), Arm A migrated at a consistent distance and as one conformer via native gel regardless of either variable (Fig. 5, top). This result indicated that neither the annealing protocol nor the Mg²⁺ concentration significantly impacts the resultant structure of Arm A. We note that this finding does not rule out more subtle conformational changes that may occur within Arm A. In contrast, for Arm B we observed an unexpected trend between Mg²⁺ concentration and annealing protocol (Fig. 5, bottom). When Arm B is snap-cooled, Mg²⁺ addition significantly alters the conformational space of Arm B, favoring the formation of two central conformations over a lower and higher conformation observed at low magnesium concentrations (e.g., 0 and 2.5 mM Mg²⁺). However, when Arm B is slow-cooled, this Mg²⁺ dependence is greatly reduced, and Arm B migrates in one major conformer throughout the magnesium titration. We note that slow-cooling of Arm B results in a minor population of more slowly migrating (higher molecular weight) band that is not observed in snap cooling, which we attribute to dimerization.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 5.

Native gel electrophoresis of Arm A and Arm B with varying annealing protocols and magnesium concentrations. RNA was snap- or slow-cooled in filtration buffer (50 mM K-HEPES pH 7.5, 150 mM KCl, 0.1 mM EDTA) before incubation with the given Mg²⁺ concentration. Asterisk denotes putative Arm B dimer favored in slow cooling. RNA was purified using the “kit-purified” approach. Arm A gel is representative of two biological replicates, and Arm B gel is representative of three biological replicates.

Next, we evaluated how snap- and slow-cooling impact the migration of Arm B when incubated with divalent cations besides Mg²⁺. We evaluated the migration of Arm B in response to 5 mM Ba²⁺, Sr²⁺, Ca²⁺, Mn²⁺, or Ni²⁺ (all as metal chlorides). We chose these divalent metals as their roles in RNA folding and stability have been assessed in previous studies, particularly for RNA G-quadruplexes and the Tetrahymena ribozyme, and they cover an array of chemical properties including radius, charge density, and hydration enthalpy (Koculi et al. 2007; Balaratnam and Basu 2015). As above, Arm B was snap- or slow-cooled in filtration buffer before incubation with the indicated divalent metal. For snap-cooled Arm B, we found that some metals, such as Ba²⁺ and Sr²⁺, resulted in Arm B migrating primarily in the lowest conformer (i.e., fastest moving band), previously observed with 0 mM Mg²⁺. However, when incubated with other metals, such as Mn²⁺ and Ni²⁺, the additional, higher apparent conformations of Arm B are preferred (Fig. 6A). Again, this trend was not observed for slow-cooled Arm B (Fig. 6B), where the lowest conformer is observed regardless of the identity of the divalent metal.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 6.

Native gel electrophoresis of Arm B with various divalent metals and annealing conditions. Arm B was prepared by snap-cooling (A), slow cooling (B), or native isolation (C) in filtration buffer (see Materials and Methods) before incubation with divalent metals. Asterisk denotes dimer favored during slow cooling. Black line through gels is for reference only; magnesium titration and incubation for other divalent metals was run on the same gel for each panel. RNA was generated using FPLC purification, and each gel was run in biological duplicate (kit-purified of snap- and slow-cooled Arm B run in triplicate, Supplemental Fig. S10B).

Given that changes in Arm B migration in the presence of divalent metal were dependent on snap- or slow-cooling, we next compared natively purified Arm B, where the cotranscriptionally folded structure is preserved. We isolated Arm B using the same native purification approach that we used for full-length SChLAP1 and thereafter incubated the RNA with the aforementioned divalent metals. We observed that the metals did not induce a significant change in migration for natively folded Arm B (Fig. 6C), that is, Arm B migrated as one major conformer regardless of the identity of the divalent metal. We additionally confirmed that the differences inherent to the native purification workflow did not impact the resultant RNA conformational landscape by snap- and slow-cooling Arm B that had previously been natively purified. The same trends were observed for these samples by native gel: snap-cooling Arm B that was previously natively purified resulted in the above-mentioned metal-dependence that was not observed if natively purified Arm B was slow-cooled (Fig. 6; Supplemental Fig. S10B).

Together, our native gel data suggested that metal recognition by Arm B is significantly altered by annealing method, wherein snap-cooled Arm B showed altered migration by divalent metals depending on their identity, but slow-cooled and natively folded Arm B showed largely consistent migration regardless of divalent metal. These results indicated that (1) a structural difference within Arm B occurs between snap- and slow-cooling, (2) slow-cooling rather than snap-cooling much more effectively mimics the native fold of Arm B, and (3) divalent metals alter the conformational landscape of a nonnative Arm B fold. These data give support for multiple discrete conformations within this region of SChLAP1 in vitro, which may in turn explain the elevated Shannon entropy observed in our ex cellulo SHAPE data (Fig 3C). As stated above, this complex structure additionally appears to function as a protein binding region.

THE1B insertion within SChLAP1 may form a G-quadruplex, but is not sufficient to explain its conformational heterogeneity in vitro

Mg²⁺ is generally thought to favor RNA compaction, but in our experiments, an increase in apparent size is observed with increasing Mg²⁺. The order of divalent cation identity in increasing the apparent size of snap-cooled Arm B is reminiscent of recent work of divalent metal recognition of G-quadruplexes, where metals such as Ni²⁺ induced an unfolded conformation of a G-quadruplex at lower concentrations over other metals such as Ba²⁺ (Balaratnam and Basu 2015). To evaluate potential G-quadruplex formation, we snap- and slow-cooled Arm B in a modified filtration buffer wherein the concentration of the G-quadruplex stabilizing cation K⁺ (in the form of KCl) was varied from 0 to 150 mM prior to separation by native gel. We observed a significant change in migration for both annealing conditions with increasing K⁺, indicating the presence of a K⁺-dependent structure regardless of annealing protocol (Fig. 7A). As a control, we studied the impact of K⁺ on Arm A migration, and we observed no change in migration in response to K⁺ for Arm A (Supplemental Fig. S11A). This result suggested that a G-quadruplex forms within Arm B in both annealing protocols and that G-quadruplex formation may be a component of the structural landscape in Arm B/THE1B.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 7.

Analysis of G-quadruplex formation of Arm B by native gel electrophoresis and SHAPE. (A) Native gel electrophoresis of kit-purified Arm B in response to varying KCl concentrations and annealing conditions. Gel is representative of two biological replicates. (B, top) Sequence of predicted G-quadruplex in Arm B by QGRS Mapper (Kikin et al. 2006). G-tracts are numbered, and guanosines predicted to participate in G-quadruplex are in bold. (Bottom) Violin plots comparing the reactivities of guanosines in predicted G-quadruplex to the remaining guanosines in each SHAPE probing experiment reported herein. Reactivities below −1 or above 2 were not plotted here but were included in statistics. A line at SHAPE = 0.4 is shown for reference as SHAPE reactivities above this value are unfavorable for base-pairing (Hajdin et al. 2013). Significance was calculated using a Kolmogorov–Smirnov test for each experiment. Nucleotides 1208 and 1209 are indicated by color coding.

To further investigate potential G-quadruplex formation, we used QGRS Mapper (Kikin et al. 2006) against the SChLAP1 Iso. 1 sequence to search for potential G-quadruplex forming regions. Searching with the default parameters, which allows G-quadruplex formation within a 30 nt window, we identified four lowly predicted potential G-quadruplex forming regions (Supplemental Fig. S11B, left). However, extending the maximum allowed G-quadruplex to 45 nt facilitated the prediction of a three-tract G-quadruplex (tract referring to the number of guanosines in a row per strand of the putative G-quadruplex) located within Arm B (nucleotides 1194–1228) (Supplemental Fig. S11B, right). This predicted G-quadruplex was the only predicted G-quadruplex with three G-tracts, which we hypothesized could confer added stability to this fold compared to the other predicted G-quadruplexes in SChLAP1, which only contained two G-tracts. In addition, this predicted G-quadruplex would have two large loops, one of which (nucleotides 1211–1225) was predicted by RNAstructure (Bellaousov et al. 2013) to form a stem–loop structure itself (Supplemental Fig. S11C).

This computational result supports our K⁺-dependent native gel studies, where evidence of G-quadruplex formation was observed in Arm B (and thus the THE1B insertion). However, we note that this G-quadruplex overlaps with the coordinates of the ScanFold/SHAPE-informed 3WJ (nucleotides 1141–1213) and adjacent stem–loop (nucleotides 1224–1242), indicating that a nonquadruplex conformation is also thermodynamically predicted.

Analysis of SHAPE data indicates a structured but nonquadruplex conformation of Arm B/THE1B insertion in cell-derived SChLAP1

Given the above indication of a G-quadruplex structure in Arm B of SChLAP1 in vitro, and its overlap with favorably predicted nonquadruplex secondary structures, we evaluated our in vitro and in/ex cellulo SHAPE data for evidence of G-quadruplex formation. Recent work from Schneekloth and coworkers using the probe 2-methylnicotinic acid imidazolide found that guanosines participating in a G-quadruplex within the NRAS mRNA 5′ UTR were generally low in SHAPE reactivity in vitro, supporting their structuredness (Balaratnam et al. 2023). We analyzed the SHAPE reactivity of the predicted G-quadruplex within Arm B identified above (nucleotides 1194–1228) to evaluate its structuredness and indications of quadruplex formation. In four of our five probing experiments, we observe a statistically significant decrease in SHAPE reactivity of the guanosines participating in the putative G-quadruplex compared to the SHAPE reactivity of all remaining guanosines (Fig. 7B), suggesting that the guanosines in this potential G-quadruplex are indeed functioning in an RNA structure that hinders SHAPE reactivity. While statistical significance was not obtained for in-cell replicate 1, the data trends in the same direction and may be consistent with heterogeneity of this structure in in cellulo RNA as indicated above. However, we observed that two residues in this predicted G-quadruplex (nucleotides 1208 and 1209, within G-tract 3) showed a large enhancement in SHAPE reactivity in our ex cellulo data as compared to our in vitro data (Fig. 7B). While our in vitro data showed low SHAPE reactivity in each G-tract, consistent with our K⁺-dependent native gels for Arm B that indicated G-quadruplex formation (Fig. 7A), the heightened SHAPE reactivity in our ex cellulo data for G-tract 3 in cell-derived SChLAP1 suggested nonquadruplex structure. This trend is also observed in the in cellulo data, although with greater variation between the replicates (Fig. 7B).

A nonquadruplex, 3WJ structure for this region is supported by our representative ex cellulo MFE structure (Fig. 4), which is in line with our ScanFold predictions (Fig. 2). In this model, the above-mentioned nucleotides 1208 and 1209 are base-paired across a 1 nt bulge, which may introduce more flexibility in these nucleotides despite their base-pairing and thus yield heightened SHAPE reactivity for these nucleotides in cell-derived RNA. Thus, while the data suggested that Arm B of SChLAP1, which encompasses the THE1B insertion, has a propensity to form G-quadruplexes, especially in vitro, this conformation appears less likely to form in cell-derived SChLAP1. Nonetheless, the combined data supports the SHAPE reactivity predictions of a structured region in the 3′-end of SChLAP1 with multiple discrete, potentially interconverting, structures.

Identification of protein binding partners for structured regions of SChLAP1

Based on our structural modeling and ΔSHAPE data, we selected four substructures of SChLAP1 for biotin pulldowns followed by mass spectrometry protein identification: the E2–E5 junction, the AU-rich helix, the 3WJ, and the CU-strand. As described above, the E2–E5 junction and the AU-rich helix were lowSS regions in our ex cellulo structure model (Fig. 3). While well-folded ex cellulo, the E2–E5 junction showed indications of in cellulo unfolding (Fig. 4; Supplemental Fig. S9), which we hypothesized could be due to protein-induced structural changes. The AU-rich helix (within Arm A) appeared to remain well-folded in cellulo as evidenced by a lack of ΔSHAPE in-cell enhancements, suggestive of a structurally stable site of RNA:protein interaction. The 3WJ (within Arm B) showed structural complexity/heterogeneity in our native gel work. While we attribute some of this architecture to G-quadruplex formation in vitro, comparison to our SHAPE data suggested that a nonquadruplex conformation was more likely for this region in cells (Fig. 7). Lastly, the CU-strand (within Arm B) is high in SHAPE reactivity in both our in vitro and ex cellulo models (Fig. 4; Supplemental Fig. S4A), indicating that it is likely single-stranded in a variety of conditions. The MFE structures of all RNA constructs used (in the absence of any SHAPE data) are shown in Supplemental Figure S12 and are largely consistent with the structures observed in our experimentally informed structural modeling of SChLAP1.

Protein pulldowns were performed for the above four constructs, which were ordered as biotinylated RNA oligonucleotides from IDT. The immobilized constructs were incubated with nuclear lysate derived from the LNCaP cell line. Silver staining of pulldown samples showed successful isolation of proteins by all constructs in comparison to a beads-only control (Fig. 8A). We submitted the eluates of the four RNA pulldown samples for analysis via qualitative proteomics. Approximately 5000 proteins were identified across all constructs in our experiment (Supplemental Data 3). We selected the 100 most abundant proteins across all four pulldowns for future analysis, which together constituted ∼70% of the signal intensity (Supplemental Fig. S13A). Gene ontology analysis revealed enrichment for functions related to RNA binding (Supplemental Fig. S13B), as expected. While we note that analysis of the highest signal proteins, particularly in a qualitative mass spectrometry experiment, may bias toward abundant proteins or limit contribution of transient interactions, we expected that this list would give a useful starting point for identifying previously undescribed SChLAP1:protein interactions with potential importance in prostate cancer.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 8.

Identification of SChLAP1 substructure binding proteins. (A) Silver stain of pulldown eluates. Gel is representative of two independent experiments. (B) Western blot for BRG1, HNRNPL, FUBP1, HNRNPK, and PCBP2 in pulldown samples. Western blots are representative of two independent experiments. (C) Comparison of previously identified binding motifs to ΔSHAPE sites identified in our work. Structures of E2–E5 junction, AU-rich helix, and 3WJ with ΔSHAPE sites indicated are reprinted from Figure 4. Purple highlights represent motifs for HNRNPL, FUBP1, or HNRNPK/PCBP2.

To identify putative structure-specific protein interactions for follow up validation by western blot (see Materials and Methods), we first normalized the abundance for a particular protein in each of the four pulldowns to a scale of 0–1 across the pulldowns. We heuristically defined putative structure-specific protein interactions as proteins with this normalized score ≤0.3 compared to the maximum pulldown (set to 1); that is, for a given protein across the four pulldowns, the three nonmaximum pulldowns had to have abundance values <30% of the maximum pulldown for a protein to be considered potentially SChLAP1 substructure-specific. Multiple proteins showed similar protein counts across multiple constructs (see Supplemental Fig. S13C for examples), which we attribute to nonspecific RNA:protein interactions. Out of the 100 top proteins, only 17 passed our selection criteria, having normalized scores for all three nonmaximum pulldowns of <0.3 (Supplemental Fig. S13C). We focused on proteins isolated by the E2–E5 junction, AU-rich helix, or 3WJ given the single-stranded nature of the CU-strand. The list of candidate proteins is shown in Supplemental Figure S13C, from which antibodies were selected.

Based on our selection criteria, the following proteins were evaluated for their presence in protein pulldown eluates by western blot: HNRNPL (most abundant for the E2–E5 junction), far upstream element binding protein 1 (FUBP1, most abundant for AU-rich helix), heterogeneous nuclear ribonucleoprotein K (HNRNPK, most abundant for 3WJ), and poly(rC) binding protein 2 (PCBP2, most abundant for 3WJ). In addition, we evaluated the binding of BRG1 to each of these constructs given the previous work on SChLAP1:SWI/SNF interactions (see Introduction).

Western blot and mass spectrometry analyses produced concordant findings. By western blot, HNRNPL was only observed to bind the E2–E5 pulldown sample (see Supplemental Text 2). Similarly, FUBP1 was only observed in the AU-helix pulldown, and HNRNPK and PCBP2 were only observed for the 3WJ pulldown (Fig. 8B). We additionally inspected each substructure for overlap with previously identified consensus sequences for these RNA-binding proteins (RBPs) (Fig. 8C). A CA-repeat within the E2–E5 junction overlaps with the ΔSHAPE site identified in our work (see Supplemental Text 2). While FUBP1 is generally considered to bind AU-rich sequences (Debaize and Troadec 2019), previous RNA Bind-n-Seq studies identified a more specific UNUUU motif for FUBP1:RNA interactions (Dominguez et al. 2018). Accordingly, we observed two concatenated UNUUU motifs overlapping the ΔSHAPE we identified in the AU-rich helix. Lastly, HNRNPK and PCBP2 are both polyC binding proteins (Makeyev and Liebhaber 2002). Five tracts containing three or more cytosines are observed in the 3WJ (Fig. 8C). In addition, analysis of RNA Bind-n-Seq as well as eCLIP data has suggested that HNRNPK preferentially binds RNA in secondary structure contexts (Dominguez et al. 2018; Lubelsky and Ulitsky 2018), which is in line with our findings of HNRNPK binding to the 3WJ. Of particular note, each of these proteins have themselves been implicated in prostate cancer progression (see Discussion), suggesting that some of their activity may come in part through their interaction with SChLAP1.

BRG1, while observed in the input sample (Fig. 8B), was not observed by western blot in any pulldown sample. Based on these data, we do not rule out BRG1 binding to these substructures in vivo. However, we believe that HNRNPL, FUBP1, HNRNPK, and PCBP2 constitute a set of preferred, selective protein interactors for these SChLAP1 substructures in vitro.