Expression of mRNAs, miRNAs, and lncRNAs in Human Trabecular Meshwork Cells Upon Mechanical Stretch.

Purpose: Intraocular pressure (IOP), the primary risk factor for primary open-angle glaucoma, is determined by resistance to aqueous outflow through the trabecular meshwork (TM). IOP homeostasis relies on TM responses to mechanical stretch. To model the effects of elevated IOP on the TM, this study sought to identify coding and non-coding RNAs differentially expressed in response to mechanical stretch.
Methods: Monolayers of TM cells from non-glaucomatous donors (n = 5) were cultured in the presence or absence of 15% mechanical stretch, 1 cycle/second, for 24 hours using a computer-controlled Flexcell unit. We profiled mRNAs and lncRNAs with stranded total RNA sequencing and microRNA (miRNA) expression with NanoString-based miRNA assays. We used two-tailed paired t-tests for mRNAs and long non-coding RNAs (lncRNAs) and the Bioconductor limma package for miRNAs. Gene ontology and pathway analyses were performed with WebGestalt. miRNA-mRNA interactions were identified using Ingenuity Pathway Analysis Integrative miRNA Target Finder software. Validation of differential expression was conducted using droplet digital PCR.
Results: We identified 219 mRNAs, 42 miRNAs, and 387 lncRNAs with differential expression in TM cells upon cyclic mechanical stretch. Pathway analysis indicated significant enrichment of genes involved in steroid biosynthesis, glycerolipid metabolism, and extracellular matrix-receptor interaction. We also identified several miRNA master regulators (miR-125a-5p, miR-30a-5p, and miR-1275) that regulate several mechanoresponsive genes. Conclusions: To our knowledge, this is the first demonstration of the differential expression of coding and non-coding RNAs in a single set of cells subjected to cyclic mechanical stretch. Our results validate previously identified, as well as novel, genes and pathways.

Glaucoma is a group of optic neuropathies characterized by progressive loss of retinal ganglion cells (RGCs), optic nerve atrophy, and visual field loss., Glaucoma affects more than 70 million people worldwide, with primary open-angle glaucoma (POAG) being the most common subtype.– Often, POAG remains undiagnosed until visual field loss is clinically severe. Known risk factors for POAG include advanced age, a positive family history of glaucoma, African or Hispanic ancestry, and/or elevated intraocular pressure (IOP)., Elevated IOP is the only clinically modifiable risk factor. Regardless of their starting IOP, lowering IOP in glaucoma patients, with pharmacological interventions or surgeries can delay the progression of vision loss.

IOP is determined by the dynamic production and outflow of the aqueous humor (AH)., The AH is secreted from the ciliary epithelium into the posterior chamber, travels through the pupil, and exits the anterior chamber primarily through the conventional pathway of the trabecular meshwork (TM) and Schlemm's canal, with a small portion draining via the unconventional pathway., The resistance to unimpeded outflow determines IOP. Excessive resistance to AH outflow through the TM causes elevated IOP, which may lead to compression and damage of RGC axons at the region of the lamina cribrosa. This elevated pressure will result in progressive peripheral vision loss and, if not treated, may result in complete, irreversible blindness. Even though conventional outflow pathway dysfunction is responsible for elevated IOP, most IOP-lowering medications target secretory processes of the unconventional outflow pathway. New medications targeting the conventional outflow pathway, such as the novel therapeutic Rhopressa (Aerie Pharmaceuticals, Durham, NC, USA), are in high demand.

Due to large fluctuations in IOP from blinking, eye movement, and ocular pulse, TM cells are constantly under mechanical stretch.,, Strain causes profound changes to cell morphology, affecting motility, stiffness, contraction, orientation, and cell alignment.– TM cells must react to this stress in order to prevent injury., Recent studies have indicated autophagy as one of the relevant stretch response adaptive mechanisms.,, Other mechanoresponsive genes include those involved in extracellular matrix (ECM) synthesis/remodeling, cytoskeletal organization, and cell adhesion.,– It remains unknown how these genes are regulated. Because microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) regulate gene expression, we aim here to identify the RNA profile changes in response to mechanical strain in TM cells and determine how coding and non-coding RNAs interact in stretched TM cells.

Our study focuses on discovering specific pathways that may induce elevated IOP. We hypothesize that stretch-responsive pathways in human TM are critical for maintaining AH outflow resistance homeostasis and thus modulating IOP. The identification of such pathways will enable the identification of novel therapeutic targets. To identify these genes and pathways, we stretched cultures of human TM (HTM) cells and measured genome-wide mRNA, miRNA, and lncRNA expression profiles and their respective signaling pathways. To our knowledge, this is the first study to examine the stretch-responsive differential expression of both coding and non-coding RNAs in the same set of cells.

Methods

Cell Culture

Primary cultures of HTM cells were obtained from cadaver eyes without a history of eye disease (Table 1)., Tissues were processed in accordance with the tenets of the Declaration of Helsinki. The study protocol was approved by the Institutional Review Board of Duke University Medical Center. Secondary cultures of cells expanded from primary isolates of HTM cells were grown at 37°C in 5% CO2 in low glucose Dulbecco's modified Eagle medium (DMEM) with l-glutamine, 110 mg/ml sodium pyruvate, 10% Fetal Bovine Serum – Premium Select (Atlanta Biologicals, Flowery Branch, GA, USA), 100 µM non-essential amino acids, 100 units/ml penicillin, and 100 µg/ml streptomycin from Invitrogen (Thermo Fisher Scientific, Waltham, MA, USA), as previously described. TM cells were characterized using established standards.

Table 1.

HTM Cells Derived from Postmortem Donors Without a History of Eye Disease (n = 5)

Donor	Age	Gender	Ethnicity	Cause of Death
TM93	54 y	Male	Unknown	Unknown
TM122	35 y	Male	Unknown	Unknown
TM126	88 y	Female	Unknown	Unknown
TM136	3 mo	Female	Caucasian	Chronic lung disease
TM141	38 y	Female	Caucasian	Respiratory

