Comparison of genome-wide gene expression in suture- and alkali burn-induced murine corneal neovascularization

Animals

Balb/c mice, 6–8 weeks old, were used in this research. All mice were purchased from Beijing Pharmacology Institute, Chinese Academy of Medical Sciences (Beijing, China). Use of animals was approved by institution and observed the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research.

Corneal neovascularization models

Mice were anesthetized with ketamine (50 mg/kg) and chlorpromazine hydrochloride (10 mg/kg) by intraperitoneal injection. Compound tropicamide eye drops (Santen, Osaka, Japan) and 0.5% proparacaine hydrochloride (Alcon-Couvreur, Puurs, Belgium) were applied topically for corneal anesthesia. CorNV was induced by suture placement or chemical burn, as previously described by others [9]. Briefly, a corneal trephine 2 mm in diameter was pressed lightly on the central cornea to make a circular mark. Three interrupted 10–0 polypropylene sutures (MANI Inc., Togichi, Japan) spanning the mark were placed through the epithelial and stromal layers, but without penetrating the endothelial layer. To cause alkali burn to the cornea, a piece of filter paper (2.0 mm in diameter) soaked with 2 μl of 1 mol/l NaOH solution was placed on the central corneal surface for 40 s, followed by immediate rinsing with 30 ml of 0.9% saline buffer.

Histology

Mice eyeballs were formalin fixed, paraffin embedded, and sectioned at a thickness of 4 μm for routine histological processing. After staining with hematoxylin and eosin, vessels and cell infiltration were examined by light microscopy.

Isolation of total RNA and microarray procedure

The procedure for total RNA isolation and microarray assay was described earlier [10]. In brief, at the desired time points after CorNV induction, corneas were excised using a 2 mm diameter trephine and placed in ice-cold TRIzol reagent (Invitrogen, Gaithersburg, MD), with five corneas from each model pooled into one sample. The untreated corneas from the same mice were used as controls. Three pairs of samples were prepared for each model. Total RNA was extracted using isopropanol precipitation and was purified using NucleoSpin RNA clean-up columns (Macherey-Nagel, Düren, Germany). Dual cRNA labeling and microarray hybridizations were performed by Capital Bio Corporation using Capital Bio cRNA labeling kits and the Capital Bio 36 K Mouse Genome Oligo Array (Capital Bio, Beijing, China) [11,12]. The array comprises 35,852 70-mer oligonucleotide probes representing approximately 25,000 genes of Mouse Genome Version 4.0 (Operon Biotechnologies, Huntsville, AL). Three replicate arrays were used for each model. After hybridization, the arrays were scanned using a LuxScan 10KA (Capital Bio), and signals were processed with LuxScan 3.0 software (Capital Bio). Intra-array normalization was done using the LOcally Weighted Scatter plot Smoothing (LOWESS) linearization method. Inter-array normalization of the whole data set was performed according to the global means of Cy5 and Cy3 signals [12].

Microarray data analysis

In each array, the probes that passed various quality check and gave signal intensity of over 1,500 were labeled as Expressed, otherwise a probe was considered Marginal (intensity 800–1500) or Absent (less than 800). Only those probes that were expressed in at least two out of three chips were taken into account for further analysis. After log² transformation of the fold values, one class t-test analysis was performed by Significance Analysis of Microarray software SAM 3.0 (Stanford University, Stanford, CA) [13] to determine their change or stability in the CorNV group compared with the control group. Those genes that gave less than a 1% false discovery rate (FDR) and no less than twofold changes were considered to be differentially expressed genes. Last, the differentially expressed genes identified in each CorNV model were compared and annotated using the Database for Annotation, Visualization and Integrated Discovery (DAVID, v6.7) with the whole murine genome as the background [14]. Gene Ontology (GO) categorization was performed using a modified Fisher exact test, and the p value for each GO category was calculated as a expression analysis systematic explorer (EASE) score [15]. An EASE score of not more than 0.01 indicated a significant enrichment. Hierarchical clusters were performed for different sets of interested genes using the Cluster 3.0 program with the Pearson correlation (uncentered) distance, average linkage [16,17], and the resulting CDT files transferred into heat maps using Java Treeview [18]. The complete sets of normalized data of this microarray assay are deposited in the NCBI Gene Expression Omnibus (GEO) with a GEO accession number GSE23347.

