Overexpression and Activation of αvβ3 Integrin Differentially Affects TGFβ2 Signaling in Human Trabecular Meshwork Cells

Studies from our laboratory have suggested that activation of αvβ3 integrin-mediated signaling could contribute to the fibrotic-like changes observed in primary open angle glaucoma (POAG) and glucocorticoid-induced glaucoma. To determine how αvβ3 integrin signaling could be involved in this process, RNA-Seq analysis was used to analyze the transcriptomes of immortalized trabecular meshwork (TM) cell lines overexpressing either a control vector or a wild type (WT) or a constitutively active (CA) αvβ3 integrin. Compared to control cells, hierarchical clustering, PANTHER pathway and protein-protein interaction (PPI) analysis of cells overexpressing WT-αvβ3 integrin or CA-αvβ3 integrin resulted in a significant differential expression of genes encoding for transcription factors, adhesion and cytoskeleton proteins, extracellular matrix (ECM) proteins, cytokines and GTPases. Cells overexpressing a CA-αvβ3 integrin also demonstrated an enrichment for genes encoding proteins found in TGFβ2, Wnt and cadherin signaling pathways all of which have been implicated in POAG pathogenesis. These changes were not observed in cells overexpressing WT-αvβ3 integrin. Our results suggest that activation of αvβ3 integrin signaling in TM cells could have significant impacts on TM function and POAG pathogenesis.


Introduction
Glaucoma is a chronic optic neuropathy that results in irreversible blindness [1]. Although the exact etiology of glaucoma remains unknown, a number of studies have now determined that elevated levels of transforming growth factor β2 (TGFβ2) and connective tissue growth factor (CTGF) and the use of glucocorticoids can lead to the development of glaucoma [2,3]. Clinically, glaucoma induced by TGFβ2, CTGF and glucocorticoids (GC) is very similar. All cause an increase in extracellular matrix (ECM) deposition and cytoskeleton changes associated with contractility that lead to a restriction in aqueous humor outflow and together contribute to the elevation of intraocular pressure (IOP) associated with glaucoma.
Recent studies have shown that the activation of αvβ3 integrin signaling plays a role in regulating these changes in ECM deposition and contractility. Activation of the integrin enhances the deposition of EDA+ fibronectin into the ECM [4] which in turn can trigger profibrotic changes in the ECM and elevation in IOP [5,6]. The αvβ3 integrin also promotes the formation of crosslinked actin networks (CLANs) that may cause the contractility changes in the TM thought to be associated with a stiffer fibrotic-like tissue [7][8][9][10]. Activation of αvβ3 integrin also increases IOP in mice and porcine organ cultured anterior segments. In contrast, knockdown of the αvβ3 integrin gene and its expression in the TM and ciliary muscle (CM) using a tamoxifen inducible Cre +/-β3 integrin flox/flox mouse model decreases IOP in mice [11]. The activity of the αvβ3 integrin in human trabecular meshwork (HTM) cells can be increased by the glucocorticoid dexamethasone via the calcineurin/NFATc1 pathway [12]. This pathway has been associated with glaucoma since it regulates the expression of myocilin, a protein involved in an early onset form of POAG and GC-induced ocular hypertension [13]. Together these results suggest that αvβ3 integrin signaling can contribute to the development of glaucoma.
αvβ3 integrin is a member of the integrin family of transmembrane receptors that mediate adhesion to the ECM. In addition to their mechanical roles in mediating adhesion, integrins transmit chemical signals into the cell that provided information on a cell's local environment [14]. These signals determine a cell's ability to respond to other signals transmitted by growth factor-or G-protein-coupled receptors. Some of the major pathways co-regulated by integrins include PI3K/Akt, MAPK/Erk, Jnk, and Rho GTPase [15]. In addition, integrins can mediate gene expression independent of their adhesion role [16,17]. Thus, activation of integrins has the potential to play several critical roles in regulating signaling pathways proposed to be altered in glaucoma.
To understand how αvβ3 integrin could affect the signaling pathways associated with glaucoma in the trabecular meshwork, RNA-Seq analysis was used to study changes in the transcriptome of TM cells when αvβ3 integrin was activated. The study used the immortalized TM-1 cell line in which either a wild type (WT) or a constitutively active (CA) form of the αvβ3 integrin was overexpressed [18]. TM-1 cells transduced with an empty vector were used as a control. The differential expression of genes between these cell lines was analyzed using various bioannotations, network and pathway analysis databases to determine specific signaling pathways and categories of genes that may be altered when αvβ3 integrin is active.

Cell Culture
Human SV40-immortalized TM-1 cells where generated from primary human TM (HTM) cells [18] that were derived from a 30-year-old donor eye as previously described [19,20]. To the best of our knowledge, the donor eye issue was free of any ocular diseases. The isolated parent HTM cells showed morphological features expected of human trabecular meshwork cells [19,20] and expressed elevated levels of myocilin in response to glucocorticoids [21]. The TM-1 cells generated from these primary HTM cells appeared similar to HTM cells and meet the criteria outlined in Keller, et al. [22]. They upregulated expression of both fibronectin and myocilin mRNA in response to glucocorticoids [23] and were capable of forming CLANs [8]. TM-1 cells overexpressing a control empty vector (EV), a wild type αvβ3 integrin (WTβ3) or a CA-αvβ3 integrin (CAβ3) that contains a T562N mutation in the β3 integrin subunit [24] were subsequently generated as described by Gagen, et al. [25]. All the cell lines were cultured in low glucose DMEM (MilliporeSigma, St. Louis, MO, USA), containing 10% fetal bovine serum (FBS, Atlanta Biologicals/R&D Systems, Flowery Branch, GA, USA), 2 mM L-glutamine (MilliporeSigma), 0.2% Primocin (Thermo Fisher Scientific, Waltham, MA) and 0.05% gentamicin (Corning, Glendale, AZ, USA). Cells were kept under selective pressure using 2 µg/mL puromycin.

RT-qPCR Analysis
cDNA was generated with the High-capacity cDNA reverse transcription kit (Thermo Fisher Scientific) according to the manufacturer's instructions. qPCR was performed using an Applied Biosystems QuantStudio 6 Pro Real-Time PCR system and PowerUp SYBR green master mix (Thermo Fisher Scientific) as we described [23]. Fold changes in gene expression were determined using the ∆∆Ct method. Data were normalized to the housekeeping gene succinate dehydrogenase complex subunit A (SDHA). Primer-BLAST (National Center for Biotechnology Information (NCBI), Bethesda, MD, USA) was used to design the primers which were made by IDT (Integrated DNA Technologies, Inc., Coralville, IA, USA).

R Analysis
R statistical computing and graphics software was used to identify unique or shared differentially expressed genes in WTβ3 and CAβ3 cells that demonstrated a ±1.5 Log 2 FC. VennDiagram, readxl, tidyverse and dplyr R packages were used for this purpose. Generation of heat maps was performed using readxl, dplyr and gplots R packages. In some instances, the FC of a gene in one of the three contrasts was missing from the dataset because that gene demonstrated a subthreshold read count in one of the three cell lines in the RNA-Seq analysis (see above). In order to display all three contrasts for these genes within the heat maps, the missing FCs were calculated in Microsoft Excel using the pre-filtered FPKM (Fragments Per Kilobase of transcript per Million) counts from the raw RNA-Seq data using the following formula: Cells 2021, 10, 1923 4 of 33 In these instances, the Log 2 FC was manually calculated for all three contrasts so the calculations would be consistent within the heat map. These genes are identified by asterisks in the heat maps.

TGFβ2 ELISA Analysis
TM-1 EV or CAβ3 cells were plated at a density of 2.0 × 10 5 cells/4 cm 2 in growth medium containing 10% FBS and allowed to attach overnight. Cells were incubated for 24 hr in medium containing 1% FBS. The conditioned media (CM) was collected and debris removed via centrifugation. The CM was acidified by addition of 1N HCL for 10 min in order to activate the TGFβ2. The CM was then neutralized by addition of NaOH in 0.5 M HEPES. ELISA analysis was performed using an R&D Systems Human TGFbeta 2 Quantikine ELISA Kit (R & D Systems, Minneapolis, MN), and the procedure was performed according to the manufacturer's instructions. Samples were read on a BioTek Synergy Neo2 plate reader (Agilent, Santa Clara, CA, USA) at 450 nm. Statistical analysis of the data was performed using analysis of variance (ANOVA).