HTM cells from five non-glaucoma donors (Table 1) were plated on collagen-coated flexible silicone bottom plates (Flexcell International Corporation, Burlington, NC, USA). Cells from donor TM126 were studied twice in separate passages (TM126-1 and TM126-2). Results from these two passages were averaged for mRNA and lncRNA analyses. Due to batch effects, only the results for the first passage were included in the miRNA analysis. After reaching confluence, TM cells were switched to serum-free DMEM for 3 hours, followed by cyclic mechanical stretch for 24 hours (15% stretching, 1 cycle/second) using the computer-controlled FX-5000 Tension System (Flexcell International Corporation). A 24-hour time point was selected based on previous mechanotransduction studies.,, Control cells were cultured on flexible plates, under the same conditions, but without mechanical stretch. Stretch and RNA isolation protocols were conducted in two batches (batch 1: TM126-2, TM136, and TM141; batch 2: TM93, TM122, and TM126-1).

Gene Expression and Analysis

We extracted total RNA from the TM cells using the mirVana miRNA Isolation Kit with phenol (Thermo Fisher Scientific) following the recommended procedures as previously described. We evaluated the RNA quality using a 2100 Bioanalyzer with RNA 6000 Pico Kit (Agilent, Santa Clara, CA, USA). Only samples with an RNA Integrity Number ≥ 6 were used for RNA sequencing (RNA-Seq). A total of 200 ng RNA per sample was used to generate the sequencing libraries as previously described using the RiboGone – Mammalian kit and TaKaRa SMARTer Stranded RNASeq Kit (TaKaRa Bio USA, Inc., Mountain View, CA, USA), followed by sequencing with a NextSeq 500 System (Illumina, San Diego, CA, USA) using High Output v2 with paired-end 75-bp reads (Integrated Genomics Shared Resource, Georgia Cancer Center, Augusta University, Augusta, GA, USA). After quality checks and quality control, all the sequencing reads were demultiplexed and aligned using TopHat with paired-end reads. Cufflinks software was used to normalize sequencing read counts to fragments per kilobases and millions reads (FPKM).

For both mRNA and lncRNA data, the average expression of transcripts was determined for each of the two passages of TM126. Transcripts were excluded if they had six or more samples with FPKM < 0.001. In this way, only transcripts present at detectable levels in at least five (i.e., half) of the samples were retained for further analysis. To allow for expression analysis between stretched and control cells, missing data (i.e., values < 0.001) for any sample were replaced with a value of 0.001. Averages for all control and stretched samples were calculated, and those genes with an average of <1 for both control and stretched samples were removed in order to compare only genes with meaningful fold changes (FCs). Two-tailed paired Student's t-tests were performed and FCs calculated. Transcripts were considered to have significant differential expression if they had an absolute FC > 2 and an unadjusted P < 0.05.

For mRNAs, to focus on protein-coding RNAs we removed miRNA precursors and small nucleolar RNAs. In order to identify enriched functions and pathways, we conducted functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis on the annotated list of differentially expressed, protein-coding mRNAs using the 2013 version of the WEB-based GEne Set AnaLysis Toolkit (WebGestalt; available in the public domain, www.webgestalt.org). In order to annotate lncRNAs and identify their sequence similarity to coding genes, their identifiers were entered into the NONCODE version 5.0 database and the resulting sequence searched using the National Center for Biotechnology Information (NCBI) nucleotide-nucleotide Basic Local Alignment Search Tool (BLASTn).

miRNA Expression and Analysis

As previously described, the nCounter Human v3 miRNA Expression Assay Kit (NanoString Technologies, Seattle, WA, USA) was used on 100 ng total RNA to measure the expression of 800 human miRNAs that were selected prior to the analysis. This technique selectively examines the expression of the miRNAs assigned by the assay kit. After specific sample preparation and overnight hybridization, digital readouts of the relative miRNA abundance were obtained and translated to miRNA expression.

The raw data Reporter Code Count files produced by nCounter were processed with the NanoString nSolver 3.0 Software using a previously described analysis pipeline. Each sample was run through several quality control checks before analysis. These checks examined the quality of imaging, binding density, positive control linearity, and positive control limit of detection. Based on RNA content, the raw miRNA counts were normalized using the trimmed geometric mean with the nSolver software. The geometric mean calculated the average number of counts for each sample using the miRNAs with the median 40% of the counts. To account for technology-associated sources of variation, the miRNA counts were normalized further using the geometric mean of the positive control probes.

After normalization, the data were exported from nSolver into comma-separated values files and imported into the R Language environment for statistical computing. Using NanoString's recommended methods, the background level probes for control and stretched samples were identified within R. The negative control probes for each sample were extracted and multiplied by the corresponding normalization factors produced by nSolver. A Welch's t-test was performed for each probe to compare the normalized negative control counts to the sample counts. If the counts for a probe did not differ significantly from the negative controls (P > 0.05), then it was considered background. Finally, the BioConductor limma package was used to perform the differential analyses with paired samples after the background level probes were identified for each sample type. Due to batch effects, instead of averaging the transcript expression of the two passages of TM126 as was done for the mRNA and lncRNA analyses, the miRNA analysis included the expression of only one passage of TM126 as a representative for that donor cell line. Because many miRNAs demonstrated high significance at low fold changes, miRNAs with an absolute FC > 1.3, P < 0.05, and a change in counts greater than five were considered to be significantly differentially expressed. Cellular functions were determined by Ingenuity Pathway Analysis (IPA; Qiagen, Hilden, Germany).

Because the laboratory culture condition may change the miRNA expression in the cultures of HTM cells isolated from TM tissues, we examined whether these differentially expressed miRNAs were present in HTM tissues using our published miRNA expression profile from seven non-glaucoma HTM tissues with miRNA-Seq.

Integrative Analysis of miRNA and mRNAs

miRNAs are known to regulate the expression of many target genes. One miRNA may target many genes, and one gene could be regulated by many miRNAs. To identify miRNAs that potentially regulate the genes differentially expressed during the cyclic mechanical stretch of HTM cells, we used the miRNA Target Filter function in IPA. Briefly, we uploaded the list of differentially expressed miRNAs and the list of differentially expressed mRNAs to IPA. The miRNA Target Filter in IPA detects target genes from multiple sources (e.g., TargetScan, TarBase, miRecords, and Ingenuity Knowledge Base)– and identifies which genes in the uploaded mRNA gene list are potential targets of the selected miRNAs. We required that the expression of interested miRNAs be negatively paired with their potential mRNA targets. For accuracy, we limited our miRNA targets to those experimentally validated or predicted at high confidence.

