The retina/RPE proteome in chick myopia and hyperopia models: Commonalities with inherited and age-related ocular pathologies

Animals and rearing

Thirty-one male chicks (Leghorn/New Hampshire), obtained from a commercial hatchery, were raised from post hatch days 0–4 under a 12 h:12 h light-dark cycle. On post-hatch day 5, chicks were randomly assigned to a lens condition (+10 or −10 diopters [D], or no lens), and lenses (polymethyl methacrylate) attached to Velcro were fixed to the periocular feathers of the right eye. Separate no lens control animals were used because monocular treatments in chicks can affect blood flow [24,25], refraction, and axial length [26] in the contralateral eye. Lensing was staggered so that the two tissue collection time points were circadian matched; the chicks in group 1 were lensed 3 h into the light cycle, and the chicks in group 2 were lensed 9 h into the light cycle. All procedures adhered to the ARVO Statement for the use of Animals in Ophthalmic and Vision Research and were conducted in accordance with approved La Trobe University Animal Ethics Committee protocols.

Ocular biometry and tissue collection

Following 6 h (group 1) or 48 h (group 2) with the lenses attached, the chicks were anesthetized (ketamine, 45 mg/kg; xylazine, 4.5 mg/kg i.m.), and the right eye refraction and axial dimensions were determined with retinoscopy (Keeler, Vista Diagnostic Instruments, Windsor, UK) and A-Scan ultrasonography (7 MHz probe; A-Scan III, TSL: Teknar, Inc., St Louis, MO). Biometric measures were collected for five chicks per condition, with the exception of the 6 h negative lens condition where biometry could not be performed on one animal due to anesthesia effects and the 48 h no lens condition where six chicks were profiled. The chicks were then euthanized via decapitation, and their right eyes were enucleated. The section of retina/RPE between the ora serrata and 1 mm temporal to the pecten was removed under a surgical microscope (OPMI 1-FC; Zeiss, Jena, Germany), frozen in liquid nitrogen, and transferred to a −80 °C freezer. The time from initial anesthesia to tissue freezing ranged from 5 to 10 min, the order of tissue collection was counterbalanced across lens groups, and all tissue was collected in a single batch within 1 h of the designated time point (6 or 48 h).

Protein extraction, digestion, and LC-ESI-MS/MS

Retina/RPE samples were resuspended in digestion buffer (8 M urea, 50 mM ammonium bicarbonate, and 10 mM DTT) and incubated for 5 h at 25 °C. Samples were then centrifuged for 15 min at 14,000 ×g, and the soluble fraction was used to determine the protein concentrations (Pierce 660 nm Protein Assay, Rockford, IL). Each sample (50 µg protein) was adjusted to 100 µl with digestion buffer, and 55 mM iodoacetamide was then added to the alkylate thiol groups at 20 °C for 35 min in the dark. The alkylated preparation was diluted to 1 M urea with 25 mM ammonium bicarbonate (pH 8.5) before sequencing-grade trypsin (Promega, Madison, WI) was added to the 5 µM final concentration. Digests were performed overnight at 37 °C. The digests were acidified with 1% (v/v) trifluoroacetic acid (TFA) and the peptides desalted on SDB-XC Empore StageTips (3M Company, St. Paul, MN), as previously described [27]. Peptides were reconstituted in 0.1% TFA and 2% acetonitrile (ACN), and 1 µg was loaded onto C₁₈ PepMap 100 µm ID × 2 cm trapping columns (Thermo Fisher Scientific, San Jose, CA) at 5 µl/min for 6 min, and washed for 6 min before the precolumn was switched in line with the analytical column (Vydac MS C₁₈, 3 µm, 300 Å and 75 µm ID × 25 cm, Grace Pty. Ltd., Columbia, MD). The peptides were separated at 300 nl/min using a nonlinear ACN gradient of buffer A (0.1% formic acid and 2% ACN) and buffer B (0.1% formic acid and 80% ACN), starting at 5% buffer B to 55% over 120 min. Data were collected on an Orbitrap Elite (Thermo Fisher Scientific) in Data Dependent Acquisition mode using m/z 300–1,500 as the MS scan range. Collision induced dissociation (CID) MS/MS spectra were collected for the 20 most intense ions per MS scan. Dynamic exclusion parameters were set as follows: repeat count 1, duration 90 s, and the exclusion list size was set at 500 with early expiration disabled. Other instrument parameters for the Orbitrap were MS scan at 120,000 resolution, maximum injection time 50 ms, AGC target 1 × 10⁶, and CID at 35% energy for a maximum injection time of 150 ms with an AGT target of 5,000. The Orbitrap Elite was operated in dual analyzer mode with the Orbitrap analyzer used for MS and the linear trap used for MS/MS.