Bioannotations, Network and Pathway Analysis
Networks were generated using the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database v. 11.0b [30] to identify known and predicted protein-protein interactions. Medium (0.4) and high (0.7) confidence networks were both examined. The PANTHER (Protein Analysis Through Evolutionary Relationships) database v. 16.0 [31] was also used for pathway analysis. A statistical overrepresentation test (Fisher's Exact test) was used to determine if any metabolic or signaling pathways were enriched in the WTβ3 or CAβ3 cell lines.

αvβ3 Integrin Alters the Expression Pattern of the TM Transcriptome
Since activation of αvβ3 integrin signaling can trigger changes linked to glaucoma [41], RNA-Seq analysis was performed to determine how activation of αvβ3 integrin signaling affects the transcriptome of TM cells. Figure 1 shows unsupervised multidimensional scaling (MDS) of the genes (1A) and transcripts (1B), respectively, in cells that overexpressed either a CA-αvβ3 integrin (CAβ3 cells), a WT-αvβ3 integrin (WTβ3 cells) or a control empty vector (EV cells) that lacked αvβ3 integrin. Both MDS plots show that, relative to the EV cell line, both the WTβ3 and CAβ3 cell lines represented distinct populations of cells with specific gene transcripts being upregulated or downregulated by the expression and/or activity of αvβ3 integrin. Similarly, the WTβ3 and CAβ3 cell lines were found to be distinct relative to each other. Furthermore, these plots showed that there was little variance within the four samples from each cell line. Together these data show that the differences between the cell lines were larger than the differences within each cell line. Unsupervised clustering was used to verify that the cell lines used were distinct from each other. Multidimensional scaling (MDS) plots at the gene (A) and the transcript (B) level were generated. Both plots show that the greatest variation observed was between the three cell lines rather than within the individual samples of each cell line.
Comparing the transcriptomes of the three cell lines found significant DGE between them ( Table 1). The transcriptome of WTβ3 cells showed a statistically significant (p < 0.05) difference in 5,186 genes compared to that of control EV cells. Of those altered genes, 2,645 were upregulated in WTβ3 cells while 2,541 were downregulated. We also found a statistically significant difference (p < 0.05) in 8,249 genes in the transcriptome of CAβ3 cells compared to EV cells. In this case, 4,094 genes were upregulated in CAβ3 cells while 4,255 were downregulated. The transcriptome of WTβ3 cells compared to CAβ3 cells also differed. The WTβ3 transcriptome showed a significant difference (p < 0.05) in 6,999 genes compared to that of CAβ3 cells (3,524 upregulated; 3,475 downregulated). Table 1. Summary of differential gene expression analysis for TM-1 cells overexpressing either a wild type β3 integrin (WTβ3) or a constitutively active β3 integrin (CAβ3). Control cells were transduced with an empty vector (EV). These data can be accessed in the supplemental file: edgeRglm_GENE_3 contrasts_low-high_logFC_no NS genes.xlsx. The complete data set has also been deposited in NCBI's Gene Expression Omnibus and is accessible through GEO Series. Accession number: GSE180407. # = number.

Contrast
Total

RNA-Seq Analysis Reveals Genes Enriched or Depleted by Activation of αvβ3 Integrin Signaling
We then identified genes that demonstrated a ± 1.5 Log2 FC in expression in the WTβ3 and CAβ3 cell lines compared to the control EV cells (Tables S1-S6). As shown in Figure  2, overexpression of the WT-αvβ3 integrin caused an enrichment in 113 genes and a depletion in 73 genes. In contrast, expression of the CA-αvβ3 integrin caused an enrichment in 273 genes and a depletion in 149 genes. When comparing the gene expression changes, we saw that 65 genes were enriched in both cell lines while 34 genes were depleted in both cell lines. This suggests that activation of αvβ3 integrin signaling differentially affected the TM-1 transcriptome compared to mere overexpression of WT-αvβ3 integrin. Unsupervised clustering was used to verify that the cell lines used were distinct from each other. Multidimensional scaling (MDS) plots at the gene (A) and the transcript (B) level were generated. Both plots show that the greatest variation observed was between the three cell lines rather than within the individual samples of each cell line.
Comparing the transcriptomes of the three cell lines found significant DGE between them ( Table 1). The transcriptome of WTβ3 cells showed a statistically significant (p < 0.05) difference in 5186 genes compared to that of control EV cells. Of those altered genes, 2645 were upregulated in WTβ3 cells while 2541 were downregulated. We also found a statistically significant difference (p < 0.05) in 8249 genes in the transcriptome of CAβ3 cells compared to EV cells. In this case, 4094 genes were upregulated in CAβ3 cells while 4255 were downregulated. The transcriptome of WTβ3 cells compared to CAβ3 cells also differed. The WTβ3 transcriptome showed a significant difference (p < 0.05) in 6999 genes compared to that of CAβ3 cells (3524 upregulated; 3475 downregulated). Table 1. Summary of differential gene expression analysis for TM-1 cells overexpressing either a wild type β3 integrin (WTβ3) or a constitutively active β3 integrin (CAβ3). Control cells were transduced with an empty vector (EV). These data can be accessed in the supplemental file: edgeRglm_GENE_3 contrasts_low-high_logFC_no NS genes.xlsx. The complete data set has also been deposited in NCBI's Gene Expression Omnibus and is accessible through GEO Series. Accession number: GSE180407. # = number.

Contrast
Total We then identified genes that demonstrated a ± 1.5 Log 2 FC in expression in the WTβ3 and CAβ3 cell lines compared to the control EV cells (Tables S1-S6). As shown in Figure 2, overexpression of the WT-αvβ3 integrin caused an enrichment in 113 genes and a depletion in 73 genes. In contrast, expression of the CA-αvβ3 integrin caused an enrichment in 273 genes and a depletion in 149 genes. When comparing the gene expression changes, we saw that 65 genes were enriched in both cell lines while 34 genes were depleted in both cell lines. This suggests that activation of αvβ3 integrin signaling differentially affected the TM-1 transcriptome compared to mere overexpression of WT-αvβ3 integrin.

Figure 2.
CAβ3 and WTβ3 cells have unique and shared sets of genes with altered expression. Venn diagram analysis was performed looking for differentially expressed genes that demonstrated a Log2 FC ± 1.5 in CAβ3 or WTβ3 cells relative to control EV cells. The datasets from the CAβ3-EV and WTβ3-EV contrasts were analyzed using R statistical computing and graphics software. Shared and unique sets of genes that were upregulated (A) and downregulated (B) relative to EV cells were identified. The complete lists of unique and shared genes that met the FC cutoff are found in Tables S1-S6.
This differential expression was detected in both the coding and non-coding genomic regions of both cell lines. A total of 447 genes were detected from the coding regions (Table  S7) while 62 transcripts were detected from non-coding regions (Table S8). This latter group included lncRNA, lincRNA, microRNA, snRNA, snoRNA, ribosomal RNA, novel transcripts, open-reading frames, pseudogenes and antisense genes.
To identify biological processes that may be enriched by overexpression of either WT-αvβ3 integrin or CA-αvβ3 integrin, gene ontology enrichment analysis and gene classification was performed on the genes listed in Tables S1-S6. As shown in Figure 3, these analyses showed that multiple biological processes were affected. The top biological process affected was transcription, followed by ECM, cell adhesion and cytokine signaling. Interestingly, these biological processes were differentially affected in CAβ3 cells compared to WTβ3 cells. More biological processes were affected by overexpressing CA-αvβ3 integrin than WT-αvβ3 integrin. Biological processes that were enriched in CAβ3 or WTβ3 cells were identified through gene ontology enrichment analysis and gene classification that was performed on the genes that met the Log2 FC ± 1.5 cutoff from each cell line using a combination of the following databases: CAβ3 and WTβ3 cells have unique and shared sets of genes with altered expression. Venn diagram analysis was performed looking for differentially expressed genes that demonstrated a Log 2 FC ± 1.5 in CAβ3 or WTβ3 cells relative to control EV cells. The datasets from the CAβ3-EV and WTβ3-EV contrasts were analyzed using R statistical computing and graphics software. Shared and unique sets of genes that were upregulated (A) and downregulated (B) relative to EV cells were identified. The complete lists of unique and shared genes that met the FC cutoff are found in Tables S1-S6.
This differential expression was detected in both the coding and non-coding genomic regions of both cell lines. A total of 447 genes were detected from the coding regions (Table S7) while 62 transcripts were detected from non-coding regions (Table S8). This latter group included lncRNA, lincRNA, microRNA, snRNA, snoRNA, ribosomal RNA, novel transcripts, open-reading frames, pseudogenes and antisense genes.
To identify biological processes that may be enriched by overexpression of either WT-αvβ3 integrin or CA-αvβ3 integrin, gene ontology enrichment analysis and gene classification was performed on the genes listed in Tables S1-S6. As shown in Figure 3, these analyses showed that multiple biological processes were affected. The top biological process affected was transcription, followed by ECM, cell adhesion and cytokine signaling. Interestingly, these biological processes were differentially affected in CAβ3 cells compared to WTβ3 cells. CAβ3 and WTβ3 cells have unique and shared sets of genes with altered expression. Venn diagram analysis was performed looking for differentially expressed genes that demonstrated a Log2 FC ± 1.5 in CAβ3 or WTβ3 cells relative to control EV cells. The datasets from the CAβ3-EV and WTβ3-EV contrasts were analyzed using R statistical computing and graphics software. Shared and unique sets of genes that were upregulated (A) and downregulated (B) relative to EV cells were identified. The complete lists of unique and shared genes that met the FC cutoff are found in Tables S1-S6.
This differential expression was detected in both the coding and non-coding genomic regions of both cell lines. A total of 447 genes were detected from the coding regions (Table  S7) while 62 transcripts were detected from non-coding regions (Table S8). This latter group included lncRNA, lincRNA, microRNA, snRNA, snoRNA, ribosomal RNA, novel transcripts, open-reading frames, pseudogenes and antisense genes.
To identify biological processes that may be enriched by overexpression of either WT-αvβ3 integrin or CA-αvβ3 integrin, gene ontology enrichment analysis and gene classification was performed on the genes listed in Tables S1-S6. As shown in Figure 3, these analyses showed that multiple biological processes were affected. The top biological process affected was transcription, followed by ECM, cell adhesion and cytokine signaling. Interestingly, these biological processes were differentially affected in CAβ3 cells compared to WTβ3 cells. More biological processes were affected by overexpressing CA-αvβ3 integrin than WT-αvβ3 integrin. Biological processes that were enriched in CAβ3 or WTβ3 cells were identified through gene ontology enrichment analysis and gene classification that was performed on the genes that met the Log2 FC ± 1.5 cutoff from each cell line using a combination of the following databases: Figure 3. More biological processes were affected by overexpressing CA-αvβ3 integrin than WT-αvβ3 integrin. Biological processes that were enriched in CAβ3 or WTβ3 cells were identified through gene ontology enrichment analysis and gene classification that was performed on the genes that met the Log 2 FC ± 1.5 cutoff from each cell line using a combination of the following databases: KEGG, DAVID Bioinformatics Resources, the Matrisome Project and UniProt. The top twelve processes are highlighted. The remaining processes that were affected can be found in Tables S7 and S8. Since integrin signaling pathways are key regulators of gene expression, it was not surprising that activation of αvβ3 integrin signaling caused a significant change in the expression of transcription factors and their co-activators. Genes for transcription factors and co-activators that were preferentially altered in CAβ3 cells represented 13% (43/323) of the differentially expressed CAβ3 cell transcriptome when compared to the EV cell transcriptome. In contrast, only 6% (5/87) of the genes for transcription factors and coactivators were altered in WTβ3 cells ( Figure 3). The transcriptome of genes for cytokines and cytokine-related proteins was also preferentially altered in CAβ3 cells representing 6% (20/323) of the CAβ3 transcriptome. In contrast, the transcriptome of WTβ3 cells showed only a 3% (3/87) shift in genes for this group. Both transcriptomes showed similar enrichments in genes for cell adhesion receptors, especially members of the cadherin family, and represented 6% (20/323) of the transcriptome in CAβ3 cells and 7% (6/87) of the transcriptome in WTβ3 cells. Likewise, the transcriptome of ECM-associated proteins appeared to be similarly affected in both cell types and represented 7% (21/323) of the transcriptome in CAβ3 cells and 8% (7/87) of the transcriptome in WTβ3 cells. Given that integrin signaling is known to regulate the organization of the cytoskeleton, it was somewhat surprising that only 3% of the transcriptome affected by the activation of αvβ3 integrin signaling involved cytoskeleton-associated proteins and GTPases. Not surprisingly, genes in the transcriptomes of both cell lines that were commonly shared showed an enrichment in all these groups as well (see Table S7).

Expression of Transcription Factors Is Differentially Regulated in CAβ3 and WTβ3cells
Hierarchical clustering of genes encoding transcription factors and other transcriptionassociated proteins in the ±1.5 Log 2 FC cutoff ( Figure 4) shows the differential expression of these genes in cells overexpressing WT αvβ3 integrin compared to cells overexpressing CA αvβ3 integrin. A major difference appears to be in the expression of ZNF (zinc finger) transcription factors. Surprisingly, a total of 11 ZNF transcription factor genes demonstrated significantly altered expression in either CAβ3 or WTβ3 cells when contrasted to control EV cells (Tables S1-S6). In CAβ3 cells, 9/11 were uniquely upregulated and 0/11 uniquely downregulated while in WTβ3 cells there was 1/11 uniquely upregulated and 1/11 uniquely downregulated. Although the CAβ3 transcriptome showed an enrichment in genes encoding members of the ZNF family of transcription factors, the significance of this is unknown since their functions are still largely unknown.
Several genes encoding transcription factors in those clusters are of particular interest. For instance, MYOCD (CAβ3 FC = 2.94), which is upregulated in CAβ3 cells, encodes a transcriptional coactivator of serum response factor (SRF) [42] that regulates smooth muscle contraction and is involved in fibrosis [43]. SOX5 [44] (CAβ3 FC = 3.49), STAT5A [45] (CAβ3 FC = 3.32), and FOSL1 [46] (CAβ3 FC = −1.98) are also associated with fibrosis as is KLF4 (FC = −2.58) which encodes an inhibitor of TGFβ signaling [47]. Additionally, the KLF4 transcription factor is responsive to shear stress, downregulates SPARC expression in the TM and plays a role in TM stem cell development and proliferation [48][49][50]. KLF2 (CAβ3 FC = −1.54), also encodes an inhibitor of TGFβ signaling [51] that can function as a flow responsive transcription factor that induces endothelial nitric oxide synthase (eNOS) expression [52]. . Hierarchical clustering of transcription-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log2 FC ± 1.5 cutoff and were classified as either transcription factors, or otherwise part of the transcription process were analyzed using R statistical computing and graphics software. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been excluded from the initial DGE analysis due to low abundance in the read count in two or more samples of one of the three cell lines. To be included in the heat map, the missing Log2 FC values were calculated from the raw count data using Excel. To be consistent, the corresponding Log2 FC values for the same gene in the other two cell lines were re-calculated from the raw data as well. These values were then converted to standardized z-scores along with the other genes using R.
Other genes of interest are PAX6, LMX1B and DDIT3. PAX6 (FC = 2.78) is significantly upregulated in CAβ3 vs. WTβ3 cells like the genes in clusters L and M. It has been shown to play a role in anterior segment dysgenesis (ASD) [53] and is linked to glaucoma ( Table 2). LMX1B (CAβ3 FC = 1.65) encodes a transcription factor that is associated with nail patella and has been shown to be associated with the development of POAG (Table  2). DDIT3 encodes a pro-apoptotic transcription factor, also known as CHOP, that participates in regulating IOP (Table 3)   . Hierarchical clustering of transcription-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff and were classified as either transcription factors, or otherwise part of the transcription process were analyzed using R statistical computing and graphics software. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been excluded from the initial DGE analysis due to low abundance in the read count in two or more samples of one of the three cell lines. To be included in the heat map, the missing Log 2 FC values were calculated from the raw count data using Excel. To be consistent, the corresponding Log 2 FC values for the same gene in the other two cell lines were re-calculated from the raw data as well. These values were then converted to standardized z-scores along with the other genes using R.
Other genes of interest are PAX6, LMX1B and DDIT3. PAX6 (FC = 2.78) is significantly upregulated in CAβ3 vs. WTβ3 cells like the genes in clusters L and M. It has been shown to play a role in anterior segment dysgenesis (ASD) [53] and is linked to glaucoma ( Table 2). LMX1B (CAβ3 FC = 1.65) encodes a transcription factor that is associated with nail patella and has been shown to be associated with the development of POAG (Table 2). DDIT3 encodes a pro-apoptotic transcription factor, also known as CHOP, that participates in regulating IOP (Table 3) and was down regulated in both WTβ3(FC = −2.90) and CAβ3 (FC = −2.25) cells.
In contrast to the numerous genes whose expression was significantly altered by constitutive activation of αvβ3 integrin in CAβ3 cells, the expression of very few transcription factors was significantly altered by overexpression of WT αvβ3 integrin in WTβ3 cells. Of those that met the FC cutoff, only ZNF626 (FC = 1.93) was upregulated while RUNX1T1 (FC = −2.51), NUPR1 (FC = −2.47) and EGR2 (FC = −1.52) were downregulated. This suggests that overexpression of WT αvβ3 integrin by itself did not have a major impact on the transcriptional activity of TM cells. Table 2. Differentially expressed genes in WTβ3 and/or CAβ3 cells that are associated with glaucoma. Except where noted by *, Log 2 FC values were reported in DGE analysis using EdgeR. * Genes that were filtered out due to low abundance did not originally have a FC determined (ND). Log 2 FC was calculated from the raw count data. The values in parentheses were originally calculated using EdgeR.

Hierarchical Clustering Reveals Differences in the Expression of ECM Protein Genes
Another group of proteins that appeared to be differentially affected by expression of αvβ3 integrin and its activation were those involved in the organization of the ECM ( Figure 5). Interestingly, there were only six core ECM proteins that showed differential expression in either WTβ3 or CAβ3 cells compared to EV cells that met the ±1.5 Log 2 FC cutoff. Some of these genes are also linked to glaucoma ( Table 2). Of those, 4/6 were significantly altered only in WTβ3 cells while 2/6 were significantly altered only in CAβ3 cells. Hierarchical clustering showed that overexpression of WT αvβ3 integrins caused an upregulation in COL8A2 (FC = 2.71), EMILIN1(FC = 1.51) and POSTN (FC = 3.88). Both COL8A2 and POSTN have been reported to be associated with glaucoma (Table 2). POSTN plays a role in ASD and disruption in its expression leads to a disorganization of collagen fibrils in the TM. Additionally, its expression has also been associated with several fibrotic diseases [73]. COL8A2 mutations have been found in some glaucoma patients [74], and its expression is elevated in HTM cells in response to dexamethasone suggesting a role in GC-induced glaucoma as well as in the development of the anterior segment [75]. In contrast, COL11A1 (FC = −1.96), which is associated with both POAG and primary angle closure glaucoma (PACG) ( Table 2) was downregulated in WTβ3cells. Regarding core ECM genes that were significantly altered only in CAβ3 cells, neither VWA5A (FC = 2.67) nor COL13A1 (FC = −1.64) have been found to be associated with glaucoma, although COL13A1 was differentially regulated in response to the steroid triamcinolone in cultured HTM cells [61].
In addition to these core ECM protein genes, we noted the differential expression of four matricellular genes. SPOCK1 (FC = 1.93), TLL1 FC = 4.19) and SULF1 (FC = 1.76) are all uniquely upregulated in CAβ3 cells. SPOCK1 encodes a member of the SPARC family and is associated with glaucoma ( Table 2). TLL1 and SULF1 have no reported associations with glaucoma or IOP regulation, however, both genes encode proteins associated with fibrosis and an epithelial to mesenchymal transition (EMT) [76] which shows similar phenotypic changes to that observed in POAG. In contrast, TNC was downregulated in both cell lines. However, it exhibited an~2-fold lower expression level in CAβ3 cells (FC = −3.59) compared to WTβ3 cells (FC = −2.25). Relative to the organization and remodeling of the ECM in the TM and in glaucoma [60], we observed a large decrease in MMP1 expression in both cell lines (CAβ3 FC = −6.35, WTβ3 FC = −3.59). TIMP3 was also downregulated in both cell lines (CAβ3 FC = −2.02, WTβ3 FC = −1.57).