Droplet Digital PCR Validation of Differentially Expressed mRNAs and miRNAs

Approximately 100 ng total RNA was used for mRNA reverse transcription using High-Capacity cDNA Reverse Transcription Kits (Applied Biosystems, Foster City, CA, USA). Targeted miRNAs were reverse transcribed using 1 to 10 ng total RNA, sequence-specific TaqMan miRNA RT primers (Supplementary Table S1), and the TaqMan MicroRNA Reverse Transcription Kit from Applied Biosystems. We used a QX200 droplet digital PCR (ddPCR) system (Bio-Rad, Hercules, CA, USA) and predesigned Bio-Rad ddPCR expression EvaGreen/Probe assays (Supplementary Table S2) to validate the differentially expressed protein-coding genes. Sequence-specific TaqMan miRNA probes (Supplementary Table S1) were used for differentially expressed miRNA validation. The absolute number of cDNA copies derived from ddPCR was normalized to the expression of the reference gene GAPDH. The fold change of normalized gene expression between mechanically stretched HTM cells versus control HTM cells was analyzed for each gene and miRNA using paired two-tailed Student's t-tests; P < 0.05 was considered significant.

Availability of Data and Materials

The miRNA dataset supporting the conclusions of this article is available in the NCBI Gene Expression Omnibus and are accessible through the GEO Series accession number GSE113755 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE113755).

Results

mRNA Differential Expression Analysis

All statistical tests were conducted in a pairwise fashion in order to take into account individual differences among cell strains. Furthermore, this approach helped eliminate any arbitrary batch effect differences. From gene differential analysis, using a cutoff with an absolute FC > 2 and a non-adjusted P < 0.05, we identified a total of 219 unique protein-coding mRNAs with significant differential expression (Supplementary Table S3). These genes included a couple of genes related to POAG, including UCP2 and TXNIP. These genes have functions in oxidative stress response (i.e., UCP2 and TXNIP) and inflammatory response (i.e., TXNIP). The top 20 downregulated and top 20 upregulated genes are listed in Table 2. To identify enriched functional categories, the list of these 219 genes was uploaded to WebGestalt. Many of these genes were identified as being involved in biological processes, including myeloid cell homeostasis, calcium ion sequestering, homeostatic regulation, cell-cycle regulation, glutamate secretion, rhythmic processes, mesenchyme morphogenesis, limb development, mammary gland epithelial cell proliferation, and sterol and lipid metabolism (Fig. 1, Supplementary Fig. S1). These genes are also involved in several molecular functions, including nucleoside binding, ribonucleoprotein complex binding, small molecule binding, hydrolase activity, nucleoside-triphosphatase activity, and nuclease activity (Fig. 1, Supplementary Fig. S2). Furthermore, the gene products were localized primarily to the following cellular components: cytoplasm, mitochondria, and other intracellular membrane-bounded organelles (Fig. 1, Supplementary Fig. S3). Pathway analysis based on the KEGG Pathway database indicated that the following KEGG pathways were enriched in response to cyclic mechanic stretch: metabolic pathways, steroid biosynthesis, arginine and proline metabolism, glycerolipid metabolism, colorectal cancer, phosphatidylinositol signaling system, small cell lung cancer, pathways in cancer, protein export, and ECM–receptor interaction (Supplementary Table S4).

Table 2.

Top 20 Upregulated and Top 20 Downregulated Protein-Coding Genes (P < 0.05) in Primary Human TM Cells in Response to Cyclic Mechanical Stretch (15%, 1 cycle/s, 24 h)

Gene Symbol	Gene Name	Fold Change	P
TEFM	Transcription elongation factor, mitochondrial	2187.24	3.45E-02
IMMP1L	Mitochondrial inner membrane protease subunit 1	1778.30	4.05E-03
UBE2C	Ubiquitin conjugating enzyme E2 C	1693.36	3.12E-02
TUBD1	Tubulin, δ1	1531.26	2.40E-02
RANBP3L	RAN binding protein 3 like	1427.57	4.11E-02
ADPRM	Manganese-dependent ADP-ribose/CDP-alcohol diphosphatase	1241.30	2.80E-02
TMEM81	Transmembrane protein 81	1220.01	1.56E-02
KDELC1	Protein O-glucosyltransferase 2	968.60	3.61E-02
DOK7	Protein Dok-7	876.85	9.90E-03
PTTG1	Securin	875.88	2.01E-02
LTB4R2	Leukotriene B4 receptor 2	749.49	2.81E-02
NME6	Nucleoside diphosphate kinase 6	745.65	3.30E-02
GNA14	Guanine nucleotide-binding protein subunit α-14	679.42	1.79E-02
DENND6B	Protein DENND6B	611.05	3.25E-02
RAD51C	DNA repair RAD51 homolog 3	550.61	3.09E-02
MYO1D	Unconventional myosin-Id	547.27	2.93E-02
TSSK2	Testis-specific serine/threonine-protein kinase 2	487.53	1.32E-02
GEMIN2	Gem-associated protein 2	404.78	3.02E-02
FAM69B	Divergent protein kinase domain 1B	318.25	4.84E-03
CLCN2	Chloride voltage-gated channel 2	500.71	5.36E-03
RBM3	RNA-binding motif protein 3	–2.78	1.44E-03
RPS6KA5	Ribosomal protein S6 kinase α-5	–2.79	5.00E-02
TTC30B	Tetratricopeptide repeat protein 30B	–2.79	4.22E-02
OSGIN1	Oxidative stress-induced growth inhibitor 1	–2.84	3.04E-02
SGCA	Alpha-sarcoglycan	–2.85	4.88E-02
ATHL1	Protein-glucosylgalactosylhydroxylysine glucosidase	–2.87	9.60E-03
MAML3	Mastermind-like protein 3	–2.88	5.09E-03
ARHGAP19	Rho GTPase-activating protein 19	–3.01	2.44E-03
FBXO4	F-box only protein 4	–3.05	1.65E-02
DUSP5	Dual-specificity protein phosphatase 5	–3.09	2.16E-02
HIST1H4E	Histone H4	–3.17	2.74E-03
MTSS1	Protein MTSS 1	–3.19	2.07E-03
HIST1H2AB	Histone H2A type 1-B/E	–3.25	4.13E-02
KIAA1407	Coiled-coil domain-containing protein 191	–3.39	1.53E-02
CRYL1	Lambda-crystallin homolog	–3.54	2.92E-03
TMEM91	Transmembrane protein 91	–3.61	3.24E-02
TXNIP	Thioredoxin-interacting protein	–5.53	4.54E-03
CXCL8	Interleukin-8	–7.44	4.80E-02
TCEB3CL	Elongin-A3 member B	–7.87	3.13E-02
RBM3	RNA-binding protein 3	–2.78	1.44E-03

