An examination of the regulatory mechanism of Pxdn mutation-induced eye disorders using microarray analysis

Introduction

Anterior segment dysgenesis (ASD), which is caused by abnormal development of the anterior eye tissues, is a collection of heterogeneous disorders which affect numerous tissues including the iris, cornea and lens (1). It is closely linked with an increased risk of inherited eye disorders such as congenital corneal opacity, cataracts and glaucoma (2–4), of which glaucoma is considered to be the second highest cause of blindness worldwide. It has been estimated that the number of patients with glaucoma will be 79.6 million in 2020, and this number will increased by 31.5% compared to 2010 (5). Moreover, developmental glaucoma accounts for approximately one-tenth of pediatric blindness (6). Although several advanced techniques, such as combined trabeculotomy-trabeculectomy, have been utilized to manage developmental glaucoma, the long-term outcomes are still not favorable (7). Thus, it is essential to find more effective approaches, such as biological therapeutic methods, for the treatment and prevention of glaucoma.

Of past research dedicated to uncovering the pathogenesis of ASD, one study discovered that mutations in the transcription factor (TF) genes such as paired box 6 (Pax6), forkhead box C1 (Foxc1) and forkhead box E3 (Foxe3) led to ASD (8). The TF gene peroxidasin (Pxdn) is closely linked with the extracellular matrix and is expressed in the corneal epithelium. Defective Pxdn may impair the integrity of the basement membrane by damaging the structure of its major constituent, collagen IV (9). Previous studies have verified that the Pxdn mutation is implicated in severe forms of ASD, including developmental glaucoma and congenital cataracts in human and mice (10,11). Moreover, the latter study further investigated the gene alterations and physiological changes which occurred after Pxdn mutation, by comparing the embryo tissues from the KTA048 mutant line and those from the control mouse model (11). However, the authors emphasized different stages, and the regulatory mechanisms of the Pxdn mutation in eye disorder-development, particularly in relation to developmental glaucoma progression, were not well elaborated.

Therefore, in the present study we used the microarray dataset GSE49704, which was established by Yan et al (11), to identify the differentially expressed genes (DEGs) between Pxdn mutation embryo and normal tissues. Additionally, functional enrichment analysis and protein-protein interaction (PPI) network analysis were also performed on these DEGs in order to explore their biological processes (BPs) and the potential correlations. We used the GeneCodis database to construct a transcriptional regulation (TR) network involving TFs and miRNAs. By comparing the overlapped genes in PPI and TR networks, crucial genes relating to the Pxdn mutation were screened with a sub-network extracted from the TR network. Using these bioinformatic methods, we aimed to provide a comprehensive elucidation of the regulatory mechanisms in Pxdn mutation-induced eye disorder development and provide novel biomarkers which may be used for the prognosis and prevention of these disorders.

Materials and methods

Microarray dataset