Hierarchical Clustering Reveals Differences in the Expression of Genes Associated with Cell Adhesion
In addition to the differential expression of ECM genes, multiple genes for proteins within the cadherin (CAD) family of cell adhesion receptors met the ±1.5 Log 2 FC cutoff. A total of 16 genes in this category were identified in either the WTβ3 or CAβ3 cell lines along with other genes associated with cell adhesion (Figure 6). Among those genes were two classical cadherin genes, CDH6 (K-cadherin) and CDH11. CDH6 was uniquely upregulated in CAβ3 cells (FC = 1.97) while CDH11 was uniquely downregulated in WTβ3 cells (FC = −3.12). CDH6 is downregulated in bovine TM cells in response to dexamethasone treatment [77], and both CDH6 and CDH11 have been suggested to play a role in regulating IOP (Table 3).
Cells 2021, 10, x FOR PEER REVIEW 11 of 34 Figure 5. Hierarchical clustering of ECM-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log2 FC ± 1.5 cutoff and were classified as encoding either core ECM proteins, ECM-associated proteins or matricellular proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4.
In addition to these core ECM protein genes, we noted the differential expression of four matricellular genes. SPOCK1 (FC = 1.93), TLL1 FC = 4.19) and SULF1 (FC = 1.76) are all uniquely upregulated in CAβ3 cells. SPOCK1 encodes a member of the SPARC family and is associated with glaucoma ( Table 2). TLL1 and SULF1 have no reported associations with glaucoma or IOP regulation, however, both genes encode proteins associated with fibrosis and an epithelial to mesenchymal transition (EMT) [76] which shows similar phenotypic changes to that observed in POAG. In contrast, TNC was downregulated in both cell lines. However, it exhibited an ~ 2-fold lower expression level in CAβ3 cells (FC = −3.59) compared to WTβ3 cells (FC = −2.25). Relative to the organization and remodeling of the ECM in the TM and in glaucoma [60], we observed a large decrease in MMP1 expression in both cell lines (CAβ3 FC = −6.35, WTβ3 FC = −3.59). TIMP3 was also downregulated in both cell lines (CAβ3 FC = −2.02, WTβ3 FC = −1.57).