Antibodies and western blotting

Corneas were harvested by cutting along the centric side of limbal line and placed in RIPA lysis buffer (Beyotime, Shanghai, China) for total protein preparation. For each sample, two corneas were combined and 100 μl buffer were used. The tissues was cut with scissors into small pieces and homogenized using a tissue tearor (Biospec Products, Inc., Bartlesville, OK). After spinning at 11,500× g for 10 min, 80 μl cleared lysate were mixed with 20 μl 5× loading buffer. After boiling at 95 °C for 5 min, 10 μl samples were resolved on 12% SDS–PAGE gel and then transferred to a polyvinylidene difluoride (PVDF) membrane (Millipore, Billerica, MA). The blots were blocked in 5% non-fat dry milk dissolved in TBST (20 mM Tris PH7.5, 0.5 mM NaCl, 0.05% Tween-20) for 1 h and then incubated with the primary antibody in TBST for 2 h, followed by incubation with HRP-conjugated secondary antibody for 1 h. All incubations were done at room temperature and three washes with 10 ml TBST were applied between each step. Primary antibodies include anti-mouse PEDF antibody (sc-25594; Santa Cruz Biotechnology Inc., Santa Cruz, CA), anti-VEGF antibody (ab-3109; Abcam Biotechnology, Cambridge, MA). Secondary antibodies include Peroxidase-Conjugated AffiniPure Goat Anti-Rabbit IgG (H⁺L; ZB-2301; Zhongshan Golden Bridge, Beijing, China) and Peroxidase-Conjugated AffiniPure Anti-Mouse IgG (ZB-2305, Zhongshan Golden Bridge).The blot was developed with SuperSignal West Pico (NCI5079,Thermo Fisher Scientific, Rockford, USA) and exposed to X-ray film (Kodak, Rochester, NY). For detecting GAPDH, the PVDF membrane was regenerated using Stripping Solution (Applygen Technologies Inc., Beijing, China) and detected using an GAPDH detection kit (KC-5J5,KangChen Biotechnology, Shanghai, China) as suggested by the manufacturer.

Corneal neovascularization in two models

Both S-CorNV and CB-CorNV have been used successfully, including in our laboratory [10,19], when studying either mechanisms of CorNV development or interference from exogenous factors in the development of CorNV. The conclusions drawn using one model are presumed to be applicable to all CorNV. We previously used both models without any specific preference; by contrast, this study tried to define any commonalities or differences between these two models. Clinically, suture placement did not cause other significant changes to corneal transparency except for the development of neovascularization, but chemical burn caused total loss of transparency of the burned area, starting from the burn (Figure 1A). Later, the haze extended from the edge of the burn to the limbus. The difference in the gross appearance of corneas in these two models correlated with the histological changes (Figure 1B). While the suture did not affect the intactness of the epithelial layer in areas other than the suture punch, the chemical burn depleted the entire epithelial layer of the central corneas. Correspondingly, the suture caused inflammatory infiltration limited to the area between the suture and limbus, but the chemical burn caused heavy infiltration and edema of the whole cornea. In spite of the significant gross difference in these two models, the appearance and progress of CorNV in chemically burned corneas was slightly slower than that in sutured corneas. S-CorNV progressed most quickly at day 5, and CB-CorNV at day 6 (Figure 1), reaching maximal length at day 10 and day 14, respectively (data not shown). This paper focused on data obtained at the fast-growth time points, namely day 5 for the suture and day 6 for the chemical burn, respectively.

Comparison of differentially expressed genes in two corneal neovascularization models