Figure 1.

Top 10 gene ontology terms and pathways from WebGestalt functional and KEGG pathway analysis of the 219 significantly differentially expressed protein-coding genes (P < 0.05).

lncRNA Differential Expression Analysis

Following a paired two-tailed Student's t-test of lncRNA data, we identified a total of 387 unique lncRNAs with significant differential expression using a cutoff of an absolute FC > 2 and a non-adjusted P < 0.05 (Supplementary Table S5). The top 20 significantly downregulated and top 20 significantly upregulated lncRNAs are listed in Table 3. Several differentially expressed lncRNAs contained sequences similar to genes related to POAG (i.e., FOXC1 and OPTN) or to genes for POAG-related proteins (i.e., TGF-β receptor 2 and TGF-β receptor associated protein 1) (Table 4). Other lncRNA sequences shared similarity to genes related to steroid metabolism (i.e., G protein-coupled estrogen receptor 1, oxysterol binding protein like 10, ergosterol biosynthesis 28 homolog, and 24-dehydrocholesterol reductase), cell adhesion (i.e., neural cell adhesion molecule 1, protocadherin 10, and protocadherin related 15), and inflammation (i.e., VEGF-B, IL-1 receptor type 1, IL-17D, and nuclear factor-κB activating protein) (Table 4).

Table 3.

Top 20 Upregulated and Top 20 Downregulated lncRNAs (P < 0.05) in Primary Human TM Cells in Response to Cyclic Mechanical Stretch (15%, 1 cycle/s, 24 h)

Gene Name	Gene Homology	Fold Change	P
NONHSAT127145.2:776-856	BAG cochaperone 4 pseudogene	2.55E+07	2.8E-02
NONHSAT148356.1:2-154	Peroxisomal biogenesis factor 14	4.16E+06	1.9E-02
NONHSAT101532.2:162-254	Core 1 synthase, glycoprotein-N-acetylgalactosamine 3-β-galactosyltransferase 1	3.90E+06	3.5E-02
NONHSAT113163.2:447-625	High mobility group box 1 pseudogene 20	2.88E+06	2.4E-02
NONHSAT108452.2:8-101	Zinc finger and SCAN domain containing 26	2.85E+06	5.8E-03
NONHSAT218524.1:84-213	ADP-ribosylation factor 4	1.85E+06	4.2E-02
NONHSAT031917.2:83-1282	Solute carrier family 8 member A3	1.52E+06	4.9E-02
NONHSAT204812.1:1513-1657	Fibroblast growth factor 1	7.76E+05	4.9E-02
NONHSAT008749.2:42-248	Pleckstrin homology like domain family A member 3	6.46E+05	4.2E-02
NONHSAT209288.1:1114-1228	Lysophospholipase I	6.42E+05	1.2E-02
NONHSAT050260.2:1067-1189	Family with sequence similarity 149 member B1 pseudogene 1	6.20E+05	2.2E-02
NONHSAT023881.2:34-216	Transmembrane protein 123	5.07E+05	3.2E-02
NONHSAT200405.1:6-227	Protocadherin 10	4.50E+05	2.3E-02
NONHSAT194954.1:286-398	H3 histone, family 3B	4.26E+05	1.2E-02
NONHSAT075289.2:0-305	Sodium voltage-gated channel α subunit 9	3.47E+05	2.0E-02
NONHSAT202847.1:15-187	Chromosome 5 clone RP11-265O6	3.21E+05	3.4E-02
NONHSAT200755.1:22-221	Teneurin transmembrane protein 3	3.01E+05	5.0E-02
NONHSAT072212.2:21-297	Ankyrin repeat domain 36	2.95E+05	2.7E-02
NONHSAT137367.2:140-443	Akirin 1 pseudogene 2	2.66E+05	2.1E-02
NONHSAT196674.1:7-171	Insulin-like growth factor 2 mRNA binding protein 2	2.45E+05	2.1E-02
NONHSAT126268.2:21-577	Adaptor-related protein complex 3 subunit µ2	–3.91	1.6E-03
NONHSAT222566.1:0-570	Spermidine/spermine N1-acetyltransferase 1	–3.95	2.4E-03
NONHSAT069848.2:940-1450	Protein phosphatase 1 catalytic subunit β	–3.98	2.2E-02
NONHSAT201204.1:5354-5936	Zinc finger protein 717	–4.11	1.3E-02
NONHSAT166432.1:292-853	Crystallin λ1	–4.11	4.4E-02
NONHSAT151711.1:23-559	TM2 domain containing 1	–4.43	1.0E-02
NONHSAT213842.1:424-990	Inner mitochondrial membrane peptidase subunit 2	–4.47	4.6E-03
NONHSAT138481.2:3-228	Dihydrofolate reductase	–4.61	6.8E-03
NONHSAT013526.2:2-248	Catalase	–5.20	1.3E-02
NONHSAT169493.1:76-676	Aldehyde dehydrogenase 6 family member A1	–5.23	3.5E-02
NONHSAT204499.1:5-497	Chromodomain helicase DNA binding protein 1, non-coding RNA	–5.23	1.9E-02
NONHSAT152229.1:54-776	Thioredoxin interacting protein	–5.54	1.7E-02
NONHSAT195514.1:225-946	Solute carrier family 4 member 7	–5.81	1.6E-03
NONHSAT162069.1:27-809	Leucine rich repeat kinase 2	–5.90	4.0E-02
NONHSAT221848.1:287-1259	Protocadherin related 15	–7.03	3.8E-02

Table 4.

Differentially Expressed lncRNAs That Contained Sequences Similar to Genes or Proteins Related to POAG, Steroid Metabolism, Cell Adhesion, or Inflammation

