The human SNHG8, ZFHX4-AS1 and RP11-552M11.4 genes were cloned by reverse transcription and polymerase chain reaction (PCR) amplification methods. The purified DNA fragment from a gel (28004, MinElute PCR Purification Kit, Qiagen, MD, USA) was cloned into the pIRES-hrGFP-1a vector (240031, Agilent Technologies, Savage, MD, USA) using the hot fusion method by CE Design V1.04 (Vazyme Biotech Co., Ltd., Nanjing, China). Plasmids were introduced into Escherichia coli strain by chemical transformation. Then E. coli strains were plated onto LB agar plates with 1 µL/mL ampicillin (7177-48-2; Sigma, St. Louis, MO, USA), and incubated overnight at 37 °C. Colony PCR (28 cycles) was used to screen lncRNA colonies with universal primers (located on the backbone vector). Finally, the inserted sequences of lncRNAs were verified by Sanger sequencing.

Primary HTM cells (YS905C, Yaji Biological, Shanghai, China) were donated by anonymous deceased individuals without detailed information, cultured in DMEM/F12 (PM150310, Procell, Wuhan, China), and were validated by STR profiling, and tested negative for mycoplasma. Cells were cultured in a humidified incubator at 37 °C and 5% CO${}_2$. Plasmid transfection was performed with Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). For the in vitro stiffness model, HTM cells were cultured on polyacrylamide hydrogel with gradually increased bisacrylamide concentrations (0.04%, 0.1%, 0.2%, 0.5%, 1.3%, and 2.08%), mimicking the increased Young’s modulus (1.1, 2.5, 4.2, 11.9, 34.4, and 50 kPa) [15, 24]. For lncRNA plasmid transfection and enhancer of zeste homolog 2 (EZH2) immunoprecipitation experiments, HTM cells were treated under normal conditions.

We used glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as an internal control for assessing the effects of SNHG8, ZFHX4-AS1, and RP11-552M11.4 overexpression, as well as the ECM-related gene expression. The cDNA synthesis and following RT-qPCR were performed on the Bio-Rad S1000 Thermal Cycler with Hieff™ qPCR SYBR® Green Master Mix (Low Rox Plus; YEASEN, Shanghai, China), followed by denaturing at 95 °C for 10 min, 40 cycles of denaturing at 95 °C for 15 sec and annealing and extension at 60 °C for 1 min. The expression level of each transcript was normalized to GAPDH by the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method [25] using three technical replicates. The detailed sequences of primers are presented in Supplementary Table 1.

HTM cells were prepared for ultraviolet (UV) cross-linking with UV irradiation type C (254 nm) at 400 mJ per cm${}^2$. Cross-linked cell lysis was performed with 1% RNase inhibitor (2313u, Takara, Beijing, China) and a 2% protease inhibitor cocktail (B14001, Bimake, Shanghai, China) for 30 min. RQ1 DNase (M6101, Promega, Madison, WI, USA) was added to the lysate to digest the DNA. Then RNA was digested by adding micrococcal nuclease (88216, Thermo Fisher Scientific, Waltham, MA, USA). For IP, cell lysate was incubated with 15 µg anti-EZH2 antibody (5246, Cell Signaling Technology, Danvers, MA, USA) or control IgG antibody (12-371, Millipore, Burlington, MA, USA) overnight at 4 ℃ with two biological replicates. The immunoprecipitated protein-RNA complexes were eluted and isolated from the gel. Then the recovered RNA was prepared for library preparation using kit for small/microRNA (K02420-L, Gnomegen, San Diego, CA, USA). The libraries were applied to the Illumina HiSeq X Ten system for 150 nucleotide (nt) paired-end sequencing.

HTM cells were cross-linked with 1% formaldehyde for 10 min at room temperature and then lysed to obtain a chromatin-protein complex. Soluble sheared chromatin was obtained by sonication (average DNA length of 200–500 base pairs). The 10 µg anti-EZH2 antibody (5246, Cell Signaling Technology, Boston, MA, USA) and anti-H3K27me3 antibody (12-565, Millipore, Burlington, MA, USA) were used for IP with 100 µL chromatin as input at 4 °C overnight with two biological replicates. The next day, 30 µL protein beads were added and the samples were further incubated for 3 h. Immunoprecipitated DNA was used to construct the sequencing libraries following the protocol of ChIP-seq Library Prep Kit (NOVA-5143-02, Bioo Scientific, Austin, TX, USA) and the libraries were applied to the Illumina HiSeq X Ten system (Illumina, Inc., San Diego, CA, USA) for 150 nt paired-end sequencing.