Biometric data analysis

One-way ANOVA tests using SPSS for Windows (Version 20, IBM Corp, Armonk, NY) were conducted to compare the effects of lens-wear (negative, positive, and no lens) on the refraction, axial length, and vitreous chamber depth at each time point. Tukey honest significant difference (HSD) post-hoc tests were used as required, and p values of less than 0.05 were considered statistically significant. All dependent variables were normally distributed (Shapiro-Wilks test p>0.05 and/or skewness and kurtosis within the acceptable range) and met the assumption of equal variances (Levene’s test p>0.05).

LC-ESI-MS/MS data preprocessing

Identification and quantification of proteins were performed on raw output files from LC-ESI-MS/MS using MaxQuant (Version 1.5.1.6 [28,29]) and the built-in Andromeda search engine with the Gallus gallus UniProt FASTA database (September 2016) and common contaminants. Carbamidomethylation of cysteines was set as a fixed modification, and acetylation of protein N-termini and methionine oxidation were included as variable modifications. Parent mass tolerance was set to 5 ppm after refinement by MaxQuant, and fragment mass tolerance to 0.5 Da. Trypsin was set as the digestion enzyme with up to two missed cleavages allowed. Peptides with a minimum of seven amino-acid length were considered, and the false discovery rate (FDR) was set to 1% at the peptide and protein levels. The match-between-runs feature was used to transfer peptide identifications from one run to another based on the retention time and the mass to charge ratio. The Max LFQ algorithm [29] was used for label-free protein quantification.

Normalized protein group ratios obtained from MaxQuant were log2 transformed. Flagged protein groups (i.e., potential contaminants, reverse proteins, and those identified only by site) and protein groups with more than 40% missing values in any experimental condition were filtered from the results. The remaining missing values (comprising 2.51% of the data) were imputed using the hybrid QRLIC/SVD approach implemented in the imputeLCMD R package (v2.0). The data set was then normalized using the Cyclic Loess function in the Linear Models for Microarray Data R package (LIMMA; v3.30.12 [30]). Protein ANalysis THrough Evolutionary Relationships (PANTHER, v12) [31] gene ontology and protein class classifications for the 1,617 proteins remaining in the data set following filtering, imputation, and normalization are provided in Appendix 1.

Identification of differentially abundant proteins

Empirical Bayes moderated t tests (LIMMA R package v3.30.12) were used to identify differentially abundant proteins in the lens groups relative to no lens groups at each time point and the negative lens groups relative to positive lens groups at each time point. Protein groups displaying a fold change of greater than or equal to 1.5 (log2 FC ≥0.585) and Benjamini-Hochberg [32] adjusted p value of less than 0.05 were considered differentially abundant.

Identification of human orthologs

To enable comparisons with previous studies in other species and the use of curated human gene set files for GSEA, high confidence human orthologs for each Uniprot Accession in the data set were identified using InParanoid (v8.0) [33]. InParanoid employs a reciprocal best hits approach that performs well in benchmarking tests of orthology inference [34]. As orthologs could not be matched for all proteins, the data set used for the subsequent GSEA contained 1,327 identifiers. As a final step, human ortholog Uniprot IDs for the differentially abundant proteins were matched to human Ensembl gene IDs using Biomart (Homo sapiens GRCh38.p7) [35].

