Transcriptional regulatory network analysis during epithelial-mesenchymal transformation of retinal pigment epithelium

Computational analysis of gene regulatory regions

Identification and analysis of the regulatory regions of genes expressed in RPE cells was performed using the Promoter Analysis and Interaction Network Toolset version 3.3 (PAINT) program [12]. The target gene set analyzed comprised sixty genes differentially expressed during EMT of RPE cells, selected from previous literature reports and additional genes under study in our laboratory [10,11] (Table 1). This set was divided into two subsets representing those preferentially expressed in the undifferentiated or differentiated state, and further refined by functional gene ontology [13]. Homologous genes in humans, mice, rats, and chickens were identified using the Ensembl database and the Ensembl gene ID for each was used in the PAINT program [14,15]. For each gene, PAINT identified the putative transcription start site (TSS) and subsequent promoter analysis to identify transcription response elements (TREs) was performed to 5,000 bp upstream of the TSS, with exclusion of complimentary strand analysis and 1.0 core similarity threshold [12]. In addition, sequence data corresponding to the first intron of each gene was retrieved from Ensembl and also entered into PAINT, as FASTA formatted sequences, for further analysis. TREs within both the promoters as well as first introns were identified using the MATCH/TRANSFAC database [16]. These TREs were entered into the Feasibility Network Builder module of PAINT (FeasNet Builder), which constructed a candidate interaction matrix (CIM), a graphic representation of the occurrence of these TREs within the gene set. Enrichment analysis was performed using PAINT to compute the Fisher's exact test p-values indicating relative over-representation of TREs within the selected gene set as compared to the larger background gene set, the 588 genes present on the Human Atlas Array (Cat. No. 7740–1; Clontech, Palo Alto, CA) used to originally identify the set of sixty differentially regulated RPE genes [10,11]. In each of the analyses, multiple testing correction was applied using a false discovery rate (FDR) estimate [17]. In all cases, the multiple testing corrections were not informative as they did not result in any over-represented TREs, since the FDR was consistently above 90% for all TREs. Therefore, we followed a discovery approach and chose a threshold of p<0.1 on the Fisher's exact test p-value to identify those TREs to be included in further filtering as described below. Models of RPE gene interaction networks, based upon results of the FeaseNet analysis, were graphically generated using GraphViz [18].

Generation of models for gene regulatory regions and selection of targets for biochemical analysis

Global regulatory models for gene sets coordinately expressed in RPE cells, as well as for individual genes differentially expressed in differentiated or undifferentiated RPE cells, were developed based upon PAINT-derived computational data and were constructed by comparing phylogenetically conserved transcriptional regulation across the four species: human, mouse, rat, and chicken. To establish criteria for selection of specific transcription factors for further analysis, we assigned values to each TRE based upon their frequency of detection across the coordinately expressed gene sets and their evolutionary conservation. An evolutionary conservation factor (ECF) of 1 to 3 points was assigned to each TF, where 1 indicates presence on human genes, 2 indicates human and mouse or rat genes, and 3 indicates a presence on human and chicken as well as either mouse or rat genes. In addition, TREs were scored according to a frequency ratio (FR) derived from the ratio of the percent occurrence of a given TRE in a specific gene subset divided by its frequency of occurrence in the background gene set, for the human genome data. Ultimate TRE selection was based upon criteria filters of a combined ECF score of 2 or greater along with a FR greater than 3.

Cell culture