Merely by checking the gross presentation or histology of corneas in these two conditions as described above, we expected CB-CorNV to affect the expression of more genes or to affect them to a greater extent than S-CorNV would. However, gene profiling using microarray showed that S-CorNV affected gene expression more than did CB-CorNV. Among all the 35,872 probes (excluding various controls) in this chip, 7,138 that passed the filter criteria were identified as being expressed in S-CorNV, and 7,109 in CB-CorNV, for a total of 7,766 probes expressed in at least one model. Namely, about 92% of all expressed probes were detected in both models, reflecting the consistency and reliability of the data sets. Statistical analysis using SAM software, with the thresholds set at FDR<1% and fold change≥2, revealed that 1,055 probes (accounting for 14.78% of all expressed genes in the S-CorNV model) were differently expressed. Among them, 586 probes were upregulated and 469 probes were downregulated. Similarly, in the CB-CorNV model, 472 probes were upregulated and 389 were downregulated, making a total of 861 probes (12.11%) differentially expressed. Among these differentially expressed genes, 530 probes in total overlapped in the two models, including 286 upregulated and 244 downregulated probes (Table 1). Some probes remained unchanged in one model but up- or downregulated in the other model. No probes manifested contradictory changes in the two models; namely, no probes were upregulated in one model while being downregulated in the other model. To allow a better overview of the changes of all 1,386 changed probes, hierarchical clustering analysis was performed, and is shown in Figure 2.

Complimentary validation of expression patterns of VEGFA and PEDF

It was interesting and unexpected to note that such classical angiogenesis modulators as VEGFA and PEDF were not among the aforementioned 1,055 genes that produced significant alterations at mRNA levels in either CorNV model. Studying the raw array data showed that VEGFA signals fell below the threshold set by analysis criteria thus regarded Absent (namely, unexpressed) in the corneas. For example, four probes with Oligo ID being M400000704, M400000705, M400000706, and M400000707, respectively, are included in the array for the VEGFA gene. The signals for these four probes in the S-CorNV group shown as “experimental sample signal/control sample signal” were 272.3±30.7/56±52.5 (mean±standard deviation for three arrays), 412.3±80.0/220.6±74.7, 159.6±43.5/31.3±161.8, and 238.3±79.1/97.6±88.8, respectively. If the above filtering threshold would be arbitrarily disregarded and the fold changes provided by the automatic analysis accompanied with the array signal results, these four probes would produce a ratio of 2.7 for VEGFA, namely VEGFA mRNA would be regarded upregulated by 2.7 fold in S-CorNV compared with control corneas. On the contrary, PEDF gave high enough signals (between 4,000 and 6,300, not shown) and were classified as “expressed” in all samples, the change in S-CorNV was not significant (the experiment/control ratios in three arrays were 0.754, 0.890, and 1.076, respectively). To double check the expression levels of these two genes, western blot was performed and confirmed that PEDF was expressed at high level in normal control corneas but did not change significantly in S-CorNV, while VEGFA showed marginal expression but significantly upregulated in S-CorNV (Figure 3). The situation of VEGFA confirmed the rationale that microarray serves an effective primary screening method, and other complementary methods would be mandatory under some situation, like when certain highly suspected genes gave Marginal or Absent signals in microarray assay.

Functional annotation of differentially expressed genes

The DAVID functional annotation tool was used to annotate genes in various groups (refer to Table 1), mainly based on the GO biological process. The enriched GO terms in each group are summarized in Table 2. While most of the jointly upregulated GO terms (G1 group) were associated with “cell migration,” “defense/immune/inflammation responses,” or “tissue development/organization/homeostasis,” the only two main GO terms enriched in the jointly downregulated group (G7 group) were “oxidation reduction” and “programmed cell death” (Table 2). Trying to reveal common pathways determining development of two apparently different CorNV models, the whole list of involved genes in several promising GO terms from Group 1 are also given. They include “chemotaxis” (Table 3), “immune response” (Table 4), “inflammatory response” (Table 5), and “regulation of cytokine production” (Table 6).

Even without going further into each gene’s function, it could be said that the genes enriched in these groups coincided the events described above (Table 2). Since blood vessel formation is the core component of neovascularization, we looked up all genes that belonged to the “blood vessel development” term (Table 7) among all changed probes. Surprisingly, none of the common proangiogenic or antiangiogenic factors, such as VEGF and PEDF, appear in this table. However, a group of lens crystallins included in the GO term of “structural constituent of eye lens,” which has recently been proven to be expressed physiologically in mammalian corneas [20], was among the upregulated probes (Table 8). Most of these lens crystallin genes were upregulated in S-CorNV at about 2–10 fold higher than in the CB-CorNV model. By contrast, the enzyme crystallins (e.g., ALDH1A1 and ALDH3A1) which have long been proposed to be involved in antioxidation in ocular tissues [21-23], were significantly downregulated (Table 9). Among other genes that were listed as downregulated in the GO term of “oxidation reduction” (Table 9), the clustering of five members of the cytochrome P450 family in this group is quite suggestive, since this family has been proposed to be critical in angiogenesis, and inhibitors of them suppress angiogenesis [24].