The microarray dataset with the accession number GSE49704 was downloaded from the GEO database (Gene Expression Omnibus; http://www.ncbi.nlm.nih.gov/geo). The mRNA expression profile comprised 4 mouse samples from embryos with Pxdn mutation and 4 samples of normal mouse tissues (control). The platform used was the GPL6885 Illumina MouseRef-8 v2.0 expression beadchip (Illumina, San Diego, CA, USA).

Data preprocessing and DEG identification

The raw data were subjected to preprocessing, e.g., background correction, normalization and conversion to gene symbol from the probe level, using the affy package of Bioconductor R, as previously described (12). When multiple probes corresponded to a single gene, the average level of probes was calculated as the expression value of the specific gene.

Subsequently, the limma (linear models for microarray and RNA-seq data) package of Bioconductor R (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (13) was used, as previously described, to select the DEGs between Pxdn mutation embryo tissues and the control tissues, based on the empirical Bayes test. The cut-off values for DEG identification were P<0.05 and |log2 (fold-change)| >2.

Functional enrichment analysis of the DEGs

GO (gene ontology; http://www.geneontology.org/) enrichment analysis (14) was carried out, as previously described. Briefly, enrichment analysis was undertaken for the screened DEGs in order to explore their potential roles in eye disorder development resulting from Pxdn mutation, using the DAVID (Database for Annotation, Visualization and Integration Discovery; http://david.abcc.Ncifcrf.gov/) (15) online tool as previously described. The over-represented terms with P-values <0.05 were considered to indicate significantly enriched processes of the DEGs.

PPI network establishment

To further explore the potential interactions of these DEGs on the protein level, the DEGs were mapped into STRING (Search Tool for the Retrieval of Interacting Genes/Proteins; http://string-db.org/), a database which provides comprehensive data such as high throughput and genome-wide association and the validated or predicted pairwise PPIs; we chose the mouse as species (16). The pairwise interactions with a combined score ≥0.4 were filtered out to construct the PPI network, and were visualized using Cytoscape software (http://cytoscape.org/), as previously described (17). In the network, a protein served as a 'node', and the 'degree' of a node, which was calculated by topological structure analysis of the PPI network, indicated the interaction numbers of that particular protein. The nodes with high degrees were selected, and the encoded genes were classified as group 1.

Constructing a TR network of the DEGs

Considering that understanding the TR network related to eye disorders induced by the Pxdn mutation may contribute to elucidating the molecular mechanisms of the pathogenesis of these disorders, it was deemed necessary to build a TR network of the DEGs by comparing information from the GeneCodis database (http://genecodis.cnb.csic.es), which was capable of integrating the annotated files from various bioinformatic analyses, and functional enrichment analyses using GO and the Kyoto Encyclopedia of Genes and Genomes, protein domain analysis (InterPro motifs), TF and miRNA regulation analysis. Finally, we merged the data into a single annotation file by computing a statistical rank score, as previously described (18). In the present study, we emphasized the TR analysis of DEGs involving miRNAs and TFs, and the TR network containing pairwise miRNA-DEG and TF-DEG interactions were constructed and visualized by Cytoscape software, as previously described (17).

In the present study, the node in the network was referred to as a DEG, an miRNA or a TF. By applying the aforementioned topological structure analysis to the PPI network, the degrees of the nodes in this network were also calculated. The crucial nodes with high degrees were screened and classified as group 2.

Integrating the analysis of two networks

By comparing the nodes in group 1 with those in group 2, from the two networks, the nodes which overlapped to a high degree were identified. Subsequently, the sub-network of the overlapped nodes was extracted from the TR network.

Results

DEGs between Pxdn mutation embryo tissues and the control tissues

According to the t-test method and the criteria of P<0.05 and |log2 (fold-change)| >2, a total of 121 DEGs between Pxdn mutation embryo tissues and the control tissues in mice were identified, of which 75 were upregulated and 46 were down-regulated. The heat map depicts the relative expression levels of these DEGs (Fig. 1).

Figure 1 Heat map of clustering analysis of differentially expressed genes. The x-axis represents the samples from the profile: cases 1–4 denote samples from peroxidasin (Pxdn) mutation embryo tissues of mice, and control 1–4 denote samples from the embryo tissues of normal mice; the y-axis represents genes.

Functional enrichment using the GO database

After taking into account the threshold for the significant GO terms, the over-represented processes which the DEGs may participate in were screened in the present study. As is clearly demonstrated in Table I, significant BP terms such as neuron development (GO: 0048666), regulation of cellular localization (GO: 0060341), sensory organ development (GO: 007423), polysaccharide metabolic process (GO: 0005976) as well as neuron differentiation (GO: 0030182) were significantly enriched for the DEGs.

Table I Significantly enriched biological process terms for the identified differentially expressed genes (ranked by P-values).

Table I

Significantly enriched biological process terms for the identified differentially expressed genes (ranked by P-values).

GO search term	P-value	Genes
GO: 0048666~neuron development	0.00416	Sema5a, Ret, Stmn3, Myo7a, Pak1, Gli3, Slit1
GO: 0060341~regulation of cellular localization	0.009975	Sncg, Rab3b, Celsr1, Tnfrsf4, Gli3
GO: 0007423~sensory organ development	0.010999	Cdkn1b, Myo7a, Mitf, Celsr1, Gli3, Klf4
GO: 0005976~polysaccharide metabolic process	0.012321	B3galt6, Il6st, Chi3l1, Ppp1cb
GO: 0030182~neuron differentiation	0.017787	Sema5a, Ret, Stmn3, Myo7a, Pak1, Gli3, Slit1

PPI network of the DEGs

By combining the information in the STRING database, a PPI network was established, comprising 48 nodes and 148 interactions (Fig. 2). The nodes with high degrees ≥5 were considered crucial to the network and as a result, 25 critical nodes, including troponin C type 2 (fast) (Tnnc2; degree = 5), troponin I type 2 (Tnni2; degree = 5), troponin T type 3 (Tnnt3; degree = 5), P21 protein (Cdc42/Rac)-activated kinase 1 (Pak1; degree = 5), titin (Ttn; degree = 6), parvalbumin (Pvalb; degree = 6), somatostatin (Sst; degree = 6), actin, alpha 1, skeletal muscle (Acta1; degree = 8), ras-related C3 botulinum toxin substrate 1 (Rac1; degree = 9), cell division cycle 42 (Cdc42; degree = 9), cyclin A1 (Ccna1; degree = 9), cyclin D3 (Ccnd3; degree = 9), cyclin E2 (Ccne2; degree = 10), cyclin-dependent kinase 6 (Cdk6; degree = 11), cyclin D2 (Ccnd2; degree = 11), mechanistic target of rapamycin (Mtor; degree = 12), cyclin A2 (Ccna2; degree = 12), retinoblastoma 1 (Rb1; degree = 13), cyclin-dependent kinase 4 (Cdk4; degree = 13), cyclin-dependent kinase 2 (Cdk2; degree = 14), cyclin D1 (Ccnd1; degree = 15), cyclin-dependent kinase inhibitor 1B (Cdkn1b; degree = 15), v-akt murine thymoma viral oncogene homolog 1 (Akt1; degree, 16), cyclin E1 (Ccne1; degree = 16), and cyclin-dependent kinase inhibitor 1A (Cdkn1a; degree = 19) were classified as group 1.

Figure 2 Protein-protein interaction network of the differentially expressed genes (DEGs). The large circle denotes the DEG-encoded proteins and the lines between the smaller circles indicates their interactions.

TR network of the DEGs

Using the GeneCodis database, a TR network consisting of miRNAs and TFs of the DEGs was constructed, encompassing 388 nodes (213 miRNAs, 91 DEGs and 84 TFs) and 1,523 pairwise interactions (Fig. 3). In total, 120 nodes (66 DEGs, 38 miRNAs and 16 TFs) with high degrees ≥7 were classified as group 2 nodes.

Figure 3 Transcriptional regulation network of the differentially expressed genes (DEGs). The red circles represents DEGs, and pink squares denote transcription factors (TFs), while green triangles represent miRNAs. The arrows indicate the targeting interactions of DEG-miRNA and DEG-TF.

By comparing the genes in groups 1 and 2, seven overlapped genes with high degrees were identified, namely Cdkn1b, Acta1, Pak1, Pvalb, Sst, Tnnt3 and Ttn. The degrees of the seven genes in the two networks are indicated in Table II. Subsequently, a sub-network involving these seven genes was extracted from the TR network, comprising 156 nodes (7 DEGs, 99 miRNA and 50 TFs) and 205 interactions (Fig. 4). Notably, Cdkn1b was predicted as the target of miRNAs such as mmu-miR-24, mmu-miR-28*, mmu-miR-34a, mmu-miR-369-5p and mmu-miR-449b and TFs including forkhead box O4 (FOXO4) and activating enhancer binding protein 4 (AP4).

Figure 4 Sub-network of the overlapped differentially expressed genes (DEGs) extracted from the transcriptional regulation network. The red circles represent DEGs, and pink squares denote transcription factors (TFs), while triangles represent miRNAs. The arrows indicate the targeting interactions of DEG-miRNA and DEG-TF.

Table II Degrees of overlapped genes in both networks.

Table II

Degrees of overlapped genes in both networks.

Genes	Degree-TRN	Degree-PPI
Acta1	33	8
Cdkn1b	27	15
Pak1	26	5
Pvalb	27	6
Sst	18	6
Tnnt3	25	5
Ttn	49	6

[i] Degree-TRN, degree in transcriptional network; degree-PPI, degree in protein-protein interaction network.

Discussion

Pxdn plays an important role in basement membrane synthesis, and its mutation is linked to ASD, which can easily develop into developmental glaucoma and congenital cataracts (19). In the present study, we selected a total of 121 DEGs (75 upregulated and 46 downregulated) between the Pxdn mutation embryo tissues and the control tissues, which were enriched in BPs such as neuron development, polysaccharide metabolic process and neuron differentiation. Seven vital overlapped genes of the PPI and TR network were identified, including Cdkn1b, Acta1 and Tnnt3, with the sub-network extracted from the TR network. Of these, Cdkn1b was predicted be the target of miRNAs such as mmu-miR-24, and TFs including FOXO4 and AP4.

CDK2 and CDK4 are two cyclin-dependent kinases (CDKs) that positively regulate cell cycle activities (20). In a previous study, it was noted that the cyclin-dependent kinase inhibitor 1B (Cdkn1b, also known as p27kip1) encoded protein binds to the cyclin E-Cdk2 or cyclin D-Cdk4 complexes to prevent the activation of the complex, and thus serves as an inhibitor of cell cycle progression at the G1 stage. The overexpression of Cdkn1b in Tenon's capsule fibroblasts after glaucoma filtration surgery resulted in the downregulation of Cdk2 and Cdk4 and the inhibition of fibroblast proliferation, which contributed to decrease the severity of the surgical outcome (21). In addition, although Cdkn1b was not the only regulator of Sertoli cell proliferation in adult mouse testes, its downregulation was linked to the elevated proliferative activity of Sertoli cells, suggesting that the protein played a role as a suppressor in Sertoli cell proliferation and regulation (22). These previous studies verified that Cdkn1b played an important role in inhibiting proliferation during cell cycle progression. Notably, fibrovascular proliferation into the anterior segment was considered a rare complication of acute angle-closure glaucoma (23). Considering that the Pxdn mutation may also influence the proliferation and differentiation during eye development (11), and Cdkn1b was among the downregulated DEGs in the Pxdn mutation line in our present study, it may be inferred that Cdkn1b serves as an inhibitor of proliferation during eye development and that the Pxdn mutation-induced Cdkn1b decrease is related to eye disorders such as glaucoma.

Extensive studies have been conducted which report the transcriptional regulation of Cdkn1b by miRNAs and TFs. miR-452 has been reported to dramatically accelerate proliferation by targeting the 3′-untranslated region (UTR) of Cdkn1b, which eventually contributed to tumorigenesis in hepatocellular carcinoma (24). miR-222 may also bond to the 3′ UTR of Cdkn1b, and it has been suggested that it is an indicator of liver fibrosis (25). Moreover, in human prostate carcinoma cell lines, it was discovered that the expression of miR-221 and miR-222 affected cell proliferation by targeting Cdkn1b (26). As shown in Fig. 4, Cdkn1b was predicted to be targeted by cohorts of miRNAs such as mmu-miR-24, mmu-miR-28*, mmu-miR-34a, mmu-miR-369-5p and mmu-miR-449b. It has previously been noted that miR-24 inhibits cell cycle progression and directly suppresses the transcription of v-myc avian myelocytomatosis viral oncogene homolog (MYC), an important TF in cell cycle progression which regulates the transcription of Cdkn1b (27). Although no further direct evidence exists to illustrate the regulatory relationship between Cdkn1b and miR-24, this study may provide a clue for the potential targeting relationship between these molecules in eye disorders development. However, further experimental validation is warranted.

With regard to the TFs that targeted Cdkn1b, it has been reported that the binding of the TF PAX6 to multiple regulators which are responsible for cell cycle progression and proliferation, including Cdkn1b, may account for the regulation of cell proliferation and determination in the cortex (28). Moreover, it has been noted that Cdkn1b is directly regulated by the TF E2F1 and suppresses baculoviral IAP repeat containing 5 (Birc5) expression via the inhibition of CDK activity, which results in cell cycle arrest. The induction of Cdkn1b by E2F1 enabled E2F1 to inhibit the expression of Birc5 in neuroblastomas (29). In the present study, Cdkn1b was linked to TFs including FOXO4 and AP4. FOXO4 is a TF that is involved in the regulation of the insulin signaling pathway (30). The peptidylprolyl cis/trans isomerase, NIMA-interacting 1 (Pin1) encoded protein acts as a catalyst in the regulation of post-phosphorylation conformation of its substrates by binding to the phosphorylated ser/thr-pro motifs (31). In relation to tumor development, it has been indicated previously that Pin1-induced deubiquitylation results in the inhibition of the nuclear translocation of FOXO4 and the suppression of FOXO4 targeting Cdkn1b (32). The TF AP4 has been noted to suppress the transcription of Cdkn1a, another CDK inhibitor, by occupying four CAGCTG motifs in its promoter. By binding to specific motifs of CACGTG in the first intron of AP4, c-MYC directly regulates AP4 expression and thus contributes to maintaining the proliferative stage in breast cancer cells (33). This evidence substantiates our hypothesis that Cdkn1b is the target of FOXO4 and AP4 and that transcriptional regulation by the TFs plays a significant role in eye disorder development caused by the Pxdn mutation.

Of the remaining six vital genes identified in the present study, Acta1, which encodes alpha actin family protein, has previously been listed in the DEGs between retinal ganglion cells (RGCs) and glial cells (34). Given that the loss of the RGCs through apoptosis contributes to the progression of glaucoma (35), we speculate that Acta1 is also linked to eye disorder development, particularly glaucoma progression. Tnnt3 has been identified as a DEG between human trabecular meshwork in patients with primary open-angle glaucoma and normal subjects (36), and in this study it was also suggested that it is involved in eye disorder regulation. Although there is not yet sufficient evidence to fully identify correlations between these genes and eye disorder development, or the regulatory relationships between these genes and TFs and miRNAs, all the predictions described in the present study provide a novel insight into the transcriptional regulators in eye disorders such as glaucoma, which are induced by the Pxdn mutation, and this link should be borne in mind for follow-up studies.

In conclusion, using the GSE49704 dataset, in the Pxdn mutation embryo tissues of the mice, we identified several crucial genes, including Cdkn1b, Acta1 and Tnnt3, which were related to the progression of eye disorders, such as glaucoma. Notably, it may be the case that Cdkn1b exerts its roles by inhibiting proliferation and is mediated by mmu-miR-24 and targeted by the TFs of FOXO4 and AP4, during the regulation of eye disorder development. However, further experimental validationis warranted.

Acknowledgments

The present study was supported by the Hubei Provincial Natural Science Foundation of China (no. 2014CFB366).

References