Comparison with previous transcriptome and proteomic studies

Differentially abundant proteins were compared with previous exploratory proteomic and transcriptome studies of optically induced refractive error as follows. Lists of genes and proteins implicated in transcriptome and proteomic studies published before July 2016 were obtained from Riddell and Crewther’s supplemental materials [23]. As part of this previous study, the gene and protein lists were mapped to human ortholog Ensembl Gene IDs to enable cross-species and cross-platform comparisons. We performed a PubMed search on June 17, 2017, to update the lists with studies published after July 2016. Of the three studies published since this time [36-38], only two provided lists of implicated genes that could be used in the present analysis [37,38]. The genes implicated in these two studies were mapped to human ortholog Ensembl Gene IDs using the methods reported in our previous publication [23]. We then compared the updated gene and protein lists with differentially abundant proteins from the present study mapped to human Ensembl gene IDs. Chord diagrams were created to display the overlapping gene and protein findings using the GOplot R package (v1.0.2) [39]. Appendix 2 summarizes the methods used in each of the previous studies included in the comparison.

Comparison with the genes and proteins associated with ocular pathologies in humans

We also compared the differentially abundant proteins from the present study with genes and proteins previously associated with ocular pathologies in humans. Author-curated lists of genes associated with primary open angle glaucoma (POAG), choroidal neovascularization, cataract, sub-clinical AMD, and AMD in humans (either via genetic linkage or association or RNA or protein expression changes) were sourced from recent publications and converted to Ensembl Gene IDs using BioMart (Ensembl release 86; Homo sapiens GRCh38.p7) as described previously [22]. The pathology-associated lists were then compared with differentially abundant proteins from the present study using the R software package GeneOverlap (v1.10.0). This package uses Fisher’s exact test to calculate the significance of overlap between each pair of lists in comparison to the genomic background. The function returns the number of intersecting findings between the two lists, the p value, and the estimated odds ratio. The null hypothesis is that the odds ratio is no larger than 1; values larger than this odds ratio indicate a positive association between lists. A contingency table was created to display the pairwise overlaps between protein and pathology lists, including the number of intersecting genes, odds ratios, and significant Benjamini-Hochberg adjusted p values. A chord diagram was created to display the overlapping gene and protein findings using the GOplot R package (v1.0.2) [39].

STRING database analyses

Lists of differentially abundant proteins of particular interest were further explored using the STRING database (v10.5) functional enrichment tool [40]. This analysis identified gene ontologies and pathways that were overrepresented in the protein lists relative to the Homo sapiens genomic background.

Gene set enrichment analyses

GSEA was used to assess Human Phenotype Ontology [41,42] gene set enrichment in the lens groups relative to the no lens groups and in the protein data ranked by correlation with refraction. For the first analysis, a Homo sapiens Human Phenotype Ontology UniProt gene set file was obtained from the Bader laboratory on June 17, 2017 [43], and imported into the javaGSEA app [44] along with the data set of normalized log transformed label-free quantification (LFQ) intensity values mapped to human ortholog UniProt IDs. GSEA was then used to assess Human Phenotype Ontology gene set expression in the lens groups relative to the no lens groups at each time point. The signal-to-noise metric was used to rank proteins by the difference of the class means scaled by the standard deviation, and gene set enrichment at the top and bottom of the ranked lists was determined using 1,000 gene set permutations and an FDR q value cut-off of 0.05.

For the second analysis, correlations between refraction and normalized LFQ protein intensities across the negative lens, no lens, and positive lens groups at each time point were calculated using Spearman’s rho (Hmisc R package). The 6 h and 48 h protein lists ranked by correlation with refraction were then imported into the javaGSEA app [44], and GSEA was conducted using the preranked metric and 1,000 gene set permutations as described above. Leading-edge subset analysis was then used to identify the subset of proteins responsible for gene set enrichment (i.e., proteins that appeared in the ranked list at or before the point at which the running sum reached its maximum deviation from zero).

