Introduction

Animals in captivity live longer than those in the wild because they are largely protected from predators, internecine competition, and resource scarcity. However, this exposes such animals to diseases of aging that are not encountered at a high frequency in wild populations, including age-related cataracts1. Such cataracts can be heritable, but with significant environmental triggers, including oxidative stress in the lens and the resulting accumulation of DNA and protein damage2,3,4. This is the most common cause of blindness in humans and long-lived companion animals5,6,7 and it now has a prevalence of ~ 20% in the aging captive population of giant pandas (Ailuropoda melanoleuca)8, which live 10–15 years longer in captivity compared to those in the wild9,10.

Genetic analysis has revealed numerous genes associated with cataracts, encoding proteins related to oxidative stress responses, the production of antioxidant enzymes and metabolites, and various DNA repair pathways11,12,13. Previous studies in several mammalian species have investigated the regulation of transcription during the development of cataracts as well as differences in DNA methylation and/or RNA profiles between normal and affected individuals14,15,16. In our previous work, we identified 110 candidate genes that are differentially methylated in giant pandas with and without cataracts (including six genes associated with age-related cataracts in humans)17 as well as more than 700 differentially expressed transcripts18. However, the vast majority of gene products that act as effectors during cataract formation are proteins rather than RNA molecules, and changes at the RNA level do not always correlate directly with changes at the protein level. The proteomic analysis of pandas with and without cataracts may thus provide greater insight into the molecular basis of age-related cataracts in the giant panda population.

We therefore took blood samples from four aged giant pandas with cataracts and six healthy controls for quantitative proteomics analysis using tandem mass tags. We identified proteins that were exclusively or primarily present in healthy pandas or those with cataracts, and assigned them to functional classes and pathways. The analysis of the panda cataract proteome will identify key proteins involved in cataract formation and the resulting pathology, providing new information that could lead to the development of preventative strategies, new diagnostic assays, and non-invasive therapies to reduce the disease burden in the aging panda population.

Materials and methods

Ethical aspects