Hierarchical Clustering Reveals Differences in the Expression of Genes Associated with Cell Adhesion
In addition to the differential expression of ECM genes, multiple genes for proteins within the cadherin (CAD) family of cell adhesion receptors met the ± 1.5 Log2 FC cutoff. A total of 16 genes in this category were identified in either the WTβ3 or CAβ3 cell lines along with other genes associated with cell adhesion (Figure 6). Among those genes were two classical cadherin genes, CDH6 (K-cadherin) and CDH11. CDH6 was uniquely upregulated in CAβ3 cells (FC = 1.97) while CDH11 was uniquely downregulated in WTβ3 cells (FC = −3.12). CDH6 is downregulated in bovine TM cells in response to dexamethasone treatment [77], and both CDH6 and CDH11 have been suggested to play a role in regulating IOP (Table 3).

Figure 5.
Hierarchical clustering of ECM-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff and were classified as encoding either core ECM proteins, ECM-associated proteins or matricellular proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for  (Table S5) to a similar extent except for PCDHB10 which was expressed at a higher level in CAβ3 cells (FC = 3.39) than in WTβ3 cells (FC = 1.93). The significance of these genes in the pathogenesis of POAG is unclear since the functions of many of these protocadherins are still being unraveled.
Other classes of adhesion molecules demonstrated differentially expressed genes besides members of the CAD family including members of the semaphorin family, the IgG superfamily (IGSF), the TIE receptor family, a CD34 family member and an axon guidance molecule. Of interest is SEMA6A (FC = 1.55) which was uniquely upregulated only in CAβ3 cells and is also found on the glaucoma locus GLC1G (Table 4). NCAM2 is also of interest as it has been found to be associated with POAG (Table 2) and is uniquely upregulated in CAβ3 cells (FC = 1.78). RGMA encodes a member of the repulsive guidance molecule family that regulates axon guidance and can act as a co-receptor for TGFβ signaling [78]. It is only upregulated in CAβ3 cells (FC = 2.16), and elevated RGMA expression has been associated with glaucoma in mice [79]. Cells 2021, 10, x FOR PEER REVIEW 12 of 34 Figure 6. Hierarchical clustering of adhesion-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to EV control cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log2 FC ± 1.5 cutoff and were classified as encoding either adhesion proteins or adhesion-associated proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for  (Table S5) to a similar extent except for PCDHB10 which was expressed at a higher level in CAβ3 cells (FC = 3.39) than in WTβ3 cells (FC = 1.93). The significance of these genes in the pathogenesis of POAG is unclear since the functions of many of these protocadherins are still being unraveled.
Other classes of adhesion molecules demonstrated differentially expressed genes besides members of the CAD family including members of the semaphorin family, the IgG superfamily (IGSF), the TIE receptor family, a CD34 family member and an axon guidance molecule. Of interest is SEMA6A (FC = 1.55) which was uniquely upregulated only in CAβ3 cells and is also found on the glaucoma locus GLC1G (Table 4). NCAM2 is also of interest as it has been found to be associated with POAG (Table 2) and is uniquely upregulated in CAβ3 cells (FC = 1.78). RGMA encodes a member of the repulsive guidance molecule family that regulates axon guidance and can act as a co-receptor for TGFβ signaling [78]. It is only upregulated in CAβ3 cells (FC = 2.16), and elevated RGMA expression has been associated with glaucoma in mice [79]. Figure 6. Hierarchical clustering of adhesion-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to EV control cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff and were classified as encoding either adhesion proteins or adhesion-associated proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4. Table 4. Differentially expressed genes in WTβ3 and/or CAβ3 cells that are located on a glaucoma locus. Except where noted by *, Log 2 FC values were reported in DGE analysis using EdgeR. * Genes that were filtered out due to low abundance did not originally have a FC determined (ND). Log 2 FC was calculated from the raw count data. The values in parentheses were originally calculated using EdgeR. Several genes that regulate adhesion were also found to be differentially expressed in both cell lines. ANGPT1 is only upregulated in CAβ3 cells (FC = 1.76), and, in addition to having a role in Schlemm canal development [82], it has been mapped to the CLC1D glaucoma locus (Table 4) (Table 2), the proteins encoded by these genes have not been reported to be associated with glaucoma. The BMPER protein, however, has been suggested to play a role in TM function under both normal and glaucomatous conditions as it antagonizes BMP signaling which, in turn, promotes TGFβ2 signaling (Table 3).

