Functional annotation of human long noncoding RNAs via molecular phenotyping

Abstract

Long noncoding RNAs (lncRNAs) constitute the majority of transcripts in the mammalian genomes, and yet, their functions remain largely unknown. As part of the FANTOM6 project, we systematically knocked down the expression of 285 lncRNAs in human dermal fibroblasts and quantified cellular growth, morphological changes, and transcriptomic responses using Capped Analysis of Gene Expression (CAGE). Antisense oligonucleotides targeting the same lncRNAs exhibited global concordance, and the molecular phenotype, measured by CAGE, recapitulated the observed cellular phenotypes while providing additional insights on the affected genes and pathways. Here, we disseminate the largest-to-date lncRNA knockdown data set with molecular phenotyping (over 1000 CAGE deep-sequencing libraries) for further exploration and highlight functional roles for ZNF213-AS1 and lnc-KHDC3L-2.

Figure 1.

Selection of lncRNA targets, their properties, and the study overview. (A) CAGE expression levels at log₂TPM (tags per million) and human dermal fibroblasts (HDFs) specificity of lncRNAs in the FANTOM CAT catalog (Hon et al. 2017) (N = 62,873; gray), lncRNAs expressed in HDFs (N = 6125; blue), and targeted lncRNAs (N = 285; red). The dashed vertical line indicates most lowly expressed lncRNA target (∼0.2 TPM). (B) Gene conservation levels of lncRNAs in the FANTOM CAT catalog (gray), lncRNAs expressed in HDFs (blue), and targeted lncRNAs (red). Crossbars indicate the median. No significant difference is observed when comparing targeted and expressed in HDF lncRNAs (Wilcoxon P = 0.11). (C) Similar to that in B but for genomic classes of lncRNAs. Most of the targeted lncRNAs and those expressed in HDFs are expressed from divergent promoters. (D) Subcellular localization (based on relative abundances from RNA-seq fractionation data) for targeted lncRNAs. Chromatin-bound (N = 98; blue); nuclear soluble (N = 76; green); cytoplasmic (N = 108; red). Black contours represent the distribution of all lncRNAs expressed in HDFs. (E) Example of ZNF213-AS1 loci showing transcript model, CAGE, and RNA-seq signal along with targeting ASOs. (F) Number of ASOs for target lncRNAs and controls used in the experiment. (G) Schematics of the study.

Figure 2.

Cell growth and morphology assessment. (A) Selected example (PTPRG1-AS1) showing the normalized growth rate estimation using a matching NC_A (negative control). (B) Correlation of the normalized growth rate for technical duplicates across 2456 Incucyte samples. (C) Density distribution of normalized growth rates (technical replicates averaged) 252 ASOs targeting lncRNAs with successful knockdown (KD) and growth phenotype (blue) consistent in two replicates (FDR < 0.05 as compared to matching NC_A; 246 ASOs inhibited growth), 627 ASOs targeting lncRNAs with successful KD (purple), 270 negative control (NC_A) samples (gray), and 90 mock-transfected cells (Lipofectamine only) samples (yellow). (D) MKI67 staining (growth inhibition validation) for four selected lncRNA targets after siRNA and ASOs suppression. (E) Incucyte cell images of selected distinct cell morphologies changes upon an lncRNA KD. (F) An overview of the cell morphology imaging processing pipeline using a novel lncRNA target, CATG000089639.1, as an example. (G) lncRNAs (N = 59) significantly (FDR < 0.05) and consistently (after adjusting for the number of successfully targeting ASOs) affecting cell growth (N = 15) and cell morphologies (N = 44).

Figure 3.

CAGE predicts cellular phenotypes. (A) RT-qPCR knockdown efficiency for 2021 ASO-transfected samples (targeted lncRNAs only). Gray dashed line indicates 50% KD efficiency generally required for CAGE selection. Purple dashed lines indicate median KD efficiency (71.5%) for 375 ASOs selected for CAGE sequencing. After quality control, 340 ASOs targeting lncRNAs were included for further analysis. (B) Distribution of significantly differentially expressed genes (up-regulated: FDR < 0.05, Z-score > 1.645, log₂FC > 0.5; and down-regulated: FDR < 0.05, Z-score < −1.645, log₂FC < −0.5) across all 340 ASOs. (C) Motif Response Activity Analysis (MARA) across 340 ASOs. Scale indicates Z-score of the relative motif activity (the range was set to abs[Z-score] = <5 for visualization purposes). (D) Correlation between normalized growth rate and motif activities across 340 ASOs targeting lncRNAs with highlighted examples. Motif sizes shown are scaled based on the HDF expression of their associated TFs (range 1 to ∼600 TPM). (E) Enriched biological pathways across 340 ASOs. Scale indicates GSEA enrichment value calculated as −log₁₀(p) × sign(NES). (F) Same as in D but for selected GSEA pathways. Pathways sizes are scaled based on the number of associated genes.

Figure 4.

ZNF213-AS1 regulates cell growth, migration, and proliferation. (A) Normalized growth rate across four distinct ASOs (in duplicate) targeting ZNF213-AS1 as compared to six negative control samples (shown in gray). (B) Enrichment of biological pathways associated with growth, proliferation, wound healing, migration, and adhesion for ASO_02 and ASO_05. (C) Most consistently down- and up-regulated transcription factor binding motifs including those for transcription factors known to modulate growth, migration, and proliferation such as for example EGR family, EP300, GTF2I. (D) Knockdown efficiency measured by RT-qPCR after wound closure assay (72 h posttransfection) showing sustained suppression (65%–90%) of ZNF213-AS1. (E) Transfected, replated, and mitomycin C (5 µg/mL)-treated HDF cells were scratched and monitored in the Incucyte imaging system. Relative wound closure rate calculated during the 24 h postscratching shows 40%–45% reduction for the two targeting ASOs (ASO_02 [N = 10] and ASO_05 [N = 13]) as compared to NC_A transfection controls (N = 33, shown in gray) and the representative images of wound closure assay 16 h postscratching.

Figure 5.

RP11-398K22.12 down-regulates KCNQ5 and CATG00000088862.1 in cis. (A) Changes in expression levels of detectable genes in the same topologically associated domain (TAD) as RP11-398K22.12 based on Hi-C analysis. Both KCNQ5 and CATG00000088862.1 are down-regulated (P < 0.05) upon the knockdown of RP11-398K22.12 by two independent ASOs in CAGE analysis (left) as further confirmed with RT-qPCR (right). (B) (Top) Representation of the chromatin conformation in the 4-Mb region proximal to the TAD containing RP11-398K22.12, followed by the locus gene annotation, CAGE, RNA-seq, and ATAC-seq data for native HDFs. (Bottom) Schematic diagram showing Hi-C predicted contacts of RP11-398K22.12 (blue) and KCNQ5 (gray) (25-kb resolution, frequency ≥ 5) in HDF cells. Red line indicates RP11-398K22.12 and KCNQ5 contact. (C) FISH image for RP11-398K22.12, suggesting proximal regulation. TUG1 FISH image (suggesting trans regulation) is included as a comparison; (bar = 10 µm). (D) GTEx atlas across 54 tissues (N = 9662 samples) shows relatively high expression levels of RP11-398K22.12 in 13 distinct brain regions samples (highlighted). (E) Expression correlation for RP11-398K22.12 and KCNQ5 in eight out of 13 distinct brain regions, as highlighted in D. (F) Expression correlation for RP11-398K22.12 and CATG00000088862.1 in eight out of 13 distinct brain regions, as highlighted in D.