All samples were collected in accordance with the Wildlife Protection Law of the People’s Republic of China (President of the People’s Republic of China No. 16), and the sampling procedure and subsequent experiments were approved by the Beijing Zoo Academic and Ethics Committee. All methods are reported in accordance with ARRIVE guidelines (https://arriveguidelines.org).

Clinical description of test subjects

Peripheral blood samples (2 ml) were collected from 10 giant pandas (seven females and three males) ranging in age from 21 to 39 (Table 1). Giant pandas 18 or years or older are defined as aged because the equivalent human age is ~ 75. Four of the females (three deceased) were diagnosed with age-related cataracts and the other three females and all three males were healthy, cataract-free controls. We carefully reviewed the case files of each subject to rule out other diseases and traumatic sources of cataract before we selected them for analysis. Importantly, the Beijing Zoo Academic and Ethics Committee insisted that we restrict the number of aged pandas used in the experiments because they are particularly vulnerable. Therefore, the control group was not matched to the age of the cataract group, but instead represented a range of ages and thus helped to rule-out general age-related (rather than specific cataract-related) effects. Blood samples were used because they are easier to obtain than eye biopsies. However, were able to obtain one lens from a deceased panda with cataracts and we analyzed the anterior capsule. The transcriptome data were consistent with the corresponding blood sample (data not shown), which suggested that blood testing is an adequate surrogate for the direct testing of eye material. This approach is much more convenient for the predictive testing of pandas in the future. In the single living female panda with cataracts, the lens capsule had shrunk, the anterior chamber had deepened, the nucleus of the lens had subsided, the cortical granules of the lens had accumulated in the anterior chamber angle and blocked the trabecular meshwork, and we observed lens capsular degeneration and cortical spillage, lens opacification, white coloration, and increased nuclear hardness, similar to human nuclear hardness grade VI. However, this grading was only noted for reference and does not represent the actual hardness of the panda lens. This is because there is no direct correspondence between the cataract disease stages in humans and giant pandas.

Table 1 Description of giant panda sample donors.
Full size table

Protein extraction

The blood samples were centrifuged to prepare serum, which was aliquoted and stored at – 80 °C. For proteomic analysis, 40-µL aliquots of serum were diluted in binding buffer and loaded onto a Pierce Albumin/IgG Removal Kit column (Thermo Fisher Scientific, Waltham, MA, USA). After two rounds of treatment according to the manufacturer’s instructions, the depleted serum was eluted, lyophilized, and stored at – 80 °C. The lyophilized samples were redissolved in 300 µL SDS lysis buffer and centrifuged at 12,000×g for 15 min at 5 °C to remove insoluble particles. After repeating this step, the protein concentration was determined using a bicinchoninic acid (BCA) kit (Shengong Biological Technology, Shanghai, China) and aliquots of each sample were used for SDS-PAGE and quantitative proteomics.

SDS-PAGE

Reconstituted depleted serum aliquots containing 10 µg of total protein were separated by SDS-PAGE on 12% polyacrylamide gels. The gels were post-stained with Coomassie Brilliant Blue (Sigma-Aldrich, St Louis, MO, USA) and the stained gels were documented using an ImageScanner (GE Healthcare, Chicago, IL, USA).

Tryptic digestion and mass tag labeling

Reconstituted depleted serum aliquots containing 10 µg of total protein were digested with trypsin according to the filter-aided sample preparation (FASP) protocol19 with some modifications. Briefly, each protein sample was reduced with 10 mM dithiothreitol, 8 M urea and 100 mM triethylammonium bicarbonate (TEAB) buffer (pH 8.0) for 1 h at 60 °C, then alkylated with 50 mM iodoacetamide at room temperature for 40 min before centrifugal filtration at 12,000×g (10-kDa cut-off) for 20 min at room temperature. The filtrate was diluted in 100 µL 300 mM TEAB buffer and the centrifugation step was repeated twice. The filtrate was transferred to a fresh tube containing 100 µL 300 mM TEAB buffer, and on-filter digestion was performed overnight at 37 °C using 0.1 μg/μL sequencing-grade trypsin (Promega, Madison, WI, USA). The peptides were eluted by centrifugation at 12,000×g for 20 min at room temperature with one change of TEAB buffer, and the final eluate was lyophilized.

For mass tag labeling, the lyophilized peptides were reconstituted in 100 μL 200 mM TEAB buffer and stored at room temperature. The tandem mass tag (TMT) reagent (Thermo Fisher Scientific) was mixed with acetonitrile and centrifuged according to the manufacturer’s instructions. This reagent consists of three chemical groups (reporter, balance and reaction groups) that have the same combined mass (producing a single peak in MS1) but split into different masses during fragmentation (producing separate peaks in MS2 that can be quantified across multiple samples). The peptides and TMT reagent were mixed and allowed to react at room temperature for 1 h before we terminated the reaction with 8 μL 5% hydroxylamine for 15 min. The labeled peptides were lyophilized and stored at – 80 °C.

RP-HPLC–MS

Samples were separated on a Zorbax Extend C18 narrow-diameter column (2.1 × 150 mm, 5 μm; Agilent Technologies, Santa Clara, CA, USA) mounted on an Agilent 1100 HPLC system at a flow rate of 300 μL/min in a gradient of mobile phases A (water containing 2% v/v acetonitrile) and B (90% v/v acetonitrile in water). The gradient elution profile was: 0–8 min, 98% A; 8–8.01 min, 98–95% A; 8.01–38 min, 95–75% A; 38–50 min, 75–60% A; 50–50.01 min, 60–10% A; 50.01–60 min, 10% A; 60–60.01 min, 10–98% A, 60.01–65 min, 98% A. The samples were freeze-dried under vacuum and stored for MS analysis.

The samples were fractionated at a flow rate of 300 nL/min on an Acclaim PepMap RSLC C18 column (75 μm × 150 mm, 2 μm, 100 Å; Thermo Fisher Scientific) mounted on an EASY-nLCTM 1200 liquid phase system (Thermo Fisher Scientific) and were eluted in a gradient of mobile phases A (0.1% formic acid in water) and B (80% v/v acetonitrile in water with 0.1% formic acid). The gradient elution profile was: 0–55 min, 8% B; 55–79 min, 30% B; 79–80 min, 50% B; 80–90 min, 100% B. The separated peptides were introduced into a Thermo Q-active mass spectrometer (Thermo Fisher Scientific) via a nano-liter spray ion source (Thermo Fisher Scientific). The mass resolution of MS1 was set to 70,000, the automatic gain control value was 1 × 106, and the scanning range was 300–1600 m/z. The 10 most intense peaks were scanned in MS2 by high-energy collisional fragmentation in data-dependent positive ion mode, with the collision energy set to 32 eV, the resolution set to 17,500, the automatic gain control set to 2 × 105, the maximum ion accumulation time set to 80 ms, and the dynamic exclusion time set to 15 s.

Proteomic data analysis

Proteomic data were analyzed using Xcalibur v2.1 and were screened against the UniProt panda database using Proteome Discoverer v2.2 (Thermo Fisher Scientific). The false discovery rate (FDR) for peptide identification was held below 1%. The screening criteria to classify differentially expressed proteins were a fold change of at least 0.5 and a significance of p < 0.05. The resulting datasets are available in the ProteomeXchange Consortium repository (http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD031039 or https://www.iprox.cn/page/project.html?id=IPX0004000000; Project ID IPX0004000000, accession number PXD031039). Quantitative bias within samples caused by differences in amino acid sequence length was eliminated by applying the IBAQ method as previously described20. Furthermore, quantitative bias across samples caused by sampling and/or loading errors was eliminated by using the total quantity method, which consisted of (i) adding the IBAQ values of each protein identified in the sample to represent the total amount of protein in the sample, (ii) dividing the IBAQ of each protein in the sample by the total amount of protein in the sample, and (iii) multiplying this value by 100,000 to obtain the normalized protein quantitative value. Variation between the cataract and control groups was compared to within-group variation by one-way ANOVA (SPSS Statistics; IBM, Armonk, NY, USA).

Functional analysis of proteins differing in abundance between cataract and control samples

The TMT quantitative proteomics data were screened for proteins differing in abundance between the pandas with cataracts (samples A1–A4) and healthy controls (samples in bands B or C or B + C). All differentially expressed proteins were then used as search queries against the Gene Ontology and KEGG pathways database, as well as databases of protein interactions, allowing the analysis of comparative expression profiles as well as the construction of heat maps, Venn diagrams, and protein interaction network maps21,22,23. Functional enrichment analysis was carried out to identify differentially expressed proteins significantly enriched in GO terms (biological processes, cellular localization and molecular functions) or KEGG metabolic pathways. Genes were mapped to the GO and KEGG databases, the number of proteins representing each term or pathway was calculated, and hypergeometric tests were applied to identify significantly enriched GO terms or KEGG pathways in the protein list. GO terms and KEGG pathways were considered significant at q < 0.05.

Results

Protein sample preparation and assessment for quality and quantity

Blood samples from 10 giant panda specimens were used to prepare depleted serum proteins for further analysis. To confirm the quality and quantity of protein in the depleted samples, which were prepared by passing serum through an albumin/IgG capture column, we measured the protein concentration using a BCA assay followed by absorbance spectrophotometry, and separated the samples by SDS-PAGE to determine the integrity of the major bands. The protein concentrations in each sample are shown in Table 2.

Table 2 Triplicate absorbance readings (and average reading) and the calculated concentration of proteins in the analyte (diluted for measurement) and the original depleted serum sample.
Full size table

Protein identification and overview of differential expression

LC–MS/MS analysis and database searching resulted in the identification of 342 proteins with qualitative data and 318 with quantitative data meeting the criterion FDR < 1%. We focused on the differences between the pandas in band A (with cataracts) and those in bands B and C combined (healthy controls) in order to maximize the likelihood of recovering differentially expressed proteins reflecting the presence or absence of cataracts rather than those reflecting other differences, such as female vs male. This comparison revealed a total of 148 modulated proteins, 69 upregulated and 79 downregulated in the pandas with cataracts (Fig. 1a). This was based on a fold-change in abundance of ≥ 0.5 and a significance of p < 0.05, as represented in the corresponding volcano plot (Fig. 1b). Among these differentially expressed proteins, 25 of the upregulated proteins and four of the downregulated proteins were characterized by a log2 fold-change > 1 (Table 3 and Supplementary Table S1). There were also 75 upregulated and 96 downregulated proteins in the comparison A vs B (Supplementary Fig. S1) and 81 upregulated and 81 downregulated proteins in the comparison A vs C (Supplementary Fig. S2). A heat map showing the clustering of the proteins in each sample is provided in Supplementary Figure S3. The variation between the cataract and control groups was significantly greater than the variation within either group (one-way ANOVA, F = 7.888, p = 0.000401).

Figure 1
figure 1

Differential protein expression when comparing band A (pandas with cataracts) to bands B and C combined (pandas without cataracts). (a) The number of upregulated and downregulated proteins in band A compared to bands B + C. (b) Volcano plot showing the most meaningful differentially expressed proteins in band A compared to bands B + C by plotting significance on the y-axis against fold-change on the x-axis. Proteins to the left (green dots) are downregulated and those to the right (red dots) are upregulated according to the statistical threshold, below which the proteins are not considered to be differentially expressed (black dots).

Full size image
Table 3 Partial list of differentially expressed proteins (cut-off at |log2 fold change| > 1).
Full size table

Functional annotation of differentially modulated proteins

The biological functions and related pathways of the 148 differentially modulated proteins (A vs B + C) were investigated by screening the sequences against the GO and KEGG databases (Fig. 2). This revealed a total of 2359 enriched biological process classifications as well as 252 cell components, 386 molecular functions and 72 KEGG pathways. When the hits were restricted by a significance of p < 0.05, the number of enriched biological processes fell to 808, as well as 72 cell components, 178 molecular functions and 14 KEGG pathways. The top 10 enriched biological processes, cell components and molecular functions based on significance are shown in Fig. 3. Interestingly, most of the top-ranking functions related to proteolysis and its regulation, which is a significant factor in the age-related deterioration of mammalian tissues and would be expected in the context of cataract formation. The top-ranking molecular functions were in broad agreement, with most of them directly involved in the regulation of peptidases or classified more generally as regulators of enzyme activity. Remarkably, the most significant cell compartment classifications indicated extracellular activity and represented more than half of all the proteins, suggesting that cataract formation in pandas is strongly associated with the extracellular regulation of protein turnover. The inclusion of wounding responses among the biological functions may indicate the response to oxidative stress associated with aging and the formation of cataracts. In line with these results, the analysis of enriched KEGG pathways revealed the proteasome in the top-ranking position by a long margin (Fig. 4a) and also the pathway with the most significant Rich factor (Fig. 4b). The proteasome is the key regulatory complex responsible for protein turnover and its gradually declining activity is known to be associated with aging, thus providing a potential explanation for its link with the formation of cataracts.

Figure 2
figure 2

Number of enriched GO biological process, cell component and molecular function classifications and KEGG pathways defined by the differentially regulated proteins in the A vs B + C comparison. The bar chart shows the number of enriched categories with no statistical cut-off (blue) and with a threshold of p < 0.05.

Full size image
Figure 3
figure 3

GO enrichment analysis showing the top-10 biological processes, cell compartments and molecular functions enriched among the proteins that differ in abundance between pandas with cataracts (band A) and those without cataracts (bands B + C). The enriched categories are arranged left to right by p-value (more significant on the left) and the bars represent the percentage of all differentially regulated proteins included in the category.

Full size image
Figure 4
figure 4

KEGG pathway enrichment analysis. (a) The top-10 pathways enriched among the proteins that differ in abundance between pandas with cataracts (band A) and those without cataracts (bands B + C). The enriched pathways are arranged left to right by p-value (more significant on the left) which is also represented by the logarithmic scale of the y-axis. The dotted lines represent cut-offs at significance values of p < 0.05 and p < 0.01. (b) Significant enrichment functional scatter plot in which the y-axis represents functional annotation, and the x-axis represents the Rich factor of that function (the number of differential genes enriched in a pathway divided by the total number of genes annotated in that pathway). The p-value is represented by the color of the dots, and the number of differentially expressed genes representing each function is shown by the size of the dots. Only the 14 highest enrichments are shown.

Full size image

Protein interaction analysis

Interactions between the differentially expressed proteins and KEGG pathways were visualized by integrating the data using the online platform Omicsbean to generate the network diagrams shown in Fig. 5. As well as showing that many of the differentially modulated proteins interact directly with each other, the diagrams revealed that the downregulated proteins were responsible for most of the predicted interactions and were almost exclusively responsible for interactions with the proteasome. This suggests that one of the key events during age-related cataract formation in pandas is the depletion or loss of proteasome components, which would contribute to the inability of the proteasome to maintain protein integrity in response to oxidative stress.

Figure 5
figure 5

Protein interaction network visualized using the Omicsbean platform. Circles around the edge of the diagram represent the differentially regulated proteins (red = upregulated, green = downregulated, with different shades representing the fold change). The squares represent biological processes, cell localization, molecular functions or signaling pathways, and the color code indicates significance (yellow = low, blue = high, with deeper shades indicating greater significance). Solid lines represent protein–protein interactions, and dotted lines represent the involvement of proteins in metabolic/molecular pathways or functions.

Full size image

Discussion

Aging humans and other mammals can develop cataracts due to the accumulation of oxidative damage in the lens, and this is associated with a declining quality of life1. In captive giant pandas, which can live 10–15 years longer than animals in the wild3,4, the prevalence of cataracts in the aged population (≥ 19 years old) has now reached ~ 20%. The development of age-related cataracts is mainly promoted by environmental factors that cause oxidative damage to DNA and proteins5,6,7,9. This triggers changes in DNA methylation and gene expression, and in previous studies we and others have identified target genes involved in apoptosis, DNA repair and oxidative stress responses that may represent an attempt to counter or reverse the damage that occurs during cataract formation17,18,24,25,26.

Epigenetic modifications such as DNA methylation can directly affect transcription, and we have previously reported correlations between the methylated loci detected in panda cataracts and the modulation of gene expression in affected vs unaffected individuals17,18. However, the regulation of gene expression at the level of protein synthesis and turnover means that transcript and protein levels often do not correlate directly, and this is particularly relevant in the lens where the proteome is known to have a very slow turnover compared to non-lens tissues in the eye, suggesting that serum markers of cataract formation may not track closely with mRNA levels27. We therefore used TMT quantitative proteomics to compare serum from affected and unaffected pandas in order to directly investigate the proteins that change in abundance during cataract formation.

The comparison of affected females (four specimens) with unaffected males and females as a single group (six specimens) revealed 148 differentially modulated proteins, 69 upregulated and 79 downregulated in the pandas with cataracts. Functional annotation of these proteins revealed a broad range of associations, but the most significant were related to the biological processes of proteolysis and its regulation, and molecular functions involving the regulation of peptidases and other enzymes. This was supported by KEGG pathway enrichment analysis, which linked many of the differentially expressed proteins to the proteasome, a protein complex that removes damaged and otherwise surplus proteins by proteolysis28. This link was further supported by protein interaction network analysis, which clearly showed the proteasome as a core interaction hub for the differentially expressed proteins, especially those that were downregulated in the affected pandas. A survey of the list of differentially expressed proteins revealed nine known proteasome subunits, eight of which were downregulated (PSMA7, PSMB1, PSMB7, PSMB4, PSMA6, PSMA4, PSMA3 and PSMB8) while only one was upregulated (PSMA2). The list also included a large number of peptidase inhibitors, most of which were downregulated (e.g., APOA2, SERPINC1, SERPING1, ITH1, ITH2, KNG1, FETUB, SERPINA1 and AHSG) which also indicates that cataract formation involves the dysregulation of the protein turnover apparatus. Abnormalities affecting the proteasome that allow the accumulation of protein aggregates have been linked to human cataracts, including genetic variants of PSMC329, but the overexpression of the proteasome activator PA28αβ did not rescue induced cataracts in mice30. This suggests the proteasome is already highly active in lens tissue and its activity cannot be enhanced using activators, but the loss of activity increases the likelihood of cataracts, as previously reported31,32.

One of the most striking results from the GO enrichment analysis was the overrepresentation of extracellular proteins, accounting for more than 50% of all the differentially regulated proteins we identified. It is well known that one of the key functions of the proteasome is to remove advanced glycation end-products (AGEs), which can form both inside and outside the cells, contributing to the morbidity associated with aging, including cataracts33. Such products also progressively inhibit the activity of the proteasome, leading to a positive feedback cycle of proteotoxic stress that is ultimately pathogenic34. The enrichment we observed for extracellular proteins may reflect the lens attempting to deal with the accumulation of extracellular AGEs during cataract formation. One well-known example of a lens protein that forms AGEs is αA-crystallin, although this accumulates inside cells35, but lens capsule proteins are also glycated upon aging, and thus contribute to cataract formation36. Several of the differentially regulated proteins we detected were annotated as proteins that bind RAGE receptors, including S100A12, S100AB, S100A9 and PSMA2, all of which were upregulated. However, we note that the preponderance of extracellular proteins in our differentially regulated dataset could have a more prosaic explanation—that extracellular proteins are more likely to enter the bloodstream than intracellular proteins are therefore more likely to be captured as serum markers. This may also explain the detection of proteins related to the coagulation pathway such as plasma kallikrein. However, several of the differentially regulated extracellular proteins we detected are components of the extracellular matrix that have previously been associated with cataract formation37. For example, collagen 1A2 and the integrin-binding protein vitronectin were both downregulated in the affected animals (along with the intracellular integrin-binding protein TLN1), supporting the reported role of the extracellular matrix and integrin signaling in lens development and cataract formation38,39. Interestingly multiple keratin proteins were upregulated in the animals with cataracts (KRT18, KRT75, KRT1 and KRT10) whereas various cytoskeletal components and their interaction partners were downregulated, suggesting a complex network of signaling inside and outside the cell as a stress response to the accumulation of protein aggregates.

Although the above results paint an intriguing picture of the panda proteome related to cataract formation, it is important to highlight some weaknesses of the study that may influence the results. First, we acknowledge that the sample number is small, which reflects the limited availability of samples and the further limitations imposed by the ethical committee. This has the potential to introduce bias into the results causes by undetected disease in the test subjects despite our careful and strict inclusion criteria. Second, almost all recorded instances of captive giant pandas with age-related cataracts are female—only two male cases have been reported, one deceased with no material available and one living specimen. This has the potential to introduce bias into the results caused by sex-dependent factors, although we addressed this to a certain extent by dividing the control group into sex categories to enable female affected vs female unaffected and male unaffected vs female unaffected comparisons. These comparisons revealed some differences between males and females but not in the same pathways that discriminated between affected and unaffected animals, suggesting there is negligible sex-dependent interference.

In conclusion, we have generated multiple lines of evidence based on quantitative proteomics showing that age-related cataracts in pandas involve the dysregulation of the lens proteome resulting in the detection of nearly 150 positive and negative protein markers in the blood. More than half of these markers are extracellular proteins, suggesting either that cataract formation causes the extensive modification of the extracellular proteome or that extracellular proteins are more likely to enter the bloodstream in sufficient quantities for detection. One of the most prominent aspects of the cataract serum profile was the depletion of proteasome components and their interaction partners, supporting previous results showing that the proteasome maintenance system deteriorates with age and is progressively less responsive to the accumulation of protein aggregates. Such aggregates may then exacerbate the problem by directly inhibiting the proteasome. Interestingly, we did not identify a single transcription factor among the differentially regulated proteins we detected, which indicates that the response to cataract formation is largely post-transcriptional. The identification of multiple protein markers correlating with the presence of cataracts could allow the development of bioassays for the early detection of cataracts in captive animals, based either on the most profound changes in abundance of key proteins such as those listed in Table 3, or the analysis of multiple biomarker profiles to define cataract-positive patterns, as shown for other diseases40,41. In terms of clinical applications, proteins that increase in abundance during the formation of cataracts could be tested as new drug targets, whereas those depleted during the formation of in cataracts could be evaluated as candidates for replacement therapy, leading to new opportunities for the prevention and/or treatment of cataracts in aging captive pandas.