Hierarchical Clustering Reveals Differences in Genes Encoding Cytokines, Cytokine Regulators and Genes Involved in TGFβ, BMP and WNT Signaling
Another group of genes whose expression was significantly affected by overexpression or activation of αvβ3 was cytokines, cytokine receptors and cytokine-associated genes (Figure 7). Since it is well established that TGFβ signaling plays a significant role in both glaucoma and IOP regulation [2] it was noteworthy that hierarchical clustering revealed that variable differences in the expression of genes associated with TGFβ signaling were observed between the two cell lines (see WTβ3-CAβ3 column). As shown in Figure 7 Figure 7. Hierarchical clustering of cytokines and cytokine-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log2 FC ± 1.5 cutoff and were classified as encoding either cytokines or cytokine-associated proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4.
Hierarchical clustering also showed several Wnt-related genes differentially expressed in the two cell lines. (Figure 7, clusters F, J). WNT8B was uniquely upregulated in CAβ3 cells (FC = 1.90) while WNT7B (FC = -1.82) and DKK1 (FC = -1.91) were downregulated in CAβ3 cells. Although Wnt signaling has been implicated in IOP regulation and glaucoma [64], none of these genes are known to be associated with glaucoma. The gene Figure 7. Hierarchical clustering of cytokines and cytokine-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff and were classified as encoding either cytokines or cytokine-associated proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Hierarchical clustering also showed several Wnt-related genes differentially expressed in the two cell lines. (Figure 7, clusters F, J). WNT8B was uniquely upregulated in CAβ3 cells (FC = 1.90) while WNT7B (FC = −1.82) and DKK1 (FC = −1.91) were downregulated in CAβ3 cells. Although Wnt signaling has been implicated in IOP regulation and glaucoma [64], none of these genes are known to be associated with glaucoma. The gene DCDC2 (cluster F) was also uniquely upregulated in CAβ3 cells (FC = 2.025). Although not a member of the Wnt family, the DCDC2 protein has been reported to inhibit β-catenin-dependent Wnt signaling [83].
Among the other genes in this group that met the ±1.5 Log 2 FC cutoff were five members of the CXC chemokine family that were downregulated only in the CAβ3 cell line (Figure 7 Table 2). PENK (Figure 7, cluster P) which is localized to the GLC1D glaucoma locus (Table 4)

Hierarchical Clustering Reveals Differences in Genes Encoding GTP-Binding and Cytoskeleton-Associated Proteins
It is well established that integrin signaling plays a major role in regulating GTPase activity and the cytoskeleton [85]. Additionally, the significance of the Rho GTPase and its downstream effectors in regulating aqueous humor outflow and the pathogenesis of POAG has been well established [86]. Intriguingly many of the genes (15/26) that were identified in this group (Figure 8) are regulated by proteins associated with or regulated by TGFβ and/or BMP signaling (DIRAS3, NHS, ARHGAP24, MYH15, MCF2L, RAPGEF4, GAS7, ARHGAP22, SRGAP3, IQGAP2, SORBS1, RUNDC3A, AGAP2, INA and MAP7) (Tables S7 and S9) while 4/26 (ADAP2, RHOD, MYBPH and S1PR5) are not. As with the other categories of genes affected by overexpression and/or activation of αvβ3 integrin, hierarchical clustering of these genes that met the ±1.5 Log 2 FC cutoff revealed differences in gene expression between the WTβ3 and CAβ3 cell lines (Figure 8, clusters A-C and I-L).
identified in this group (Figure 8) are regulated by proteins associated with or regulated by TGFβ and/or BMP signaling (DIRAS3, NHS, ARHGAP24, MYH15, MCF2L, RAPGEF4, GAS7, ARHGAP22, SRGAP3, IQGAP2, SORBS1, RUNDC3A, AGAP2, INA and MAP7) (Tables S7 and S9) while 4/26 (ADAP2, RHOD, MYBPH and S1PR5) are not. As with the other categories of genes affected by overexpression and/or activation of αvβ3 integrin, hierarchical clustering of these genes that met the ± 1.5 Log2 FC cutoff revealed differences in gene expression between the WTβ3 and CAβ3 cell lines (Figure 8, clusters A-C and I-L). gure 8. Hierarchical clustering of GTP-binding and cytoskeleton-associated genes finds differential gene expression in Tβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for E that met the Log2 FC ± 1.5 cutoff and were classified as encoding GTP-binding and/or cytoskeleton-associated prons were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by ters on the right side of the map were determined manually. Genes within each of these clusters are listed in the disayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4.
In clusters A and B, five genes were uniquely upregulated in WTβ3 cells relative to EV control cells. In this group we see several members of the Ras GTPase family (ADAP2, DIRAS3 and RND2) and proteins that control cytoskeleton remodeling (NHS, and TUBA4A.) (Table S3). Of those genes, DIRAS3 (FC = 1.935) is the only gene with any association with glaucoma or IOP regulation [87].
In contrast to the genes in clusters A and B, the genes in cluster C and D were all downregulated in CAβ3 cells relative to WTβ3 cells. Of the three genes found in cluster C, only ARHGAP24 (cluster C, FC = −2.04), which encodes a protein that participates in Figure 8. Hierarchical clustering of GTP-binding and cytoskeleton-associated genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff and were classified as encoding GTP-binding and/or cytoskeleton-associated proteins were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4.
In clusters A and B, five genes were uniquely upregulated in WTβ3 cells relative to EV control cells. In this group we see several members of the Ras GTPase family (ADAP2, DIRAS3 and RND2) and proteins that control cytoskeleton remodeling (NHS, and TUBA4A.) (Table S3). Of those genes, DIRAS3 (FC = 1.935) is the only gene with any association with glaucoma or IOP regulation [87].
In contrast to the genes in clusters A and B, the genes in cluster C and D were all downregulated in CAβ3 cells relative to WTβ3 cells. Of the three genes found in cluster C, only ARHGAP24 (cluster C, FC = −2.04), which encodes a protein that participates in cell-ECM adhesion complexes in cultured HTM cells [88], has any relevance to TM function or glaucoma as it has also been identified in cultured ONH astrocytes as being phosphorylated in response to increases in hydrostatic pressure [89]. LPXN (FC = −1.51; cluster D) encodes a focal adhesion protein that can also function as a transcription factor [90]. It is a member of the paxillin family along with HIC5 which has been shown to play a role in TGF-β2-induced fibrogenic activity and dexamethasone-induced MYOC expression in HTM cells [91].
All six genes in clusters I-L were upregulated only in CAβ3 cells. Of these genes, however, only AGAP2 (FC = 1.77) is likely to play an important role in regulating TM function as it relates to IOP and POAG pathogenesis. AGAP2 is a TGFβ-responsive gene that encodes a positive regulator of TGFβ signaling as it relates to fibrosis [92].