Refraction and axial dimension changes

The ocular biometry measurements are illustrated in Figure 1A–C. Refractive shifts were evident in both lens groups at 6 h, and by 48 h compensation to the negative and positive 10 D lenses was almost complete (Figure 1A). One-way ANOVAs indicated that refraction was statistically significantly different between the lens groups at 6 h (F(2,11) = 68.215, p<0.001) and at 48 h (F(2,13) = 637.470, p<0.001). Post-hoc tests revealed that the refractive state of all lens groups was significantly different at both time points (p<0.01 for all comparisons). There were no statistically significant differences in the axial dimensions between the lens groups at 6 h (Figure 1B). By 48 h, however, the lens groups differed statistically significantly in axial length (F(2,13) = 8.413, p = 0.005) and vitreous chamber depth (F(2,13) = 8.539, p = 0.004). Post-hoc tests revealed that the axial length was statistically significantly shorter in the positive lens group relative to the negative lens group (p = 0.003) and that the vitreous chamber depth was statistically significantly shorter in the positive lens group relative to the no lens (p = 0.039) and negative lens (p = 0.004) groups. As expected, there was a strong negative correlation between the refraction and ocular axial dimensions (Figure 1C).

Protein abundance changes

The negative and positive lens groups were first compared with age-matched no lens controls to identify differentially abundant proteins at each time point. At 6 h, 144 proteins were differentially abundant in the negative lens group, and 163 proteins were differentially abundant in the positive lens group. At 48 h, 149 proteins were differentially abundant in the negative lens group, and 173 were differentially abundant in the positive lens group. Of these differentially abundant proteins, 72 were either up- or downregulated in both lens groups at 6 h (Figure 2A), and 67 were either up- or downregulated in both lens groups at 48 h (Figure 2B). To identify abundance changes that were sensitive to the sign of defocus, we next compared the negative and positive lens groups directly at each time point. No proteins were differentially abundant between the negative and positive lens groups at 6 h, and only 13 were differentially abundant at 48 h (Figure 2C). Appendix 3 provides a detailed overview of the proteins implicated in each comparison.

Overlap with previous proteomic and transcriptome studies of refractive error

Comparison of these differentially abundant proteins with the findings of previous exploratory transcriptome and proteomic studies of refractive error revealed several commonalities. Eleven of the 390 proteins identified as differentially abundant in the present study have been implicated in previous proteomic studies (Figure 3A). This relatively low number of commonalities may reflect the use of gel electrophoresis in most previous studies, as well as differences in the species, tissues, and time points profiled (see Appendix 2). In contrast, 55 of the 390 differentially abundant proteins have been implicated in previous transcriptome studies where the methodological conditions often matched those of the present study more closely (Figure 3B,C). We also compared our findings with the results of a previous transcriptome meta-analysis, which analyzed the subset of refractive error transcriptome studies conducted in the chick model. This comparison identified another ten unique commonalities with the present proteomic findings (and 12 commonalties in total; Figure 3B). Thus, 65 of the 390 differentially abundant proteins from the present study have previously been implicated at the transcriptome level (53 in individual transcriptome studies, ten in the chick transcriptome meta-analysis, and two in individual studies and the meta-analysis).

We tested each set of commonalities for over-representation of gene ontologies and pathways relative to the Homo sapiens background. No ontologies or pathways were significantly over-represented in the commonalities with previous proteomic studies. A large number of ontologies related to structural processes were over-represented in the transcriptome commonalities (Appendix 4), with the most significant over-representation occurring for the ‘axonogenesis’ ontology. Additionally, several ontologies related to melanin were over-represented in the transcriptome meta-analysis commonalities (Appendix 5). We also conducted an over-representation analysis for the remaining list of differentially abundant proteins that were not implicated in any previous transcriptome studies. Three hundred and ninety-eight ontologies and pathways were over-represented (Appendix 6). The most significant gene ontologies were related to extracellular vesicles, and the most significant pathways were ribosome and spliceosome.