Gene Name	Gene Homology	Fold Change	P
POAG-related genes and proteins
NONHSAT106523.2:1153-1748	Forkhead box C1	2256.23	4.54E-02
NONHSAT148308.1:1169-1667	Optineurin	1147.99	4.80E-02
NONHSAT193696.1:1-739	TGF-β receptor 2	–2.09	3.38E-02
NONHSAT089825.2:124-841	TGF-β receptor associated protein 1	–8.22	1.55E-02
Steroid metabolism
NONHSAT148979.1:0-1359	24-Dehydrocholesterol reductase	2.03	9.74E-03
NONHSAT169506.1:0-838	Ergosterol biosynthesis 28 homolog	2.78	1.56E-03
NONHSAT213073.1:0-1897	G protein-coupled estrogen receptor 1	5.23	2.64E-02
NONHSAT195547.1:0-1485	Oxysterol binding protein like 10	3.25	3.35E-02
Cell adhesion
NONHSAT024227.2:247-402	Neural cell adhesion molecule 1	204,369.47	3.02E-02
NONHSAT200405.1:6-227	Protocadherin 10	449,931.40	2.30E-02
NONHSAT200404.1:2-362	Protocadherin 10	39076.04	3.81E-02
NONHSAT128959.2:135-645	Protocadherin 10	6607.54	1.21E-02
NONHSAT060854.2:68-438	Protocadherin 10	5069.50	3.21E-02
NONHSAT221848.1:287-1259	Protocadherin related 15	–7.03	3.76E-02
Inflammation
NONHSAT021960.2:14-890	Vascular endothelial growth factor B	4054.53	2.00E-02
NONHSAT182278.1:5-557	Interleukin-1 receptor type 1	–2.19	4.15E-02
NONHSAT166437.1:6-363	Interleukin-17D	–13.84	2.88E-02
NONHSAT123398.2:8765-9729	Nuclear factor-κB activating protein	19,639.21	4.49E-02

miRNA Differential Expression Analysis

Using the NanoString nCounter Human miRNA Assay for 800 preselected human miRNAs, we detected the differential expression of 350 miRNAs in HTM cell culture (Supplementary Table S6). Using a cutoff of P < 0.05, Δ counts > 5, and absolute FC > 1.3, 42 miRNAs were differentially expressed in cyclic mechanically stretched TM cells (Table 5). The most significantly differentially expressed miRNAs were miR-4286, miR-29a-3p, miR-100-5p, miR-21-5p, and miR-151a-3p (P = 3.6E-5, 1.6E-3, 5.3E-3, 6.9E-3, and 7.3E-3, respectively). The miRNAs with the greatest fold change in response to stretch included miR-4286, miR-29a-3p, miR-100-5p, miR-32-5p, and miR-151a-3p, with fold changes ranging from 1.84 to 2.63. Consistent with previous reports, miR-100-5p (an absolute FC = 1.96, P = 3.6E-5), miR-27b-3p (an absolute FC = 1.76, P = 2.3E-2), miR-24-3p (an absolute FC = 1.62, P = 3.0E-2), miR-27a-3p (an absolute FC = 1.60, P = 2.7E-2), and miR-22-3p (an absolute FC = 1.54, P = 4.0E-2) were significantly upregulated in stretched compared with control TM cells. IPA showed that these miRNAs were associated with cellular functions related to cell cycle, inflammation, fibrosis, cytoskeleton, cell adhesion, endocytosis, cell contraction, migration, Wnt/β-catenin signaling, and sterol, hormone, and lipid metabolism signaling. We found that 22 of these 42 miRNAs were expressed in non-glaucomatous HTM tissues with normalized sequencing counts greater than 10 (Fig. 2).

Table 5.

Differentially Expressed miRNAs (an absolute FC > 1.3, P < 0.05) in Primary Human TM Cells in Response to Cyclic Mechanical Stretch (15%, 1 cycle/s, 24 h)

miRNA ID	Fold Change	P	miRNA ID	Fold Change	P
hsa-miR-4286	2.63	3.6E-05	hsa-miR-29b-3p	1.63	2.3E-02
hsa-miR-29a-3p	2.18	1.6E-03	hsa-miR-140-5p	1.63	3.4E-02
hsa-miR-100-5p	1.96	5.3E-03	hsa-miR-136-5p	1.62	2.4E-02
hsa-miR-32-5p	1.93	1.1E-02	hsa-miR-24-3p	1.62	3.0E-02
hsa-miR-151a-3p	1.84	7.3E-03	hsa-miR-27a-3p	1.60	2.7E-02
hsa-miR-4284	1.83	8.4E-03	hsa-miR-222-3p	1.59	2.1E-02
hsa-miR-93-5p	1.82	2.3E-02	hsa-miR-3615	1.58	2.6E-02
hsa-miR-21-5p	1.82	6.9E-03	hsa-miR-376c-3p	1.58	3.4E-02
hsa-miR-25-3p	1.82	1.4E-02	hsa-miR-642a-3p	1.58	1.8E-02
hsa-miR-27b-3p	1.76	2.3E-02	hsa-miR-34a-5p	1.57	1.8E-02
hsa-miR-127-3p	1.74	1.4E-02	hsa-miR-31-5p	1.56	2.6E-02
hsa-miR-377-3p	1.73	2.7E-02	hsa-miR-22-3p	1.54	4.0E-02
hsa-miR-15a-5p	1.72	3.7E-02	hsa-miR-125b-5p	1.53	4.7E-02
hsa-miR-181a-5p	1.67	9.5E-03	hsa-miR-99b-5p	1.51	4.0E-02
hsa-miR-30a-5p	1.67	1.5E-02	hsa-miR-185-5p	1.49	4.2E-02
hsa-miR-376a-3p	1.66	1.6E-02	hsa-miR-378i	1.49	3.4E-02
hsa-miR-574-3p	1.66	9.5E-03	hsa-miR-337-5p	1.49	4.3E-02
hsa-miR-125a-5p	1.65	9.1E-03	hsa-miR-3690	1.49	4.8E-02
hsa-miR-191-5p	1.65	3.9E-02	hsa-miR-1275	–1.55	2.4E-02
hsa-miR-190a-5p	1.65	1.6E-02	hsa-miR-187-3p	–1.64	1.7E-02
hsa-miR-181c-5p	1.64	1.7E-02	hsa-miR-1323	–1.66	1.1E-02

Figure 2.