Hierarchical Clustering of Miscellaneous and Unclassified Genes
Hierarchical clustering was also performed on the remaining genes that were identified using the ±1.5 Log 2 FC cutoff ( Figure 9). As with the other hierarchical clustering maps, we observed genes that exhibited clear differences in expression between WTβ3 and CAβ3 cells (clusters A, H, I, J, M and O). Of the genes found within these clusters, however, only five have any reported association with IOP regulation or glaucoma. In cluster O, CYP1B1 was uniquely downregulated in WTβ3 cells (FC = −2.28). CYP1B1 is associated with glaucoma (Table 2), and mutations in this gene in mice are associated with increases in IOP, oxidative stress and TM dysgenesis [55]. Four other genes worth noting, TXNIP, CA9, RBP7 and MBP are found within clusters H, I and J, respectively. TXNIP expression was uniquely upregulated in CAβ3 cells (FC = 3.64). It is part of the thioredoxin (TRX) system that protects cells from damage due to oxidative stress and functions in suppressing EMT induced by high glucose [93] or TGFβ (Table S9). Its expression has not been studied in the TM, but TXNIP is upregulated in retinal ganglion cells in response to experimentally induced increases in IOP [94]. CA9 (FC = 2.43), which was upregulated only in WTβ3 cells, encodes for carbonic anhydrase 9 which appears to play a role in IOP regulation rather than in the pathogenesis of POAG [95]. RBP7 which is located on the glaucoma locus GLC3B (Table 4) was uniquely upregulated in WTβ3 cells (Log 2 FC = 1.925). MBP was uniquely downregulated in CAβ3 cells (FC = −2.62). Elevated MBP serum levels as well as autoantibodies against the protein have been reported in patients with POAG [96]. . Hierarchical clustering of miscellaneous and unclassified genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Miscellaneous and unclassified genes from each of the three contrasts analyzed for DGE that met the Log2 FC ± 1.5 cutoff were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for Figure 4.