Overlap with pathology-associated genes and proteins

Our previous meta-analysis of chick transcriptome data sets demonstrated overlap between differentially expressed genes in the chick model of early (i.e., ≤3 days) optically induced refractive errors and the genes and proteins previously linked to a range of ocular pathologies for which myopia or hyperopia are risk factors (POAG, choroidal neovascularization, cataract, sub-clinical AMD, and AMD [22];). To test whether a similar association occurs at the protein level, we compared the same pathology-associated gene and protein lists with the differentially abundant proteins from the present study.

As illustrated in Figure 4A, the genes and proteins linked to AMD in humans overlapped statistically significantly with differentially abundant proteins following 6 and 48 h of negative and positive lens-wear in the chick model. The overlap with AMD-associated genes and proteins was greatest at late relative to early time points, and in the negative relative to the positive lens-group (negative 6 h odds ratio [OR] = 8.36, p<0.001; negative 48 h OR = 10.19, p<0.001; positive 6 h OR = 4.49, p = 0.002; positive 48 h OR = 4.79, p = 0.001). The differentially abundant proteins following 48 h of negative lens-wear also overlapped statistically significantly with the genes and proteins linked to choroidal neovascularization (OR = 14.20, p = 0.008). A range of ontologies were over-represented in the protein-pathology overlaps as a whole (Appendix 7). The most significant over-representation was blood microparticles (FDR<0.001). Ontologies related to physiologic stress, cholesterol homeostasis, and melanin were also implicated (Figure 4B).

Correlation of human disease phenotype gene sets with refraction

In addition to comparing the overlap of differentially abundant proteins with specific ocular-pathology associated gene/protein lists, we used GSEA to assess whether any gene sets from the Human Phenotype Ontology database were enriched in response to lens-wear. This analysis tested whether the genes from any Human Phonotype Ontology gene sets were over-represented near the top or bottom of the ranked protein lists. Two methods were used to rank the proteins. In the first approach, categorical comparisons between the lens and no lens conditions at each time point were made by ranking proteins by the difference of class means scaled by the standard deviation. No significant enrichments were identified using this approach. In the second approach, proteins were ranked by Spearman’s correlation with refraction across groups at 6 h and 48 h. No significant enrichments were found in the 48 h analysis; however, at 6 h, abnormal electroretinogram, photophobia, and nyctalopia gene sets were enriched in the proteins negatively correlated with refraction, and the arrhythmia gene set was enriched in the proteins positively correlated with refraction (Table 1).

The leading-edge subsets of the abnormal electroretinogram, nyctalopia, and photophobia gene sets were highly similar and were primarily comprised of genes associated with retinal dystrophies and mitochondrial diseases (Figure 5). Consistent with these disease associations, a range of photoreceptor-related proteins were implicated, including those involved in phototransduction (OPN1LW, GUCA1A, and CNGB1), vitamin A metabolism (RPE65 and RLBP1), photoreceptor disc morphogenesis (FSCN2), protein trafficking within photoreceptors (UNC119), and photoreceptor-controlled retinal organization (RS1). The remaining leading-edge proteins were primarily related to mitochondrial function (TIMM8A, SLC25A4, and MT-CO2) and pre-mRNA splicing (PRPF4 and PRPF8). The final disease phenotype showing a significant correlation with refraction was arrhythmia. Two of the four leading-edge proteins from this phenotype gene set (GNAQ and GNA11) are associated with melanocyte development and melanopsin signaling. Mutations in these genes and NRAS (the third leading-edge protein) have been linked to melanoma in humans. The final protein, A2ML1, is a protease inhibitor involved in innate immunity.