Expression level of cyclic mechanic stretch-responsive miRNAs (an absolute FC > 1.3, P < 0.05) in non-glaucomatous human TM tissue (n = 5). Error bars represent SEM.

Using the miRNA Target Finder in IPA, we identified 18 miRNAs that target 24 mRNA genes in response to cyclic mechanical stretch in HTM cell culture (Table 6). Most of these miRNAs targeted more than one gene, and a couple of genes were targeted by more than one miRNA (Fig. 3). For example, increased expression of miR-30a-5p inhibited the expression of SLC7A11, ATP8A1, and PFN2, whereas reduced expression of PFN2 was associated with increased expression of miR-30a-5p, miR-93-5p, and miR-151a-3p. Of the differentially expressed miRNAs, three appeared to be master switches affecting the expression of three or more genes: miR-125a-5p, miR-30a-5p, and miR-1275 (Fig. 4). This integrative miRNA target analysis identified potential regulatory expression networks in response to cyclic mechanical stretch of primary HTM cells.

Table 6.

Differentially Expressed miRNAs (an absolute FC > 1.3, P < 0.05) and Their mRNA Target Genes

miRNA ID	Fold Change	Gene Symbol	Fold Change	Function
hsa-miR-125a-5p	1.65	TTC30B	–2.79	Transporter
hsa-miR-125a-5p	1.65	ID2	–2.66	Transcription regulator
hsa-miR-125a-5p	1.65	HK2	–2.47	Kinase
hsa-miR-125a-5p	1.65	LYRM9	–2.47	Unknown
hsa-miR-30a-5p	1.67	SLC7A11	–2.47	Transporter
hsa-miR-30a-5p	1.67	ATP8A1	–2.42	Lipid transporter
hsa-miR-30a-5p	1.67	PFN2	–2.11	Cytoskeletal component
hsa-miR-1275	–1.55	CADM4	291.13	Cell adhesion
hsa-miR-1275	–1.55	PURG	3.04	Transcription regulator
hsa-miR-1275	–1.55	PPT2	2.36	Lysosomal enzyme
hsa-miR-29a-3p	2.18	PCSK5	–2.75	Endopeptidase
hsa-miR-29a-3p	2.18	TPK1	–2.32	Ion channel
hsa-miR-32-5p	1.93	DUSP5	–3.09	Phosphatase
hsa-miR-32-5p	1.93	ASPH	–2.21	Endoplasmic reticulum enzyme
hsa-miR-93-5p	1.83	CXCL8	–7.44	Chemokine
hsa-miR-93-5p	1.83	PFN2	–2.11	Cytoskeletal component
hsa-miR-27b-3p	1.77	RPS6KA5	–2.79	Kinase
hsa-miR-27b-3p	1.77	KITLG	–2.54	Kinase activator
hsa-miR-15a-5p	1.72	SGCA	–2.85	Cytoskeletal component
hsa-miR-15a-5p	1.72	KITLG	–2.54	Kinase activator
hsa-miR-34a-5p	1.57	TNS2	–2.18	Phosphatase
hsa-miR-34a-5p	1.57	MYC	–2.08	Transcription regulator
hsa-miR-4286	2.63	GRASP	–2.41	Endosome regulator
hsa-miR-151a-3p	1.84	PFN2	−2.11	Cytoskeletal component
hsa-miR-181a-5p	1.67	CEP83	–2.03	Cytoskeletal component
hsa-miR-136-5p	1.62	TPK1	–2.32	Ion channel
hsa-miR-24-3p	1.62	MYC	–2.08	Transcription regulator
hsa-miR-3615	1.58	MAFB	–2.13	Transcription regulator
hsa-miR-22-3p	1.54	ATP8A1	–2.42	ATPase
hsa-miR-337-5p	1.49	HIST1H4E	–3.17	Transcription regulator
hsa-miR-185-5p	1.49	LYRM9	–2.47	Unknown

Figure 3.

Gene network of stretch-responsive miRNAs (an absolute FC > 1.3, P < 0.05) and their validated target genes (an absolute FC > 2, P < 0.05). Upregulated RNAs appear in green and downregulated RNAs appear in orange.

Figure 4.

Three miRNAs with a negative correlation with three or more differentially expressed genes were identified and considered to be master regulators. Data are mean fold change ± SEM.

ddPCR Validation of Differentially Expressed mRNAs and miRNAs

In order to validate the data obtained via RNA-Seq and miRNA arrays, we analyzed the differential expression of select protein-coding genes (i.e., ACAT2, ACSS2, DHCR7, EBP, NSDHL, GPER1, and PFN2) and miRNAs (i.e., miR-27a-3p, miR-29a-3p, miR-181c-5p, and miR-4286). The protein-coding genes were selected based on their appearance in the mRNA–miRNA interaction analysis and/or their role in steroid metabolism. The miRNAs analyzed were selected due to their appearance in the mRNA–miRNA interaction analysis, their large fold change, and/or their differential expression in previous studies. Each of these had reached the cutoff for significant differential expression in the RNA-Seq analysis (an absolute FC > 2, P > 0.05) or the miRNA array (an absolute FC > 1.3, P > 0.05). Of the seven protein-coding genes examined, five (i.e., ACAT2, ACSS2, DHCR7, EBP, and NSDHL) were validated by ddPCR as being differentially expressed (an absolute FC > 2, P > 0.05) (Fig. 5A; Table 7). Meanwhile, of the four miRNAs examined, three (i.e., miR-27a-3p, miR-29a-3p, and miR-181c-5p) were validated by ddPCR as being differentially expressed (an absolute FC > 1.3, P > 0.05) (Fig. 5B; Table 8).

Figure 5.

(A) The expression of seven stretch-responsive protein-coding genes was examined with ddPCR. These genes were related to steroid metabolism and/or were central hubs in the mRNA–miRNA interaction analysis. The differential expression of five of the seven mRNAs was validated by ddPCR. (B) The expression of four stretch-responsive miRNAs was analyzed with ddPCR. These miRNAs have been identified previously, had large fold changes, and/or were identified in the mRNA–miRNA interaction analysis. The differential expression of three of the four miRNAs was validated by ddPCR. Data are mean fold change ± SEM, where *P < 0.05, **P < 0.01, and ***P < 0.001 were found by paired two-tailed Student's t-test (n = 5).