RT-qPCR Validation of Changes in Gene Expression
We then used RT-qPCR to validate selected genes that were altered in these categories. Table 5 shows the primer sequences used for the analysis. Overall, the data obtained through RNA-Seq analysis and RT-qPCR analysis were consistent with each other. As shown in Figure 10, RT-qPCR showed that TGFB2 and ADCY1 were both expressed at significantly higher levels in both WTβ3 (TGFB2, p < 0.0035; ADCY1, p < 0.002) and CAβ3 cells (TGFB2, p < 0.0015; ADCY1, p < 0.002) relative to control EV cells. Analysis of expression of CDH6 (p < 0.015), SPOCK1 (p < 0.002) and PAX6 (p < 0.0015) found that all three genes were significantly upregulated in CAβ3 cells relative to both EV and WTβ3 cells as predicted by the RNA-Seq studies. As predicted by RNA-Seq analysis, GREM1 (both, p < Figure 9. Hierarchical clustering of miscellaneous and unclassified genes finds differential gene expression in WTβ3 and CAβ3 cells relative to control EV cells and to each other. Miscellaneous and unclassified genes from each of the three contrasts analyzed for DGE that met the Log 2 FC ± 1.5 cutoff were analyzed as in Figure 4. A heat map displaying standardized z-scores was generated. The clusters identified by letters on the right side of the map were determined manually. Genes within each of these clusters are listed in the displayed table. Those genes highlighted by an asterisk represent genes that had been dealt with as described for

RT-qPCR Validation of Changes in Gene Expression
We then used RT-qPCR to validate selected genes that were altered in these categories. Table 5 shows the primer sequences used for the analysis. Overall, the data obtained through RNA-Seq analysis and RT-qPCR analysis were consistent with each other. As shown in Figure 10, RT-qPCR showed that TGFB2 and ADCY1 were both expressed at significantly higher levels in both WTβ3 (TGFB2, p < 0.0035; ADCY1, p < 0.002) and CAβ3 cells (TGFB2, p < 0.0015; ADCY1, p < 0.002) relative to control EV cells. Analysis of expression of CDH6 (p < 0.015), SPOCK1 (p < 0.002) and PAX6 (p < 0.0015) found that all three genes were significantly upregulated in CAβ3 cells relative to both EV and WTβ3 cells as predicted by the RNA-Seq studies. As predicted by RNA-Seq analysis, GREM1 (both, p < 0.0015), MMP1 (both, p < 0.002), DKK1 (both, p < 0.0015), BMPER (WTβ3, p < 0.002, CAβ3, p < 0.0015), DDIT3 (both, p < 0.0011) and PENK (both, p < 0.0012) were all downregulated in both cell lines compared to EV cells. Finally, RT-qPCR found VEGFA expression to be significantly downregulated in CAβ3 cells relative to both EV (p < 0.0015) and WTβ3 cells (p < 0.002) which is consistent with RNA-Seq analysis. However, by RT-qPCR there was no significant difference in VEGFA expression between EV and WTβ3 even though the RNA-Seq analysis found a statistically significant downregulation of this gene in WTβ3 cells relative to EV cells (FC = −0.325, p < 0.006). Table 5. Primer sequences used for RT-qPCR validation of select genes identified by RNA-Seq analysis. All sequences are for the human version of each gene.
Type Sequence We also sought to determine if, in addition to the upregulation of TGFB2 mRNA expression observed in CAβ3 cells relative to EV cells, there was increased production of TGFβ2 protein by CAβ3 cells. As shown in Figure 11, ELISA analysis indicated that there was significantly more TGFβ2 in CAβ3 conditioned medium compared to EV conditioned medium (p < 0.05).   We also sought to determine if, in addition to the upregulation of TGFB2 mRNA expression observed in CAβ3 cells relative to EV cells, there was increased production of TGFβ2 protein by CAβ3 cells. As shown in Figure 11, ELISA analysis indicated that there was significantly more TGFβ2 in CAβ3 conditioned medium compared to EV conditioned medium (p < 0.05). Figure 11. CAβ3 cells secrete more TGFβ2 than EV cells. ELISA data for total TGFβ2 in the medium of EV and CAβ3 cells. TGFβ2 was significantly (* p < 0.05) increased in the medium of CAβ3 cells compared to EV cells. Figure 11. CAβ3 cells secrete more TGFβ2 than EV cells. ELISA data for total TGFβ2 in the medium of EV and CAβ3 cells. TGFβ2 was significantly (* p < 0.05) increased in the medium of CAβ3 cells compared to EV cells.

Panther Pathway Analysis Finds Gene Overrepresentation Only in CAβ3 Cells
We next determined which signaling pathways were enriched with changes in gene expression. For each cell line, genes identified as being upregulated or downregulated using the Log 2 ± 1.5 FC cutoff (Figure 2, Tables S1-S6) were entered into the PANTHER pathway database. Although no pathways were found to be statistically overrepresented in WTβ3 cells, three pathways were significantly (p < 0.05) overrepresented in CAβ3 cells ( Table 6). The pathways overrepresented in CAβ3 cells are involved in TGFβ2, Wnt and cadherin signaling. This was a striking finding given the roles that each of these pathways, either individually or combined, appear to play in both normal IOP regulation and POAG pathogenesis [64]. These results support the idea that it is the activation, rather than simple overexpression, of the αvβ3 integrin that is important in significantly altering gene expression in TM cells. Table 6. Pathways enriched by activation of αvβ3 integrin compared to control cells. Genes identified as being upregulated or downregulated ± 1.5 Log 2 FC were entered into the PANTHER database (pantherdb.org) and analyzed using a statistical overrepresentation test (Fisher's Exact test). A total of 422 genes were entered. 365 genes were successfully mapped to a pathway. 57 genes were left unmapped. The genes associated with each pathway are as follows: TGFβ2 signaling (TGFB2, BMP3,  BMP4, SMAD6, TLL1,   Since all three of these pathways have been implicated in IOP regulation and/or POAG, we examined which of the genes found to be differentially expressed in response to αvβ3 integrin activation could be targeting one or more of these pathways. As shown in Table 7, we identified genes encoding seven transcription factors (DDIT3, ISL1, KLF4, KLF2, MYOCD, SOX5 and STAT5A), three ECM-associated proteins (SULF1, TLL1 and TNC), two GTP-binding or cytoskeleton proteins (AGAP2, LPXN), two adhesion-associated proteins (BMPER and PCDH1) and three cytokines (BMP4, GREM1 and IL18) that have been reported to target the TGFβ2, BMP and/or Wnt signaling pathways. Additionally, we identified many genes from these classes that were targets for these same signaling pathways (Table S9)

PPI Network Analysis Shows Enrichment Only in CAβ3 Cells
To further examine the relationship between the proteins expressed by the genes that met the ±1.5 Log 2 FC cutoff for each cell line, we submitted the gene lists from each cell line for analysis into the STRING database in order to generate protein-protein interaction (PPI) networks. Input lists consisted of the uniquely altered genes from each cell line as well as genes that were shared by both cell lines. Surprisingly, analysis of the PPI networks generated using the genes from the WTβ3 cells found only a modest number of interactions between the genes ( Figure 12A,B) and statistical analysis of the interactions at two different confidence levels (p = 0.143, medium confidence; p = 0.259, high confidence) found that these groups of interactions were no different from those of a random collection of input genes. Thus, the limited number of potential interactions in WTβ3 cells was unlikely to result in significant changes to the signaling network in these cells. In contrast, analysis of the networks calculated using the input genes from CAβ3 cells showed many possible interactions at both confidence levels (Figures 13 and 14). Statistical analysis of the CAβ3 networks at both medium and high confidence levels resulted in PPI enrichment p-values of 1 × 10 −16 indicating that the interactions of the gene products from the CAβ3 cells were significantly more than those expected for a random collection of input genes. This suggests that activation of αvβ3 integrin results in a high number of both physical and functional interactions between the proteins expressed by the input genes that could result in significant changes to the signaling network in these cells.  Figure 12B: permanent link: https://version-11-0b.string-db.org/cgi/network?networkId=b9wf05GlAE9T) (accessed on 20 June 2021). The confidence score is the approximate probability that a predicted link exists between two proteins in the same metabolic map in the KEGG database. Both PPI enrichment scores indicate that the calculated networks are not significantly different from networks generated by a random collection of input genes at the same confidence level. For both networks, the proteins for which there were no connections to be mapped were not displayed. The location of β3 integrin (INTB3), members of the TGFβ2, Wnt and Cadherin signaling pathways as well as those gene products that are also differentially expressed in CAβ3 cells are all indicated. * Genes that were validated by RT-qPCR.  High confidence CAβ3 PPI network shows an extensive network of potential interactions. Differentially expressed genes resulting from overexpression of a CA-αvβ3 integrin were entered into the STRING database. The input consisted of uniquely altered genes as well as altered genes that were shared with cells overexpressing a WT-αvβ3 integrin. The calculated network (statistical confidence = 0.7) has a PPI enrichment value = 1 × 10 −16 indicating that the network has significantly more interactions than expected for a random collection of input genes at the same confidence level. The We then took the list of altered genes in CAβ3 cells that PANTHER identified as being overrepresented for the TGFβ2, Wnt and Cadherin pathways and generated medium High confidence CAβ3 PPI network shows an extensive network of potential interactions. Differentially expressed genes resulting from overexpression of a CA-αvβ3 integrin were entered into the STRING database. The input consisted of uniquely altered genes as well as altered genes that were shared with cells overexpressing a WT-αvβ3 integrin. The calculated network (statistical confidence = 0.7) has a PPI enrichment value = 1 × 10 −16 indicating that the network has significantly more interactions than expected for a random collection of input genes at the same confidence level. The location of β3 integrin (INTB3), members of the TGFβ2, Wnt and Cadherin (CAD) signaling pathways as well as those gene products that are also differentially expressed in WTβ3 cells are all indicated. * Genes that were validated by RT-qPCR. The proteins for which there were no connections to be mapped where not displayed. (Interactive link for Figure 14: permanent link: https://version-11-0b.string-db.org/cgi/network?networkId=bZECX4AuPQ32) (accessed on 20 June 2021).
We then took the list of altered genes in CAβ3 cells that PANTHER identified as being overrepresented for the TGFβ2, Wnt and Cadherin pathways and generated medium and high confidence PPI networks ( Figure 15A,B). In addition to the genes identified by PANTHER, we included four genes that demonstrated significantly altered expression in CAβ3 cells, DKK1 [114], BMPER [70], RGMA [78] and GREM1 [69], that have also been connected to one or more of these three signaling pathways. Despite consisting of only 31 genes in total, the medium and high confidence networks that were generated had PPI enrichment p-values of 1 × 10 −16 and 5.72 × 10 −12 , respectively, suggesting that these networks could be significantly altered in CAβ3 cells. and high confidence PPI networks ( Figure 15A,B). In addition to the genes identified by PANTHER, we included four genes that demonstrated significantly altered expression in CAβ3 cells, DKK1 [114], BMPER [70], RGMA [78] and GREM1 [69], that have also been connected to one or more of these three signaling pathways. Despite consisting of only 31 genes in total, the medium and high confidence networks that were generated had PPI enrichment p-values of 1 × 10 −16 and 5.72 × 10 −12 , respectively, suggesting that these networks could be significantly altered in CAβ3 cells. Figure 15. PPI networks consisting only of differentially expressed genes from CAβ3 cells that were overrepresented in the TGFβ2, Wnt and CAD signaling pathways demonstrate a significant number of potential interactions. Genes from Table 6 along with four additional genes (DKK1, BMPER, RGMA and GREM1) that have been shown to participate in one or more of these pathways, were entered into the STRING database. Both PPI enrichment scores indicate that these networks have significantly more interactions than expected for a random collection of input genes at the same confidence level. The proteins for which there were no connections to be mapped are also displayed.

Discussion
In this paper we show that activation of αvβ3 integrin signaling in trabecular meshwork cells causes a differential expression of multiple genes that contribute to TGFβ2 signaling and fibrotic-like changes that could lead to the development of glaucoma [2]. Focusing on genes that had altered expression of ± 1.5 Log2 fold change (FC), we found that overexpressing a CA-αvβ3 integrin resulted in the differential expression of genes compared to cells overexpressing a WT-αvβ3 integrin. In particular, we found an enrichment for genes encoding proteins found in the TGFβ2, Wnt and cadherin signaling pathways all of which have been implicated in the pathogenesis of POAG [2,64]. Among the significant changes detected was an upregulation in TGFB2 which is found to be elevated in 50% of patients with POAG. CDH6, BMP4, MMP1 and TIMP3 are other genes found to be differentially expressed and associated with TGFβ2 signaling, remodeling of the TM ECM and/or increased IOP (Table 2). In contrast, no enrichment for genes in any specific signaling pathway was detected in TM-1 cells overexpressing WT-αvβ3 integrin. Many of the differentially expressed genes detected belong to several classes that may contribute to the development of POAG. These include genes encoding transcription factors, adhesion and cytoskeleton proteins, ECM proteins, cytokines and GTPases. These findings, together with our earlier studies [4,6,11,41], support the idea that changes in αvβ3 integrin signaling are likely to contribute to the changes in trabecular meshwork associated with glaucoma. Figure 15. PPI networks consisting only of differentially expressed genes from CAβ3 cells that were overrepresented in the TGFβ2, Wnt and CAD signaling pathways demonstrate a significant number of potential interactions. Genes from Table 6 along with four additional genes (DKK1, BMPER, RGMA and GREM1) that have been shown to participate in one or more of these pathways, were entered into the STRING database. (A) Network calculated using a medium (0.4) statistical confidence; PPI enrichment p-value = 1 × 10 −16 . (Interactive link for Figure 15A: permanent link: https://version-11-0b.string-db.org/cgi/network?networkId=bxHIv7XuhKtR) (accessed on 20 June 2021). (B) Network calculated using a high (0.7) statistical confidence; PPI enrichment p-value = 5.72 × 10 −12 . (Interactive link for Figure 15A: permanent link: https://version-11-0b.string-db.org/cgi/network?networkId=bWKleDuHyVqt) (accessed on 20 June 2021). Both PPI enrichment scores indicate that these networks have significantly more interactions than expected for a random collection of input genes at the same confidence level. The proteins for which there were no connections to be mapped are also displayed.

Discussion
In this paper we show that activation of αvβ3 integrin signaling in trabecular meshwork cells causes a differential expression of multiple genes that contribute to TGFβ2 signaling and fibrotic-like changes that could lead to the development of glaucoma [2]. Focusing on genes that had altered expression of ± 1.5 Log 2 fold change (FC), we found that overexpressing a CA-αvβ3 integrin resulted in the differential expression of genes compared to cells overexpressing a WT-αvβ3 integrin. In particular, we found an enrichment for genes encoding proteins found in the TGFβ2, Wnt and cadherin signaling pathways all of which have been implicated in the pathogenesis of POAG [2,64]. Among the significant changes detected was an upregulation in TGFB2 which is found to be elevated in 50% of patients with POAG. CDH6, BMP4, MMP1 and TIMP3 are other genes found to be differentially expressed and associated with TGFβ2 signaling, remodeling of the TM ECM and/or increased IOP (Table 2). In contrast, no enrichment for genes in any specific signaling pathway was detected in TM-1 cells overexpressing WT-αvβ3 integrin. Many of the differentially expressed genes detected belong to several classes that may contribute to the development of POAG. These include genes encoding transcription factors, adhesion and cytoskeleton proteins, ECM proteins, cytokines and GTPases. These findings, together with our earlier studies [4,6,11,41], support the idea that changes in αvβ3 integrin signaling are likely to contribute to the changes in trabecular meshwork associated with glaucoma.
One of the most striking changes observed in this study was the increase in TGFB2 expression along with the expression of genes associated with TGFβ2 signaling. TGFβ2 levels are known to be elevated in the aqueous humor of POAG patients and increased expression of TGFβ2 causes an elevation in IOP in rodent models of ocular hypertension [2]. However, it is not known why the levels are upregulated. The upregulation of TGFB2 expression in response to αvβ3 integrin signaling in our study would suggest that increases in αvβ3 integrin signaling may be involved. This is not entirely surprising. αvβ3 integrin has been shown to be associated with the upregulation of TGFB2 in several cell types including a breast cancer epithelial cell line [115] and human skin fibroblasts [116].
Recent studies have also highlighted the importance of αvβ3 integrin in TGFβ2mediated senescence [117] further suggesting that there is a link between αvβ3 integrin activity and TGFβ2 expression. Interestingly, the expression of αvβ3 integrin in human skin fibroblasts and in the liver of C57BL/6J mice correlated with aging [116] and senescence [117]. Elevated αvβ3 integrin levels have also been found to be associated with apoptosis [118]. This raises the interesting question of whether increases in αvβ3 integrin expression in the aged TM could also be contributing to cell loss or senescence in the TM which is known to be associated with POAG and thought to be one of the causes of POAG [119] Notably, Yu, et al. [120] found that TGFβ2 induced senescence-associated changes within HTM cells and activated a senescence-related signaling pathway similar to that reported by Rapisarda, et al. [116]. This suggests an additional mechanism whereby dysregulated αvβ3 integrin signaling together with increased TGFB2 expression could contribute to POAG pathogenesis.
In addition to TGFB2, many of the genes that were altered belonged to classes of proteins such as ECM proteins, cell adhesion molecules and transcription factors, that are known to be altered in glaucoma and either regulate TGFβ2 signaling (Table 7) or are regulated by TGFβ2 signaling (Table S9). The major gene class that was affected by activation of αvβ3 integrin signaling by far, however, was transcription-associated genes. Several of these genes are known to be associated with glaucoma (Table 2). For instance, PAX6 is considered a master regulator of TGFβ2 expression in ocular tissues [104]. FOSL1 [46], SOX5, STAT5A, MYOCD, ISL1, KLF4 and KLF2 all play a role in fibrosis, myofibroblast differentiation, EMT, TGFβ signaling, BMP signaling, and/or Wnt signaling ( Table 7). The significant up and down regulation of these types of transcription factors suggest that αvβ3 integrin signaling is altering multiple signaling pathways involved in the fibrogenic transition of the TM during glaucoma. The significance of these changes is unknown, however, the downregulation of both KLF2 and KLF4 expression would potentially attenuate two significant inhibitors of TGFβ2 signaling and thus enhance TGFβ-induced EMT [47,51].
The finding that overexpression of a constitutively active αvβ3 integrin also led to an enrichment in genes associated with TGFβ, Wnt and cadherin signaling pathways further supports the hypothesis that αvβ3 integrin is playing a role in the transcriptional regulation of the biological processes associated with glaucoma. In contrast, overexpression of wild type αvβ3 integrin did not have a significant effect on any of these signaling pathways suggesting that it is the active signaling conformation of αvβ3 integrin which is more important than its expression levels or adhesion to the ECM.
It is noteworthy that several other genes that were upregulated by αvβ3 integrin signaling are associated with various types of glaucoma and IOP regulation (Tables 2-4) further supporting the idea that αvβ3 integrin and the concomitant increase in TGFB2 expression together may be playing a role in glaucoma. Among the genes upregulated and associated with glaucoma were several components of the ECM (POSTN, COL8A2, SEMA6A) as well as ANGPT1 and BMP4. We also found a decrease in expression of several genes associated with the ECM especially MMP1 and TIMP3. These genes all have one thing in common and that is their expression would affect the remodeling of ECM.
We also saw an effect on genes associated with the development of the anterior segment. For instance, PAX6 and LMX1B are both associated with ASD which, in turn, is associated with an increased risk of glaucoma [53,54,121]. Although it does not encode a transcription factor, alterations in CYP1B1 expression have also been shown to cause ASD [55] and are associated with glaucoma. LMX1B, independent of its ASD association, and GAS7 are also associated with elevated IOP and POAG [54,56,121,122]. Thus, it is interesting that these genes exhibited DGE when αvβ3 integrin is overexpressed and/or activated. This suggests that αvβ3 integrin may play a role in TM development as well as glaucoma.
One surprising finding was the downregulation of several pro-inflammatory cytokines especially several members of the CXC family of chemokines, several members of the interleukin family and VEGFA. Although it is not clear why this would happen, the upregulation of these cytokines is usually associated with acute early inflammatory response [123] and their levels tend to wain during the later fibrotic phase when TGFB2 is upregulated. Hence, downregulation of chemokine expression along with elevated levels of TGFβ2 could be important determinants in the pathogenesis of fibrosis [124].
Finally, previous studies in our lab have also found that two phenotypes associated with glaucoma, cytoskeletal rearrangements into CLANs and decreased phagocytosis, can be triggered in response to activation of αvβ3 signaling. Both phenotypes are dependent upon signaling pathways regulated by the Rac GTPase [8][9][10]125,126]. Several genes (ARHGAP22, ARHGAP24, SRGAP3 and IQGAP2) that regulate Rac activity were also differentially expressed in CAβ3 cells and are targets of the TGFβ, BMP and/or Wnt signaling pathways (Table S9). Whether these genes play a role in CLAN formation or phagocytosis remains to be determined, however, it does further support a role for αvβ3 integrin signaling in the development of these changes in glaucoma.
In summary, to the best of our knowledge this study is the first to show that increased αvβ3 integrin expression and signaling could be responsible for enhanced TGFβ2 expression in POAG. It also shows that changes in the expression of the normal integrin repertoire on the cell surface of TM cells can have significant effects on the pathophysiology of TM cells in glaucoma and that activation of αvβ3 integrin signaling is not simply involved in mediating adhesion. Rather, significant downstream signaling events would result from this activation subsequently causing major changes in the TM transcriptome and TM function. Hence, these studies, together with our previous studies, show that αvβ3 integrin is an important regulator of the organization of the actin cytoskeleton, phagocytosis and signaling of GTPases such as Rac in TM cells thus suggesting that alterations in cell-matrix signaling are likely to contribute to the dysregulation of trabecular meshwork function that leads to glaucoma.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cells10081923/s1, Table S1: Genes and transcripts uniquely upregulated in cells overexpressing a constitutively active avβ3 integrin (CAβ3), Table S2: Genes and transcripts uniquely downregulated in cells overexpressing a constitutively active avβ3 integrin (CAβ3), Table S3: Genes and transcripts uniquely upregulated in cells overexpressing a wild type avβ3 integrin (WTβ3), Table S4: Genes and transcripts uniquely downregulated in cells overexpressing a wild type avβ3 integrin (WTβ3), Table S5: Commonly upregulated genes and transcripts in cells overexpressing either a wild type (WTβ3) or a constitutively active (CAβ3) avβ3 integrin, Table S6: Commonly downregulated genes and transcripts in cells overexpressing either a wild type (WTβ3) or a constitutively active (CAβ3) avβ3 integrin, Table S7: Classification of Protein-Coding Transcripts Altered in WTβ3 and/or CAβ3 cells, Table S8: Classification of Non-Coding Genes and Novel Transcripts exhibiting DGE in WTβ3 and/or CAβ3 cells, Table S9: Genes differentially expressed in CAβ3 cells that are targeted by TGFβ2, BMP and/or Wnt signaling . Supplementary file: edgeRglm_GENE_3 contrasts_low-high_logFC_no NS genes.xlsx.