Primary cultures of chick embryo RPE cells (cRPE) were established from RPE tissues obtained from fertile White Leghorn chicken eggs maintained in a humidified atmosphere at 37 °C until embryonic day 10, corresponding to Hamburger and Hamilton stage 36 [19,20]. Eyes were enucleated, the anterior segment and vitreous were removed, and the posterior eyecup divided in half and incubated in HBSG (HEPES buffered saline with glucose containing 1 mg/ml glucose, 10 mM HEPES, pH 7.4, 3 mM KCl, and 0.15 mM NaCl). The neural retina was removed and the eyecup was incubated in HBSG containing 20 mM EDTA (ethylene diamine tetraacetic acid) for 1 h. The eye cups were then rinsed for an additional 30 min in HBSG and the RPE was dissected from the choroid. RPE tissue was collected by gentle centrifugation and resuspended in MEM (Eagle's minimum essential medium; Sigma-Aldrich, St. Louis, MO) supplemented with 10% FBS (Fetal Bovine Serum; Sigma-Aldrich), 0.22% sodium bicarbonate and 1% penicillin, streptomycin and amphotericin (Gibco/Invitrogen, Carlsbad, CA). The RPE tissue was mechanically dissociated into single cells and plated onto plastic six well tissue culture plates at a density of 1 eye/well. The plates were previously coated with 1.1 μg/cm² mouse laminin (Sigma-Aldrich) 12 h before cell plating. The cultures were maintained in a humidified atmosphere of 5% CO₂ at 37 °C for 48 h at which time all non-adherent cells were removed by rinsing with fresh medium [21].

The human RPE-derived cell line ARPE-19 [22] was maintained cultured in DMEM-F12 (Catalog number D8900; Dulbecco's modified eagles minimum nutrient mixture F12 Ham; Sigma-Aldrich) supplemented with 5% FBS, 2 mM L-glutamine, 0.348% sodium bicarbonate and 1% Antibiotic-Antimycotic in T25 culture flasks in a humidified atmosphere of 5% CO₂. To maintain cells in an undifferentiated state, they were passaged before obtaining confluence. To obtain differentiated cells, cells were grown to confluence and then maintained in DMEM-F12 as above except that the serum was reduced to 1% [23]. These cultures reach confluence 2–3 weeks after passaging and differentiate within 4–6 weeks, though the cultures can be kept in a differentiated state for extended culture periods. After 4–6 weeks in culture the cells exhibit hexagonal packing of pigmented, polarized epithelia and expression of CRALBP and RPE65 typical of morphological and biochemical markers, respectively, of RPE cells in vivo [22].

RNA extraction and reverse transcriptase polymerase chain reaction amplification

Total RNA was extracted from all cell types using the Micro-to-Midi Total RNA Purification System (Invitrogen, Carlsbad, CA) following manufacturers specifications for isolation of RNA from animal tissues. Purity and concentration of RNA from each sample was assessed on a spectrophotometer at 260 nm wavelength (A₂₆₀). All primers used in this study are listed in Appendix 1 and were designed using the web-based tool GeneFisher Primer Design program, and all ranged between 20 and 22 nucleotides in length, with melting temperatures of 50–65 °C and G-C content between 40%–60% [24]. Primer specificity was determined by the nucleotide-nucleotide basic local alignment search tool (BLASTn) set to the specific species genome, with acceptance of a primer pair based upon expect-value of less than 1 (e-value <1) [25]. RT–PCR was performed using the SuperScript III One-Step RT–PCR System with Platinum Taq (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. The reaction parameters were (a) cDNA synthesis, 1 cycle at 55 °C for 30 min, (b) denaturation, 1 cycle at 94 °C for 2 min, (c) amplification, 40 cycles at 94 °C for 15 s, (d) melting at primer-specific temperatures for 30 s, (e) extension, 68 °C for 1 min, and (f) final extension, 1 cycle at 68 °C for 5 min. RT–PCR products were analyzed by agarose gel electrophoresis on 1.25% gels containing ethidium bromide (1 μg/ml). The resultant bands were visualized and recorded under ultraviolet light using the Kodak 1D photo system.

Transcriptional regulatory network analysis of differentially expressed retinal pigment epithelium genes reveals gene set-specific clusters of transcription response elements

To identify target genes for the analysis of transcription response elements by PAINT, a gene set containing two clusters representing those genes whose expression is preferentially associated with the undifferentiated versus differentiated state of the ARPE-19 cells was established (Table 1). Thus each gene cluster would be predicted to be preferentially associated with one or more TREs involved in their respective coordinate regulation of expression. The genes within the differentiated and undifferentiated clusters were separately analyzed using PAINT to identify and statistically analyze the occurrence of TREs within the promoter region, including 5,000 bp of sequence upstream from the transcription start site, as well as the entire first intron sequence. The complete results of this analysis by the FeasNet Builder module of PAINT are listed in Appendix 2 for the promoter regions and Appendix 3 for the first introns. The analysis identified those TREs that are statistically over-represented on the promoters and first introns of genes within either of the two gene clusters, as compared to their occurrence with the larger background set of promoter and intron regions. TREs with a p-value of less than 0.1 were deemed significant and included for further analysis. These data were then converted into graphic representations of the TRE network or CIM, using the FeasNet Viewer module of PAINT, as illustrated in Figure 1 for the promoter regions and in Figure 2 for the first introns. Inspection of the matrix patterns for the respective differentiated and undifferentiated gene clusters indicated that distinct subsets of TREs were associated with each of these two clusters, suggesting that these TREs could represent components of the networks regulating coordinated expression of genes within each cluster. Figure 1 and Figure 2 represent a subset of the complete CIM, and include only those TREs with a p-value of less than 0.1. The complete CIMs representing the results of the analysis of the full 60 gene set are tabulated in Appendix 2 and Appendix 3, as well as illustrated in Appendix 4 and Appendix 5. The results of the CIM analysis indicate that several specific TREs are statistically overrepresented within the regulatory regions of genes within each gene subset, and furthermore that several of these TREs differ between the two subsets. These distinctly represented TREs are thus candidates for further analysis regarding their potential regulatory role in coordinating gene expression during RPE cell differentiation.

To further visualize the potential interrelationships of the candidate TREs and genes within each cluster, the data generated by the FeasNet Viewer module was graphically presented using GraphViz to generate a regulatory network diagram for the human gene set, as illustrated in Figure 3 for the promoter regions and Figure 4 for the first introns. These results identify nodes within each network and further identify those TREs with the potential to coordinately regulate subsets of genes within each cluster. The TREs that comprise the nodes of this visualization represent all those TREs which are overrepresented among the genes of each cluster, and were found to be associated with as few as one gene within each cluster to as many as 28 genes within a cluster. Those TREs exhibiting high levels of interconnectivity are among the best candidates for coordinate regulation of RPE genes.

To identify those TREs of particular significance in the differential regulation of RPE genes during EMT, we next performed an analysis comparing all the overrepresented TREs from the CIM analysis with respect to their relative frequency of occurrence between the differentiated, undifferentiated and background gene clusters. While the previous analyses identified TREs overrepresented in either or both of the two gene clusters, the frequency analysis further distinguished those TREs preferentially associated with either the differentiated or undifferentiated gene cluster. The results, which are shown in Figure 5 for the promoter regions and Figure 6 for the intron regions, indicate that a select subgroup of the TREs can be assigned as potential regulators of one of each of the two cell states. For example, in the promoter regions, Hand1:E47 and COMP1 show increased relative frequency among down-regulated genes, whereas Oct-1 and SREBP-1 do so among upregulated genes. In the first intron regions, CDP and IRF show increased relative frequency among down-regulated genes, whereas GATA-3 and MAZR do so among upregulated genes.

Phylogenetic comparison of PAINT analyses identifies evolutionarily conserved transcriptional regulatory elements across divergent species

The above studies of gene regulatory networks regulating RPE cell differentiation were focused on human genome sequence analysis. However, evolutionary conservation analysis of regulatory regions can also be useful for identifying functionally important sites. Therefore, to further identify those elements of these hypothetical networks that would most likely be of functional significance due to their evolutionary conservation, we performed parallel PAINT analyses similar to those described above for human genomic elements for three additional model species: mouse (Mus musculus), rat (Rattus norvegicus), and chicken (Gallus gallus). The results of the FeasNet Builder analysis in the three additional species gene sets indicate that, as observed for the analysis of human genes, several overrepresented TREs are again identified (Appendix 6, Appendix 8, and Appendix 10 for the promoter regions and in Appendix 7, Appendix 9, and Appendix 11 for the first introns) for each species. The graphical candidate interaction matrix for each of these three species is shown in Appendix 12, Appendix 13, Appendix 14, Appendix 15, Appendix 16, and Appendix 17, while the GraphViz output of the FeasNet Viewer analyses illustrating the regulatory network diagrams for these three species is shown in Appendix 18, Appendix 19, Appendix 20, Appendix 21, Appendix 22, and Appendix 23, which both again indicate the presence of a complex regulatory network within each of these species. As for the human genome, we also compared the frequency of occurrence of each identified TRE between each gene cluster and the background gene set for each of the three model species. These frequency analyses are shown in Appendix 24, Appendix 25, Appendix 26, Appendix 27, Appendix 28, and Appendix 29. As for the human genome analysis discussed above, these results of the CIM, GraphViz and frequency analyses identify TREs for each of the three species which are also candidates for regulation of gene expression during RPE differentiation. Comparison of the results from the human as well as model system analyses identified both phylogenetically conserved as well as species-specific TREs, which were incorporated into criteria tables and gene regulatory region models as described below (Table 2 and Figure 7).

Compilation of human and cross-species transcriptional regulatory network analyses generates global as well as gene-specific regulatory models

To generate a global working model of gene expression regulation in the RPE, we compared the above results across the four species analyzed and identified those TREs of highest evolutionary conservation, as well as highest frequency of occurrence, by calculation of an evolutionary conservation factor (ECF) score, as well as a frequency ratio (FR) score, for each TRE identified in the computational analysis, as described in Methods. The results of these calculations are shown in Table 2, with the table containing only those results for TREs that passed the three criteria of a p-value >0.1, an ECF greater or equal to 2, and an FR greater or equal to 3. The full data set of calculated ECFs and FRs for all the TREs identified in the computational analysis is shown in Appendix 30, for the promoter regions and Appendix 31 for the first intron regions. These calculations provide a scheme for ranking TREs, and by extension their corresponding TFs, for inclusion into models of RPE gene regulation, as well as for further studies probing their expression and function in RPE differentiation.

The TREs passing all these criteria, and hence representing potentially evolutionarily conserved regulators of RPE gene expression, were incorporated into archetypal gene models for the coordinated up- and down-regulation of genes during RPE cell differentiation, as shown in Figure 7. Although evolutionary conservation is a potential indicator of conserved function, in recognition of the likelihood that species-specific aspects to coordinated gene regulation also exist, we also developed an additional model with modified ECF criteria focused on regulation in the human genome, such that incorporation of TREs into this second model was dependent on the TRE occurring within the human gene set and any one additional species, while maintaining the same p-value and FR criteria filters. The resulting models are shown in Figure 8, and share many elements of the trans-species model, in some cases eliminating TREs not found in the human analysis, while adding some others due to the less stringent ECF criteria. Finally, additional gene-specific regulation models were constructed for selected genes representing paired members of multigene families, or other markers, with reciprocal expression during RPE cell differentiation. These paired models were constructed for the cell adhesion proteins N- and R-cadherin, for the lactate transporters MCT3 and −4, and for α-smooth muscle actin and RPE65 (Figure 9). These gene regulation models identify those TREs most likely to coordinate expression of specific genes as well as broader sets of up- and down-regulated genes during RPE cell differentiation, and provided targets for validation of these models through direct biochemical analysis as described below [22,23,26].

Validation of gene models by identification of transcription factors whose expression is dependent upon the state of retinal pigment epithelium cell differentiation

The above gene regulatory models identified TREs that could serve as elements of the gene regulatory network during RPE cell differentiation. For these TREs to play a role in regulation of their associated genes, the activity of the corresponding transcription factors would be expected to be dynamically regulated at appropriate times to effect such control. One common mechanism of regulation for TF activity is at the transcriptional level resulting in differential steady-state mRNA expression levels. Thus, to test the computationally derived models, we used RT–PCR to assay for the presence of mRNA encoding transcription factors (TFs) predicted to play a regulatory role in the RPE. The results obtained from comparison of mRNA extracted from undifferentiated and differentiated RPE cells, using both the human ARPE-19 cell line and primary embryonic chicken RPE cells, are shown in Figure 10, Figure 11, Figure 12. The analyses of ARPE-19 cells were performed using both undifferentiated cells that been in culture for one month, exhibiting morphology of fusiform, unpigmented mesenchymal cells, as well as differentiated ARPE-19 cells exhibiting morphology of polygonal, darkly pigmented, epithelial cells (Figure 10A,B). RT–PCR was first performed on each respective cell population targeting mRNAs for α-SMA and RPE65, respective markers of the mesenchymal undifferentiated and epithelial differentiated state of RPE cells [22]. These results indeed demonstrated that mRNA encoding α-SMA, but not RPE65, was expressed in undifferentiated cells, whereas RPE65 mRNA was readily detected among differentiated cells with a reduced level of α-SMA mRNA (Figure 10C,D). Semi-quantitative RT–PCR was used to further distinguish these levels of α-SMA mRNA, which more clearly distinguished these two cell states (Figure 10E). Having verified that these two markers were distinctly expressed between these two cell populations, we then used RT–PCR to determine the levels of mRNA encoding the specific transcription factors previously identified by the computational analysis, which were expected to fall into three categories, exhibiting either quantitative, qualitative, or no differences between the two test cell populations. Of these TFs, mRNAs encoding four were found to be reciprocally expressed in differentiated versus undifferentiated ARPE-19 cells, with Oct-1 and TFE3 detected only in differentiated cells, and Core Binding Factor and NKX3A detected only in undifferentiated cells (Figure 11). mRNAs encoding additional transcription factors, including GATA-1, IRF-1, and SMAD3, were detected in both cell states (Figure 11A). Semi-quantitative RT–PCR was then used to further analyze differences in the expression patterns of these factors, with GATA-1 detected at higher levels in differentiated cells, whereas IRF-1 and SMAD3 were detected at higher levels in undifferentiated cells (Figure 11B).

Similar analyses were performed using freshly isolated chicken RPE tissues and primary chick RPE cultures. When cultured, chick RPE cells re-enter the cell cycle and de-differentiate [27]. RNA was prepared from freshly isolated RPE cells as well as from cells cultured for five days in vitro, and both cell populations were probed for TFs corresponding to those TREs identified in association with both the global as well as gene-specific regulatory models. The results indicated that, similar to ARPE19 cells, and consistent with model predictions, the primary chick RPE cells also revealed reciprocal TF expression between differentiated and undifferentiated cells. mRNAs encoding AML-1 and HNF-1 were detected only in the differentiated chicken RPE cells, whereas mRNAs encoding HNF-3 and SREBP-1 were detected only in the undifferentiated cells, while two additional mRNAs encoding TFs DP-1 and TFII-I were detected in both cell populations (Figure 12). Finally, as an adjunct to the PAINT-derived analyses, an additional series of RT–PCR reactions were performed to determine whether other TFs, not identified through PAINT, but known from prior studies to be involved in EMT of epithelial cells other than RPE, were expressed in chick or ARPE-19 cells. As shown in Figure 13, RT–PCR amplification of mRNA from undifferentiated and differentiated primary chick RPE cells generated similar levels of amplicons for Slug, Twist and SIP1, whereas Snail was detected at higher levels in undifferentiated cells, and LEF1 was detected only in differentiated cells. When similar analyses were performed with total RNA isolated from differentiated or undifferentiated ARPE-19 cells, Slug, Snail, Twist, and SIP1 were not detected, while SMAD2 was detected at equal levels in both samples. Interestingly, while LEF-1 was also detected in both cell populations, a distinct additional amplicon was detected in differentiated RPE cells, indicating that differential splicing of this gene transcript occurs during the course of RPE differentiation. Our inability to detect expression of certain classical mediators of EMT such as Snail, Slug, Twist or SIP1 in ARPE-19 cells, while we were able to detect them in primary cultures of embryonic chick RPE cells, may be related to the different stages of development represented by these two model systems (embryonic versus adult), to the unique properties of RPE cells compared to other epithelial cell types that may exhibit species-specific differences, or to some specific phenotypic property of ARPE-19 cells that arose during their derivation [28,29]. Overall, the results of the RT–PCR analyses indicate that the computational biology approach was successful at identifying transcription factors whose expression is regulated during RPE cell differentiation, and thus may play a role in control of differential gene expression and modulation of RPE cell phenotype.