The RNA-seq data files of GSE123100 from our previous study [15] were re-analyzed. The raw files were converted to fastq file by the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) Tool fastq-dump. Low-quality bases were removed from the raw reads using the FASTX-Toolkit (version 0.0.13, The Hannon Lab, New York, NY, USA; http://hannonlab.cshl.edu/fastx_toolkit/). The short reads less than 16 nt were also dropped. Then the quality filtered reads were evaluated using FastQC (version0.12.1, Babraham Institute, Cambridge, MA, USA, http://www.bioinformatics.babraham.ac.uk/projects/fastqc).

Raw sequenced reads were filtered to remove low quality reads using the above methods. Then, clean reads were aligned onto the human GRCH38 genome using the HISAT2 (version2.2.1, University of Texas Southwestern Medical Center, Dallas, TX, USA) alignment program [26] with no more than four mismatches. Uniquely mapped reads were used for gene expression calculation by FPKM values (fragments per kilobase of transcript per million fragments mapped) [27]. The statistical power of this experimental design (one biological replicates), calculated by online implementation (https://rodrigo-arcoverde.shinyapps.io/rnaseq_power_calc/) was 0.93.

StringTie (version2.1.6, Johns Hopkins University, Baltimore, MD, USA) [28] was used to assemble and predict transcripts. Final transcripts were screen out by eliminating the transcripts with FPKM 1, and then were combined into one general transfer format (GTF) file. We then used four software to predict coding potential of transcripts: CPC2 (version2.0, Peking University, Beijing, China) [29], LGC (version1.0, Chinese Academy of Sciences, Beijing, China) [30], CNCI (version2.0, Chinese Academy of Sciences, Beijing, China) [31] and CPAT (version3.4.0, Mayo Clinic College of Medicine, Houston, TX, USA) [32]. The filtering steps for lncRNA selection were followed one published study [33]. Briefly, Transcripts without coding potential, longer than 200 nt, and more than 1000 bp distance from the adjacent coding genes, were considered lncRNAs.

We used edgeR (version3.32, The Walter and Eliza Hall Institute of Medical Research, Parkville, Victoria, Australia) [34] to identify DEGs with p-value 0.01 and fold change $>$1.5 (up) or 0.67 (down) as the cut-off criteria according to our previous study using the same dataset [15]. Another reason for this threshold was the small amplitude of transcriptome alternation caused by the increased stiffness.

After alignment, uniquely mapped reads were used to predict EZH2-bound regions (peaks) by the “ABLIRC” strategy (version1.0, ABLife Inc., Wuhan, China) [35] and CIMS (version1.0, The Rockefeller University, New York, NY, USA) [36]. After stimulation and calculation, all the observed peaks with p-value 0.05 were selected. For ChIP-seq, reads were aligned to the reference using Bowtie2 (version2.4.4, University of Maryland, College Park, MD, USA) with one mismatch [37]. ChIP-seq peak calling was performed using MACS14 (version1.4.2, Dana-Farber Cancer Institute and Harvard School of Public Health, Boston, MA, USA) with default settings and p-value 1 $\times$ 10${}^{-5}$ [38]. The binding motifs of EZH2 protein were identified using HOMER software (version3.12, University of California, San Diego, CA, USA) [39] with the default setting.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified to sort out the functional categories of DEGs with the KOBAS server (version2.0, Peking University, Beijing, China) [40]. Principal component analysis (PCA) was performed using R package factoextra (https://cloud.r-project.org/package=factoextra). The pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R software (version4.2, https://www.r-project.org) was used to perform clustering based on Euclidean distance. Student’s t-test was used to compare the two groups with p-value 0.05 considered statistically significant. The Pearson’s correlation coefficient (PCC) was used to calculate the correlation between differentially expressed lncRNAs (DElncRs) and DEGs; p-value 0.01 and absolute PCC $>$0.9 were used as the criteria. Venn diagram and associated hypergeometric test were performed in R software (version 4.2).

Our previous study systematically investigated the mRNA expression profiles in response to the tunable stiffness and found extensive dysregulation of ECM-associated genes [15]. Here, we explored the expression profile of lncRNAs, which are important regulatory RNAs that may influence the global transcriptome profile. We re-analyzed the GSE123100 dataset and performed lncRNA prediction and expression analyses. The investigation pipeline used in this study is shown in Fig. 1A. Biological replicates were HTM cells cultured for 3 days (D3) and 5 days (D5) with the same Young’s modulus. After predicting novel lncRNAs and calculating their expression levels, we divided the expressed genes into mRNA and lncRNA groups. The PCA results demonstrated that D3 and D5 samples were separated by the first component, and that obvious distances for D5 samples with high Young’s modulus were observed (Fig. 1B). Then we assessed the expression patterns of several glaucoma-associated lncRNAs, including ANRIL/CDKN2B-AS1, MALAT1, and LOXL1-AS1 [18]. ANRIL/CDKN2B-AS1 showed lower expression levels at 1.1 kPa and 50 kPa; the expression level of LOXL1-AS1 was decreased at 50 kPa; whereas MALAT1 showed specifically higher expression level at 50 kPa (Fig. 1C). These results demonstrated that their expression patterns were dynamically changed with the increase of Young’s modulus. Together, these results indicated that in HTM cells, lncRNA expression levels are extensively regulated by increased substrate stiffness and that lncRNAs may have important regulatory functions. We then analyzed the differentially expressed lncRNAs (DElncRs) using the same method for mRNAs [15]. By analyzing DElncRs between adjacent samples with increased stress, we found that DElncR number was the largest in 34.3 kPa vs. 11.9 kPa, following by 11.9 kPa vs. 2.5 kPa, and that DelncR number was very small in samples with lowest or highest stress (Fig. 1D). According to a previous study, the average Young’s modulus of HTM was 80.0 $\pm$ 32.5 kPa in glaucomatous tissue [41], consistent with our discovery that the largest DElncR number was detected in the 34.3 kPa vs. 11.9 kPa group. The PCA result demonstrated the high discreteness of all of the DElncRs (Fig. 1E), suggesting that the expression patterns of lncRNAs vary with the increased stiffness.

Fig. 1.

Global analysis of long noncoding RNAs (lncRNAs) from human trabecular meshwork (HTM) cells. (A) The work flow of this study. (B) Principal component analysis (PCA) results showing the diversity of expressed lncRNAs in the HTM cell model. (C) Line plot showing the expression patterns of CDKN2B-AS1, MALAT1, and LOXL1-AS1 in HTM cells in response to stiffness. (D) Bar plot showing the numbers of differentially expressed lncRNAs (DElncRs) in different groups. (E) PCA result showing the diversity of DElncRs in the HTM cell model. mRNA, messenger RNA; CLIP, cross-linking immunoprecipitation; EZH2, enhancer of zeste homolog 2; DE, differentially expressed; Chip, chromatin immunoprecipitation; DEGs, differentially expressed genes; OE, over expression.

Based on the above results, we explored the potential functions of DElncRs in HTM cells. Focusing on the two groups with the most DElncRs (Fig. 1D), we found that 97 DElncRs were detected in both groups (Fig. 2A, p-value = 3.43347 $\times$ 10${}^{-50}$, hypergeometric test). Then we analyzed the potential functions of these shared DElncRs by co-expression analysis with the identified DEGs. After the co-expressed DElncR-DEG pairs were filtered with correlation coefficient $>$0.9 and p-value 0.01, 1621 and 4031 pairs were detected from the 11.9 kPa vs. 2.5 kPa and 34.4 kPa vs. 11.9 kPa groups, respectively. Then, we performed functional enrichment analyses of these co-expressed DEGs. From the 11.9 kPa vs. 2.5 kPa group, the immune response and cytokine-mediated signaling pathways were enriched (Fig. 2B). From the 34.4 kPa vs. 11.9 kPa group, co-expressed DEGs were significantly enriched in response to mechanical stimulus, negative regulation of chondrocyte differentiation, and integrin-mediated signaling pathways; ECM organization was also enriched (Fig. 2C). The cellular component of these co-expressed DEGs also revealed their enrichment in ECM-related pathways (Supplementary Fig. 1A,B). Fig. 2A shows the co-expressed DEGs of overlapped DElncRs. Although insignificant, the extracellular space and extracellular region were the top two cellular components (Supplementary Fig. 1C). These results indicate that DElncRs may regulate the expression of ECM-related genes and participate in the pathogenesis of glaucoma. To have a closer view of the DElncRs and their co-expressed DEGs, we plotted the potentially regulatory network between DElncRs and DEGs derived from ECM-related pathways (Fig. 2D). Heatmap analysis of DElncRs associated with ECM-related pathways revealed that they were highly dysregulated due to the increased substrate stiffness (Fig. 2E).

Fig. 2.

DElncRs are tightly associated with extracellular matrix (ECM) genes by co-expression analysis. (A) Venn diagram showing the overlapped DElncRs in the 11.9 kPa vs. 2.5 kPa and 34.4 kPa vs. 11.9 kPa group. (B) Bubble plot showing the top enriched Gene Ontology (GO) Biological Process (BP) pathways of co-expressed DElncRs-DEGs in the 11.9 kPa vs. 2.5 kPa group. (C) Bubble plot showing the top enriched GO BP pathways of co-expressed DElncRs-DEGs in the 34.4 kPa vs. 11.9 kPa group. (D) Network presentation of DElncR-DEG pairs enriched in ECM-related pathways. (E) Heatmap presentation of DElncRs co-expressed with ECM-related DEGs. DEG, Differentially Expressed Gene.

According to these results, we proposed that lncRNA expression levels were extensively changed between 11.9 kPa vs. 34.4 kPa, and were associated with ECM organization. Next, we selected several lncRNA-mRNA pairs with a negative correlation to further explore their relationships. The following 10 lncRNAs were selected to explore their potential functions: ZFHX4-AS1, NNF-AS1, CRNDE, TP73-AS1, SNHG8, DLEU2, CTD-3252C9.4, LINC00707, RP11-552M11.4, and MEG8. Most of these lncRNAs were located in intergenic regions, and may regulate other genes or RNAs in a cis or trans manner. Three lncRNAs, ZFHX4-AS1, NNF-AS1, and TP73-AS1, were at the antisense strand of protein coding genes, functioning in a cis manner. These known lncRNAs were classified into two clusters: increased expression levels by elevated stiffness, and the contrary cluster. The RT-qPCR experiment was performed to validate their expression patterns during substrate stiffness elevation. The RT-qPCR results of the eight lncRNAs showed a consistent pattern with the RNA-seq data (Fig. 3, Supplementary Table 1). These results demonstrated that the dysregulated lncRNAs may have important functions in the HTM cells of glaucoma pathogenesis.

Fig. 3.

Reverse Transcription and Quantitative Polymerase Chain Reaction (RT-qPCR) validation of lncRNAs under the effect of substrate stiffness in HTM cells. Box plot showing the RT-qPCR validation of selected lncRNAs for D3 and D5 samples. RNA-seq results (brown box) and RT-qPCR results (light green box) showed the consistency of these two methods. RT-qPCR analysis was performed with three biological replicates.

To further explore the functions of these selected lncRNAs, we attempted to construct lncRNA overexpression cell models by transfecting overexpression plasmids into normal HTM cells. Plasmid without lncRNA sequences was used as the common negative control (NC). RT-qPCR has validated the overexpressed transcriptional levels of several lncRNAs (Fig. 4A). Finally, three lncRNAs with significant expression changes were selected for overexpression experiments, including ZFHX4-AS1, SNHG8, and RP11-552M11.4. Meanwhile, ZFHX4-AS1 itself had control (NC_ZFHX4-AS1_OE) because the experiment was done in another batch, and SNHG8 and RP11-552M11.4 had the same control (NC_OE). Then RNA-seq experiments were performed to explore the transcriptional influence of these three lncRNAs. The expression analyses showed that they showed distinct higher levels in the corresponding overexpression samples, and were not changed in other samples (Fig. 4B). Sample correlation analysis revealed that overall gene expression levels were highly correlated with PCCs $\ge$0.98 except the NC_ZFHX4-AS1_OE sample (Fig. 4C). These results indicated that lncRNA dysregulation did not trigger overall transcriptional alternation. The EdgeR [34] method was employed to perform DEG analysis with false discovery rate (FDR) 0.05 and fold change (FC) $>$2. The DEG number between lncRNA overexpression vs. NC group ranged from 296 to 749 (Fig. 4D), accounting for 1–3% of total detected genes. We then performed functional enrichment analyses of these DEGs. Extracellular space- and extracellular region-related cellular component terms were significantly enriched in both the up- and down-regulated genes (Fig. 4E–G). The DEGs involved in ECM were regulated by these lncRNAs. ZFHX4-AS1 especially increased the expression levels of 24 ECM-related genes, and decreased the expression levels of 46 genes, thus having the most significant influence. Next, we investigated the global expression patterns of all ECM-related genes, and found that they were extensively regulated after treatment (Fig. 4H). Finally, we analyzed the overlapped DEGs among these three lncRNAs and found that most of the regulated DEGs were specific (Fig. 4I). We randomly selected three ECM-related genes, including APOE, EMILIN3, and F12, and validated their changed expression levels in SNHG8-OE and NC samples (Fig. 4J). Collectively, these results show that aberrant lncRNA expression can affect the ECM environment, which may contribute to the pathological condition of HTM cells.

Fig. 4.

LncRNA overexpression extensively regulates the expression of ECM remodeling-related genes. (A) Bar plot showing the overexpressed lncRNA levels in HTM cells using RT-qPCR method. (B) Bar plot showing the expression levels of three overexpressed lncRNAs and two negative controls. (C) Sample correlation heatmap of RNA-seq samples. (D) Volcano plot showing the DEGs between overexpressed lncRNAs and negative controls. (E) Bubble plot showing the enriched GO Cellular Component (CC) terms of DEGs from SNHG8-OE samples compared with negative control. (F) The same with (D) but for the DEGs from ZFHX4-AS1-OE. (G) The same with (D) but for the DEGs from RP11-552M11.4-OE. (H) Hierarchical clustering heatmap showing the expression patterns of ECM-associated and ZFHX4-AS1-regulated DEGs. (I) Venn diagram showing the overlapped and specific DEGs among these three DElncRs. (J) Bar plot showing the changed expression levels of APOE, EMILIN3, and F12 between SNHG8-OE vs. NC samples using RT-qPCR. One NC sample was removed due to its outlier feature. ** p-value 0.01; *** p-value 0.001; Student’s t-test. NC, negative control; OE, over expression.

LncRNAs regulate their targets in multiple functional types with the help of RNA binding proteins (RBPs), including cis- and trans-acting [42]. A well-studied partner of lncRNAs is EZH2, the core component of the PRC2 complex [43]. We performed CLIP-seq analyses for EZH2 to explore the interacting RNA partners of EZH2 in HTM cells. After alignment of sequenced tags, genomic region distribution analysis showed that compared to IgG, EZH2 preferred to bind to coding sequence (CDS) and 3${}^{\prime }$UTR regions (Fig. 5A). ABLIRC [35] and CIMS [44] tools were used to detect EZH2 binding peaks. Several known noncoding RNAs interacting with EZH2 were found to be significantly enriched in EZH2 immunoprecipitation (EZH2-IP) samples compared with IgG control samples (p-value 0.05), including XIST [23] with 3 binding peaks, and MEG3 [43] with 10 binding peaks (Fig. 5B). We then searched the binding profiles of the above-checked lncRNAs, and found that only SNHG8 was bound by EZH2 (Fig. 5C). SNHG8, a small nucleolar RNA host gene, is known to affect several gastric cancer-specific pathways and the target genes of Epstein–Barr virus (EBV) [45]. The above results showed that SNHG8 expression levels decreased under higher stiffness (Fig. 3); and that SNHG8 overexpression could also cause the expression of many ECM-related genes to be dysregulated (Fig. 4D). We analyzed the preferred binding sequences (motifs) of EZH2, and found that AGGUAAG and GAA-repeats were the two most significant motifs for both peak sets (Fig. 5D). We checked the enriched functional pathways of EZH2-bound transcripts in HTM cells. EZH2-bound transcripts identified by CIMS were significantly enriched in axon guidance, gene expression, and ECM organization (Fig. 5E). Overlapping analyses between EZH2-bound transcripts and SNHG8-OE-DEGs identified the overlapped genes (Fig. 5F,G). This indicates that SNHG8 and EZH2 may coordinately regulate the expression of ECM-related genes.

Fig. 5.

EZH2 interacts with mRNAs and lncRNAs in HTM cells. (A) Bar plot showing the genomic distribution of mapped reads in EZH2-IP and IgG samples. (B) Table showing the detailed information of EZH2-bound peaks from XIST and MEG3 lncRNAs using ABLIRC method. (C) Reads density plot for lncRNA SNHG8 in EZH2-IP and IgG samples. (D) Motif analysis showing the top three enriched motifs from EZH2-bound peaks using ABLIRC and CIMS methods. (E) Bubble plot showing the top enriched GO BP pathways of EZH2-bound gene transcripts using CIMS method. (F) Venn diagram showing the overlapped genes between EZH2-bound transcripts (replicate 1) and SNHG8-OE-DEGs. (G) Venn diagram showing the overlapped genes between EZH2-bound transcripts (replicate 2) and SNHG8-OE-DEGs.

The physical interactions between EZH2 and SNHG8 have been demonstrated, but it remains unclear how this protein-RNA complex plays a role in HTM cells remains unclear. Increasing evidence has shown that EZH2-RNA interactions prevent EZH2 from catalyzing H3K27 methylation to regulate gene expression [46, 47, 48]. Thus, we performed ChIP-seq analyses of EZH2 and H3K27me3, with two biological replicates for each. After MACS2 software was used to call peaks [38], significantly enriched binding peaks of EZH2 and H3K27me3 were obtained (p-value 1 $\times$ 10${}^{-5}$) (Fig. 6A). After analyzing the overlapped genes between two biological replicates, 1168 and 4251 genes were found to be bound by EZH2 and H3K27me3, respectively (Fig. 6B,C), with 266 genes being bound by both EZH2 and H3K27me3 (Fig. 6D). Functional enrichment analyses of these genes revealed diverse functions (Fig. 6E). The SNHG8-OE-regulated genes and EZH2- and H3K27me3-bound genes were overlapped to make the functional impacts of the interactions between EZH2 and SNHG8. A total of 49 SNHG8-OE-regulated genes were found to be bound by EZH2 or H3K27me3 (Fig. 6F). Functional analyses of these 49 genes showed enriched GO terms extracellular space and extracellular matrix (Fig. 6G), including HAPLN3, TACSTD2, COMP, TAC1, GDF10, DLK1, and EGF. These results demonstrate that SNHG8 and EZH2 could coordinately regulate the expressions of a subset of genes related to ECM-remodeling via the methyltransferase activity of EZH2. However, a large number of DEGs by SNHG8-OE were not bound by EZH2 or H3K2me3, indicating that other mechanisms could be involved in the SNHG8-regulated transcriptome in HTM cells, which needs to be further investigated in the future.

Fig. 6.

ChIP-seq analysis reveals the DNA interaction profiles of EZH2 and H3K27me3 in HTM cells. (A) Peaks and genes identified by EZH2 and H3K27me3 ChIP-seq data are shown in the bar plot. (B) Venn diagram showing the overlapped genes between the two replicates of EZH2 ChIP-seq- associated genes. (C) Venn diagram showing the overlapped genes between the two replicates of H3K27me3 ChIP-seq-associated genes. (D) Venn diagram showing the overlapped genes between EZH2 and H3K27me3 ChIP-seq- associated genes. (E) Bubble plot showing the enriched GO BP pathways of overlapped genes in (D). (F) Venn diagram showing the overlapped genes among SNHG8-OE-DEGs, EZH2-associated genes, and H3K27me3-associated genes. (G) Bubble plot showing the enriched GO CC pathways of 49 overlapped genes in (F).