Table 10 further lists some genes that were significantly changed during CorNV and deserve further studies. For example, since little information about the gene Corneal Endothelial-specific Protein-1 (NM_026358) is available and its function remains unclear [25], the fact that it was downregulated in both models suggested that this gene might be functionally involved in CorNV or that corneal endothelial cells might be also involved in CorNV development. Besides the genes of same family and that showed concerted changes (e.g., serine peptidase inhibitors), those genes that belong to same family but manifested opposite changes in one or two CorNV models also deserve attention. Detailed comparison of each member of such families will help to dissect the modulation of CorNV development. For instance, while Col3a1 and Col5a2 expression increased and Col4a4 expression decreased in both models, change of keratin 12 was also opposite to that of other genes in its family. Opposite changes were also observed for interferon-induced transmembrane protein 3 and other interferon-induced proteins. Some hypothetical genes that have not been confirmed for any biologic functions by experimental study are also listed for comparison.

Microarray analysis has been used successfully in many studies to screen for potential key molecules during a specific process, such as infectious keratitis [26-28]. Aiming at providing a panorama of gene changes during the fast-growing phase of CorNV, this paper does not try to provide experimental confirmation of any proposed genes. However, we had been successful in identifying certain potential targets for manipulating CorNV based on our discovery using microarray profiling. For example, S100A8 and S100A9 were first observed in this study (Table 3). Considering that S100A family members are important for neutrophil functions, later of which are found important in CorNV development [29], we were able to design experiments and confirm that depleting S100A8 inhibited S-CorNV [10]. We also previously reported that the so-called lens crystallins, for example, the α-, β-, and γ-crystallins, were physiologically expressed in murine corneas [20]. On the other hand, αB-crystallin has been reported to act as a chaperone for VEGFA in angiogenesis [30] and to promote tumor angiogenesis by increasing vascular survival during tube morphogenesis [31]. Since the current study showed that they were significantly upregulated during CorNV, it will be worth investigating the exact role of lens crystallins in CorNV, such as by investigating lens crystallin-deficient animals.

Some genes in the jointly upregulated “immune response” GO terms deserve special attention. For example, the complement pathway, especially the alternative activation pathway, has been proven to be involved in several angiogenic conditions, including ocular neovascularization [32-34], the most famous example being the involvement of complement factor H in the pathogenesis of age-related macular degeneration [35]. The fact that all three types of C1q subunits (α, β, and either γ or C) were upregulated in a similar pattern (Table 4) may suggest that complement pathways were involved in the pathogenesis of CorNV. Besides, C-type lectins, including domain family 4, member d (Clecsf8) and member n (Clecsf10; Table 4), have been proven to be involved in the recognition of pathogens by macrophages . The significance of the upregulation of these genes, as well as of those of two Fc genes (Fc receptor, IgE, high affinity I, gamma polypeptide [Fcer1g] and Fc receptor, IgG, low affinity III [Fcgr3]) in CorNV might be due to the infiltration of various inflammatory cells. Similarly, the study of some genes in other groups, such as semaphorin, 5A (Sema5a) under the “blood vessel development” GO term (Table 7) will also be of interest. Sema5a was recently reported to promote angiogenesis by increasing the proliferation and migration of endothelial cells . However, we have no explanation as to why this gene was downregulated in such CorNV models. Obviously, the genes that are concertedly upregulated or downregulated in both models might suggest common pathways for the pathogenesis of these two models. However, the genes that show differential changes in these two models also deserve attention since they might determine the differential presentation of these two models. Furthermore, this study also provided clues to dissecting the functions of the “predicted genes” that were up- or downregulated.

In summary, we compared the changes of genes expression in two commonly used CorNV models and found that while significant differences existed at the levels of gross presentation and histology, some gene-expression change patterns were shared by these two CorNV models. Further studies are needed to help define the role of some promising genes in CorNV, thereby supplementing our current understanding about the pathogenesis of neovascularization diseases. Taking a step further, such studies might reveal new targets that could be used for manipulation of diseases accompanied with CorNV.