Table 7.

Expression of Seven Stretch-Responsive Protein-Coding Genes Examined with ddPCR

		RNA-Seq		ddPCR
Gene Symbol	Gene Name	Fold Change	P	Fold Change	P
ACAT2	Acetyl-coenzyme A acetyltransferase	2.43	9.42E-03	2.51	3.41E-03
ACSS2	Acetyl-coenzyme A synthetase	2.07	4.27E-02	2.78	1.97E-02
DHCR7	7-Dehydrocholesterol reductase	2.21	4.45E-03	2.65	1.46E-02
EBP	3-β-Hydroxysteroid-δ(8),δ(7)-isomerase	2.20	2.48E-03	2.97	3.74E-03
NSDHL	Sterol-4-α-carboxylate 3-dehydrogenase	2.24	3.62E-03	2.72	4.09E-03
GPER1	G-protein coupled estrogen receptor 1	3.17	2.73E-02	0.96	9.89E-01
PFN2	Profilin-2	–2.11	1.69E-04	1.07	9.52E-01

These genes were related to steroid metabolism and/or were central hubs in the mRNA–miRNA interaction analysis. The differential expression of five of the seven mRNAs was validated by ddPCR.

Table 8.

Expression of Four Stretch-Responsive miRNAs Analyzed with ddPCR

	NanoString Assay		ddPCR
miRNA ID	Fold Change	P	Fold Change	P
miR-27a-3p	1.60	2.69E-02	2.80	1.47E-02
miR-29a-3p	2.18	1.57E-03	2.38	2.72E-02
miR-181c-5p	1.64	1.75E-02	2.47	4.87E-02
miR-4286	2.63	3.61E-05	1.74	2.78E-01

These miRNAs have been identified previously, had large fold changes, and/or were identified in the mRNA–miRNA interaction analysis. The differential expression of three of the four miRNAs was validated by ddPCR.

Discussion

Overview

We identified 219 protein-coding mRNAs, 387 lncRNAs, and 42 miRNAs differentially expressed in primary HTM cells in response to cyclic mechanical stretch. Differentially expressed mRNAs were involved in cell cycle regulation and sterol and lipid metabolism and included genes related to POAG. Meanwhile, differentially expressed lncRNA sequences showed sequence similarity to POAG-related genes as well as genes involved in cell adhesion, inflammation, and steroid metabolism. Our miRNA analysis identified alterations in miRNAs potentially related to cell cycle, inflammation, fibrosis, cytoskeleton, cell adhesion, endocytosis, cell contraction, migration, Wnt/β-catenin signaling, and sterol, hormone, and lipid metabolism signaling. Pathway analysis of coding RNAs indicated the potential involvement of steroid biosynthesis, glycerolipid metabolism, and ECM–receptor interaction pathways. Although our findings validate those of other mechanotransduction studies,,,,, to our knowledge this is the first time that the steroid biosynthesis pathway has been associated with cyclic mechanical stretch of HTM cells. This finding is significant given the history of epidemiological evidence implicating steroids (i.e., estrogen, androgens, and corticosteroids) in risk for POAG development.

Differentially Expressed POAG-Associated Genes

A couple of differentially expressed mRNAs (i.e., UCP2 and TXNIP) have been shown previously to play a role in the pathophysiology of POAG.– These differentially expressed mRNAs have functions related to oxidative stress response. Overexpression and knockout of UCP2 have been shown to affect mitochondrial function, mitophagy, and RGC survival.– The thioredoxin-interacting protein gene (TXNIP) may also play a role in oxidative stress response. Thioredoxin is known to be involved in sustaining a healthy redox state in order to minimize the effects of oxidative stress. In addition, TXNIP has been shown to play a role in inflammatory response., Knockout of TXNIP has been shown to attenuate cytokine expression and release, inflammasome activation, and glial cell activation in a neurotoxicity model of retinal degeneration, thus suggesting that TXNIP is a critical component of inflammatory response.

lncRNA Sequence Similarity with POAG-Associated Genes

lncRNAs are non-coding RNAs that are greater than 200 bases in length. Although they share similar properties with mRNAs, they have regulatory functions more similar to those of miRNAs. They may enact their regulatory roles in gene transcription through epigenetic mechanisms or through interactions with transcription factors or other gene activating/suppressing complexes.– lncRNAs may also act post-transcriptionally or post-translationally to modulate gene expression through alternative splicing, miRNA interactions, transport, or protein modification.,–

A few differentially expressed lncRNAs shared sequence similarity to genes associated with POAG, including FOXC1 and OPTN. The gene for Forkhead Box C1 (FOXC1) is positioned on chromosome 6 next to the gene for GDP-mannose 4,6-dehyrdratase (GMDS). Variants in and around both genes have been associated with POAG, as well as developmental and primary congenital glaucoma.– The mechanism of the contribution of these variants to glaucoma pathogenesis is unknown. However, the expression of a lncRNA RP11-157J24.2 (ELF2P2, E74-like factor 2 pseudogene 2) has been associated with the FOXC1 single nucleotide polymorphism rs2745572, suggesting that POAG-associated variants near FOXC1 may contribute to POAG pathogenesis through lncRNA regulation.

Similarly, mutations in the chromosome 10 optineurin gene (OPTN) have been identified as contributing to the development of POAG by linkage analysis.,, OPTN has been linked specifically to normal tension glaucoma, a POAG subtype in which patients demonstrate the retinal characteristics of glaucoma while having phenotypically normal levels of IOP.,, The normal function of OPTN includes inhibition of TNFα-induced inflammation, cell division, vesicular transport, pathogen defense, and autophagy, including autophagy of defective mitochondria.,– Although the exact mechanism for the contribution of OPTN to POAG development is under continued investigation, its role as an autophagy receptor has been of particular interest.,

Other lncRNAs shared sequence similarity with genes related to cell adhesion, steroid metabolism, inflammation, and TGF-β, a well-known role-player in POAG pathogenesis.–

TGF-β Signaling Pathway and ECM Response

Cyclic mechanical stretch has been shown to potentially activate the TGF-β/bone morphogenetic protein pathway through the SMAD-mediated Runx2 pathway or the non-SMAD-mediated mitogen-activated protein kinase pathway., TGF-β and its related pathways have been strongly implicated in the pathogenesis of glaucoma and IOP regulation.– TGF-β may increase outflow resistance by altering ECM homeostasis and cell contractility in the TM through interactions with other proteins and signaling molecules.,,

Our pathway analysis of differentially expressed protein-coding mRNAs implicated changes in the ECM–receptor interaction pathway. A couple of the differentially expressed mRNAs were located in the ECM or extracellular regions or had functions related to ECM homeostasis. Furthermore, one of the differentially expressed lncRNA sequences showed similarity to TGF-β receptor 2, and two other sequences showed similarity to TGF-β receptor associated protein 1. In addition, we observed an upregulated response of lncRNAs with sequence similarity to genes related to the ECM and cell adhesion.

TGF-β treatment results in the regulation of many miRNAs, including the upregulation of miR-21, miR-181, miR-494, miR-10b, miR-27a, miR-183, miR-182, miR-155, and miR-451 and the downregulation of miR-200, miR-34a, miR-203, miR-584, and miR-450b-5p. Most members of the TGF-β pathway may be targeted by a number of miRNAs, including miR-18a, miR-24, let-7, miR-744, miR-30, miR-200, miR-128a, miR-21, miR-17, miR-148a, miR-99a/b, and miR-92b., Our study identified the differential expression of miR-21-5p, miR-181a-5p, miR-181c-5p, miR-27a-3p, miR-34a-5p, miR-24-3p, miR-30a-5p, miR-21-5p, and miR-99b-5p, all of which potentially target genes involved in TGF-β pathway-mediated ECM homeostasis in TM cells.

Steroid Metabolism and Estrogen Signaling Pathway

Corticosteroid use has been known to induce elevated IOP and glaucoma.– In fact, treatment with the steroid dexamethasone has been used to create cell and animal models of POAG.– Dexamethasone treatment results in increased ECM deposition in the TM, thereby reducing outflow., More recently, there is evidence suggesting that not only corticosteroids but also other steroids, including sex hormones, may play a role in glaucoma pathophysiology. Sterol carrier protein 2 was identified by Vittal et al. in their study examining changes in TM gene expression following mechanical stretch, and androgen receptor has been implicated in glaucoma pathogenesis. Furthermore, several epidemiological studies have demonstrated differences in POAG risk between males and females as well as among females at different stages of reproductive capacity, suggesting that estrogen may play a role in POAG pathophysiology.–

Our pathway analysis of the differentially expressed protein-coding genes suggests that steroid metabolism pathways are responsive to cyclic mechanical stretch. Moreover, the mRNA for G protein-coupled estrogen receptor 1 was differentially expressed in response to cyclic stretch in the RNA-Seq data (an absolute FC = 3.17, P = 2.7E-2) (Fig. 5A; Table 7); however, this gene was not validated by ddPCR analysis (Fig. 5A; Table 7). Nevertheless, several other differentially expressed genes (i.e., ABCA5, ACAT2, ACSS2, CCNA2, DHCR7, EBP, FHL2, HN1, LIPG, NSDHL, PCSK9, PRKAA2, SC5D, TAF11, TMEM97, and TXNIP) share lipid, sterol, or hormone metabolism as a biological process gene ontology term (https://www.uniprot.org/). The differential expression of five of these (i.e., ACAT2, ACSS2, DHCR7, EBP, and NSDHL) was validated by ddPCR (Fig. 5A; Table 7). In addition, several differentially expressed lncRNAs shared sequence homology to steroid hormone metabolism genes, including G protein-coupled estrogen receptor 1, oxysterol binding protein like 10, ergosterol biosynthesis 28 homolog, and 24-dehydrocholesterol reductase. These results suggest that steroid metabolism, and estrogen signaling in particular, are responsive to mechanical stretch of the TM and may therefore be involved in IOP homeostasis.

miRNA Regulation of Stretch Responsive Genes

Several of our identified miRNAs (i.e., miR-100, miR-27a, mir-27b, miR-22, and miR-24) have been reported to be differentially expressed in TM cells after 3 hours of cyclic mechanical stretching. Meanwhile, miR-15a has been reported to be downregulated in stress-induced senescent HTM cells. We also identified miR-93, which has been found to increase expression and induce apoptosis in glaucomatous TM cells. Our integrative miRNA–mRNA expression analysis identified several miRNA–mRNA interactions in the TM cells, highlighting several miRNA master regulators, such as miR-125a-5p, miR-30a-3p, and miR-1275. Additional master miRNA regulators may be identified if we expand our analysis to include predicted targets with moderate confidence, although this would also potentially increase the number of false-positive findings.

Conclusions and Future Work

Our study validated the results of previous studies by using an integrative approach to examine the coding and non-coding RNAs expressed in the same set of cells subjected to cyclic mechanical stretch. In addition to replicating previously reported differentially expressed RNAs, we also identified new RNAs responsive to cyclic mechanical stretch in TM cells, including mRNAs, lncRNAs, and miRNAs. Despite the strength of our experimental and bioinformatics analyses, our study has a few limitations. First, having more donor-derived TM cells with similar ethnic backgrounds and age distributions could provide our study with greater statistical power. Additionally, having more eyes from both males and females could help elucidate differences in steroid metabolism and estrogen signaling gene response to mechanical stretch. Second, the cyclic mechanical stretch experiment was done at a single time point for 24 hours; therefore, including a time-series design would be beneficial in identifying time-dependent molecular responses. Third, all of the TM cells used in our study were derived from unaffected non-glaucomatous postmortem donors. Including TM cells derived from glaucoma-affected postmortem donors will be necessary to further explore disease-specific mechanical stretch-induced changes in expression. These factors will be considered and included in our future experiments.

In summary, we conducted a genome-wide analysis of mRNA, lncRNA, and miRNA expression in the same set of cells in response to mechanical stretch. Our RNA expression profiling has provided valuable foundational data, identifying a large number of differentially expressed genes and miRNAs in response to cyclic mechanical stretch in primary HTM cells, validating previous reports, and identifying novel targets. Our analysis has identified several important signaling pathways involved in this response, such as ECM–receptor interaction and steroid metabolism. The miRNA–mRNA integrative analysis found several miRNA master regulators, suggesting their potential role in TM cellular function in response to cyclic mechanical stretch. Additional functional studies will help further validate the role of these RNAs in relation to TM cellular function.