Comparative Analysis Reveals Novel Changes in Plasma Metabolites and Metabolomic Networks of Infants With Retinopathy of Prematurity

Abstract

Purpose: Advances in mass spectrometry have provided new insights into the role of metabolomics in the etiology of several diseases. Studies on retinopathy of prematurity (ROP), for example, overlooked the role of metabolic alterations in disease development. We employed comprehensive metabolic profiling and gold-standard metabolic analysis to explore major metabolites and metabolic pathways, which were significantly affected in early stages of pathogenesis toward ROP.

Methods: This was a multicenter, retrospective, matched-pair, case-control study. We collected plasma from 57 ROP cases and 57 strictly matched non-ROP controls. Non-targeted ultra-high-performance liquid chromatography–tandem mass spectroscopy (UPLC-MS/MS) was used to detect the metabolites. Machine learning was employed to reveal the most affected metabolites and pathways in ROP development.

Results: Compared with non-ROP controls, we found a significant metabolic perturbation in the plasma of ROP cases, which featured an increase in the levels of lipids, nucleotides, and carbohydrate metabolites and lower levels of peptides. Machine leaning enabled us to distinguish a cluster of metabolic pathways (glycometabolism, redox homeostasis, lipid metabolism, and arginine pathway) were strongly correlated with the development of ROP. Moreover, the severity of ROP was associated with the levels of creatinine and ribitol; also, overactivity of aerobic glycolysis and lipid metabolism was noted in the metabolic profile of ROP.

Conclusions: The results suggest a strong correlation between metabolic profiling and retinal neovascularization in ROP pathogenesis. These findings provide an insight into the identification of novel metabolic biomarkers for the diagnosis and prevention of ROP, but the clinical significance requires further validation.

First reported by Theodore Terry in the 1940s, retinopathy of prematurity (ROP) is characterized by the presence of retinal ischemia and retinal neovascularization, which may lead to a retinal detachment. It is a potentially blinding disease that mainly occurs in very low birth weight (VLBW) preterm infants.^1,2 Despite improvement in neonatal intensive care and increased survival rates of VLBW infants in low- and middle-income countries,³ the incidence and severity of ROP are still high.^4,5 A recent epidemiological study reported that approximately 100,000 children with ROP have been diagnosed worldwide.^2,6 

The diagnostic and treatment strategies of ROP have a number of disadvantages. Clinically, the diagnosis of ROP relies on indirect ophthalmoscopy and color fundus imaging, both methods that provide information on pathological changes in the eye. The current therapeutic approaches for ROP include retinal laser photocoagulation and intravitreal injection of anti-vascular endothelial growth factor (VEGF) agents. However, reported risks include damage to the iris and lens and choroidal rupture due to retinal laser photocoagulation, which may also lead to the constriction of visual field and refractive errors in the future.⁷ Injection of anti-VEGF agents in premature infants has raised concerns related to the long-term suppression of systemic VEGF expression, which could cause developmental delays and even inhibit growth of the central nervous system (CNS).^8,9 Furthermore, the current therapeutic methods only aim at those who have already developed ROP, and damage to the eyes at that stage is irreversible. Therefore, it is vital to identify high-risk groups at the early stages and to develop a method to prevent disease progression. 

In humans, most retinal vasculature development is completed by birth. After preterm birth, ROP commences with suppression in the growth of immature retinal vasculature, secondary to oxygen supplementation, loss of growth factors provided in utero, and metabolic dysregulation.¹⁰ In the two-phase hypothesis of ROP, such an acute rise in oxygen tension can stimulate apoptosis of vascular endothelial cells (ECs) and may cause vaso-obliteration via generation of reactive oxygen species (ROS) (phase I). In a subsequent phase, an infant's vaso-obliterated retina undergoes hypoxic/ischemic stress, which triggers a series of events, such as stabilization of hypoxia-inducible factor-1α (HIF-1α) and production of various proangiogenic factors, resulting in neovascularization (phase II).¹¹ Most cells respond to hypoxia by increasing the stability of HIF-1, which in turn induces the transcription of HIF-1–dependent VEGF-A.¹² VEGF-A is critical for blood vessels and neural growth and stability.¹³ In oxygen-induced retinopathy (OIR) models of proliferative retinopathy, relative hypoxia induces high levels of VEGF-A production, causing uncontrolled vessel growth.¹⁴ The metabolic needs of photoreceptors regulate this process. Photoreceptor metabolic alterations can control pathological angiogenesis.¹⁵ It is also possible to modulate the neovascular response in OIR indirectly by manipulating hypoxia sensing or metabolic homeostasis in the entire retina.¹⁶ This can be achieved by interfering with the hypoxia sensing mechanisms in retinal cells in an OIR model. Several compounds that can antagonize HIF-1α have also been shown to reduce neovascularization after intraocular injection in an OIR model,^17–20 associated with Müller cell-specific HIF-1α deletion.²¹ Deletion of HIF-2α, rather than HIF-1α deletion, in retinal astrocytes also reduces neovascularization.^22,23 Apart from measuring oxygen levels, retinal cells can also use the concentration of metabolic intermediates to respond to changes in vascular supply. During the second, hypoxic phase in an OIR model, aerobic respiration is reduced and the Krebs cycle intermediate, succinate, builds up. Retinal ganglion cells (RGCs) respond to succinate build-up by upregulating angiogenic factors, including VEGF, contributing to neovascularization.²⁴ Thus, manipulating metabolic sensing or activity in the retina could also be a promising strategy to reduce hypoxia-induced neovascularization. 

Metabolomics is the qualitative and quantitative assessment of the metabolites in body fluids. The metabolites are downstream of the genetic transcription and translation processes, as well as being downstream of the interactions with environmental exposures; thus, they could be closely associated with the phenotype, especially for multifactorial diseases. In the last decade, metabolomics has been increasingly used to identify biomarkers in disease, and it is currently recognized as a very robust tool with a great potential for clinical translation. Paris et al.²⁵ revealed metabolic dysregulation in ischemic retinopathy. Liu et al.²⁶ reported that purine metabolism and glycolysis were identified as the major disturbing pathways in polypoidal choroidal vasculopathy. Lambert et al.²⁷ implicated the key role of the pyruvate dehydrogenase kinase (PDK)/lactate axis in age-related macular degeneration pathogenesis, and revealed that the regulation of PDK activity has potential therapeutic value in this ocular disease. Recent advances in multiomics approaches should allow the mechanism underlying ROP to be better understood.²⁸ 

Not surprisingly, the retina, with energy consumption over redouble energy expenditure of tumor cells, is one of the most metabolically active tissues.^8,9 Numerous studies have verified that the occurrence and development of retinal diseases are closely correlated with the results of metabolic abnormalities.²⁸ However, we recognized the limitations of simulating human clinical diseases using animal models. In the present study, we performed metabolomics analysis on plasma from ROP and control infants and comprehensively compared alterations in the metabolic profiles between the groups. Our aim was to provide insight into early diagnosis and treatment by identifying specific potential biomarkers and relevant pathways. 

Methods

Study Design

In the present multicenter case-control study, data were collected between April 2018 and October 2019 from neonatal intensive care units (NICUs) of 22 hospitals in Shenzhen, China. This study adhered to the tenets of the Declaration of Helsinki, has been reviewed and approved by the medical ethics committee of all hospitals, and has been registered at ClinicalTrials.gov (ChiCTR1900020677). All participants’ families signed written informed consent forms. Infants were enrolled if they were born preterm, at gestational age (GA) of 20 to 37 weeks. A total of 114 patients fulfilled the inclusion criteria. The presence of ROP was determined by specialized ophthalmologists from Shenzhen Eye Hospital. The classification criteria for the ROP group included severe ROP cases requiring treatment, such as aggressive posterior retinopathy of prematurity (AP-ROP); zone I ROP requiring treatment; or zone II ROP requiring treatment (posterior pole). For the control group, the criterion was infants who had no ROP until 40 weeks’ postmenstrual age (PMA). The control group was strictly matched with the ROP group for birth weight (BW), with no greater than a 200-g difference and GA and PMA with no more than a 2-week difference. Participants were excluded from the study if (1) the preterm infants had received any specific ocular treatment for ROP (including photocoagulation, vitrectomy, and intravitreal injections); (2) the participants’ families requested to withdraw from the study; (3) the participants were diagnosed with congenital metabolic diseases; (4) severe complications occurred, including but not limited to sepsis, necrotizing enterocolitis, neonatal respiratory distress syndrome, and severe metabolic disorders; (5) the participants’ mothers had a history of medication and/or serious diseases, such as gonorrhea, syphilis, or acquired immunodeficiency syndrome during pregnancy; or (6) ROP occurred in a control group participant in subsequent follow-up visits. To reduce interference in preterm infants, all of the blood samples in the ROP group were drawn on the treatment day, before treatment. 

Diagnosis of ROP

The diagnostic and therapeutic criteria of ROP followed the international classification^29,30 and the screening guidelines for retinopathy of prematurity in China published in 2014.³¹ All infants were screened using a binocular indirect ophthalmoscope (Heine Optotechnik, Bavaria, Germany) and a RetCam 3 (Natus Medical, Pleasanton, CA, USA) until the end of the follow-up. Each infant was examined by two experienced retina specialists independently, and the eligibility of participation was confirmed by both specialists. In our study, the severity of ROP was graded according to the criteria extracted from reliable guidelines (Supplementary Table S1). 

Metabolic Profiling

Sample preparation was carried out at Metabolon (Morrisville, NC, USA) as follows:³² Recovery standards were added prior to the first step in the extraction process for quality control purposes. To remove proteins, to dissociate small molecules bound to proteins or trapped in the precipitated protein matrix, and to recover chemically diverse metabolites, proteins were precipitated with methanol under vigorous shaking for 2 minutes (2000 Geno/Grinder; Spex SamplePrep, Metuchen, NJ, USA), followed by centrifugation. The resulting extract was divided into five fractions: two for analysis by two separate reverse-phase (RP)/ultra-performance liquid chromatography–tandem mass spectrometry (UPLC-MS/MS) methods with positive-ion electrospray ionization (ESI) mode, one for analysis by RP/UPLC-MS/MS with negative-ion ESI mode, one for analysis by hydrophilic interaction liquid chromatography (HILIC)/UPLC-MS/MS with negative-ion ESI mode, and one sample reserved for backup. Samples were placed briefly on a TurboVap (Zymark, Hopkinton, MA, USA) to remove the organic solvent. The sample extracts were stored overnight under nitrogen before preparation for analysis. 

Three types of controls were analyzed in concert with the experimental samples: (1) a pooled matrix sample generated by taking a small volume of each experimental sample, which served as technical replicate throughout the dataset; (2) extracted water samples, which served as process blanks; and (3) a cocktail of quality control standards carefully chosen to avoid interference with the measured endogenous compounds, which was spiked into each analyzed sample, allowing instrument performance to be monitored and aiding in chromatographic alignment. Instrument variability was determined by calculating the median relative standard deviation (RSD) for the standards that were added to each sample before injection into the mass spectrometers (median RSD = 6%). Overall variability was determined by calculating the median RSD for all endogenous metabolites (i.e., non-instrument standards) present in 100% of the pooled human plasma samples (median RSD = 8%). Experimental samples and controls were randomized across the platform run. 

UPLC-MS/MS

The UPLC-MS/MS platform utilized an Acquity UPLC (Waters Corp., Milford, MA, USA) and a Q-Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) interfaced with a heated electrospray ionization source and Orbitrap mass analyzer operating at 35,000 mass resolution. The sample extract was dried and then reconstituted in solvents compatible with each of the four methods. Each reconstitution solvent contained a series of standards at fixed concentrations to ensure injection and chromatographic consistency. One aliquot was analyzed using acidic positive-ion conditions and was chromatographically optimized for more hydrophilic compounds. In this method, the extract was gradient eluted from a C18 column (Waters UPLC BEH C18-2.1 × 100 mm, 1.7 µm) using water and methanol, containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% formic acid (FA). Another aliquot was also analyzed using acidic positive-ion conditions and the extract was gradient eluted from the same aforementioned C18 column using methanol, acetonitrile, water, 0.05% PFPA, and 0.01% FA, and was operated at an overall higher organic content. Another aliquot was analyzed using basic negative ion optimized conditions via a separate dedicated C18 column. The basic extracts were gradient eluted from the column using methanol and water, with 6.5-mM ammonium bicarbonate (pH 8). The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Aquity UPLC BEH Amide column, 2.1 × 150 mm, 1.7 µm) through a gradient consisting of water and acetonitrile with 10-mM ammonium formate (pH 10.8). The MS analysis alternated between MS and data-dependent MSⁿ scans using dynamic exclusion. The scan range varied slightly between methods and covered a range from 70 to 1000 m/z. 

Compound Identification, Quantification, and Data Curation

Raw data were extracted, peak identified, and quality control processed using Metabolon hardware and software.³³ Compounds were identified by comparing them to library entries of purified standards or recurrent unknown entities. Metabolon maintains a library based on authenticated standards that contains the retention time/index (RI), mass-to-charge ratio, and chromatographic data (including MS/MS spectral data) on all molecules present in the library. Furthermore, biochemical identifications were carried out based on three criteria: RI within a narrow window of the RI of the proposed identification, accurate mass match to the library (±10 ppm), and the MS/MS forward and reverse scores between the experimental data and authentic standards. The MS/MS scores were obtained based on a comparison of the ions present in the experimental spectrum to the ions present in the library spectrum. 

Peaks were quantified using the area under the curve (AUC) method. The raw area counts for each metabolite in each sample were normalized to correct for variation resulting from instrument day-to-day tuning differences by the median value for each run day; the median was set to 1.0 for each run. This preserved variations among samples allowed metabolites of widely different raw peak areas to be compared on a similar graphical scale. 

Primary Data

The non-targeted metabolomics analysis of 831 metabolites in 114 eligible subjects was performed by the UPLC-MS/MS method (Metabolomic data project ID: CALI-020-0001). Each biochemical index in OrigScale was rescaled to set the median equal to 1. Then, missing values were imputed with the minimum value imputation method. Exogenous substances (e.g., EDTA, molecular constructs) were removed. To minimize the effects of extreme numbers on the results, values outside the mean ± 5 SDs were excluded.³⁴ The original data, code, and statistical values are described in detail in the Supplementary Materials. 

Multivariate Data Analysis and Pathway Enrichment Analysis

Multivariate data analysis was conducted using SIMCA 14.1 (Umetrics, Umeå, Sweden). Orthogonal partial least-squares discriminant analysis (OPLS-DA) was employed to increase the class separation, flatten the dataset, and locate variable importance in projection (VIP).³⁵ The quality of models was validated using two parameters: cumulated R²Y (R²Ycum, goodness of fit) and cumulated Q² (Q²cum, goodness of prediction). A threshold of 0.5 is widely accepted in model classification to classify good (Q²cum ≥ 0.5) or poor (Q²cum < 0.5) predictive capabilities. 

We validated the OPLS-DA model using a permutation test (200 times) to reduce the risk of overfitting and possibilities of false-positive findings. Plots showed correlation coefficients between the original Y and the permuted Y versus the cumulated R²Y and Q². Fitted regression lines were also displayed, which connected the observed Q² to the centroid of the permuted Q² cluster. The model was considered valid³⁶ if (1) all Q² values from the permuted dataset to the left were lower than the Q² values on the actual dataset to the right, and (2) the regression line had a negative value of intercept on the y-axis. To spot moderate and strong outliers, the DModX test and Hotelling’s T-squared test were performed using SIMCA software. VIP values were then singled out by supervised investigations. VIP is a readout of the contribution of each variable on the x-axis to the model. It is summed over all components and weighted to the Y accounted for by every single component.³⁷ Therefore, VIP ranking reflects the contribution of each metabolite to the model. VIP > 0.5 indicates a higher variation when a metabolite was included in the model; VIP < 0.5 indicates that a metabolite played a less important role.³⁸ In the present study, we selected metabolites with VIP > 0.5 as major discriminant metabolites to expand our selection and pathway enrichment. 

Pathway analysis and statistical analysis were undertaken using MetaboAnalyst 4.0 software³⁹ with available Human Metabolome Database identifiers. For pathway analysis, metabolites were mapped to the Kyoto Encyclopedia of Genes and Genomes⁴⁰ metabolic pathways, and quantitative pathway enrichment and pathway topology analysis were carried out. To select the most discriminant pathways, we chose levels of impact and extremely significant difference (impact > 0.3 and P < 0.05, respectively). Locally weighted scatterplot smoothing (LOWESS) was applied to plot an association between metabolites and the severity of ROP. Also, the raw data were transformed to Z scores using the mean and SD. 

Machine Learning

Our simulation experiment was implemented in Python 3.7.7. We first used χ² tests to select important features. Random forest (RF) was further performed for classification. RF is a supervised learning algorithm that can perform tasks related to classification and regression. It makes decisions using a collective approach by creating multiple trees (i.e., a forest).⁴¹ We employed fivefold cross-validation to determine the optimal tuning parameters. The result was the average accuracy of the fivefold cross-validation procedure. 

Statistical Analysis

Demographic characteristics were first determined by QQ- and PP-plots. Normally distributed variables were presented as mean ± SD; otherwise, data were expressed as median and mean ranks. The two-tailed paired Student's t-test and the paired-samples Wilcoxon test were then used to compare variables that were normally and abnormally distributed, respectively, between ROP and non-ROP groups. Following log transformation and imputation of missing values, if any, with the minimum observed value for each compound, the two-tailed paired Student's t-test was utilized to identify biochemicals that differed significantly between groups. Differences were considered significant when P < 0.05, as well as those approaching significance (0.05 < P <0.10). The statistical analysis was performed using SPSS Statistics 25.0 (IBM, Armonk, NY, USA). For all conditions, fold changes in the levels of metabolites compared with the control group were calculated. Further statistical details can be found in the figure legends. 

Results

Demographic Characteristics

The demographic characteristics of the study cohort are presented in Table 1. The GA, BW, and PMA of both groups were strictly matched according to the matching criteria. Statistical differences were identified in delivery mode (cesarean delivery for 33.30% in ROP group vs. 57.90% in control group; P = 0.008), feeding strategy on the day of the blood draw (breast milk for 43.90% in ROP group vs. 87.70% in control group; P <0.001), and the number of leukocytes in the blood (8.61 ± 3.19 × 10⁹/L in ROP group vs. 9.83 ± 3.08 × 10⁹/L in control group; P = 0.047), which showed statistically significant differences. There was no significant significance in BW on the day of the blood draw (P = 0.140), maternal age (P = 0.416), length of hospitalization (P = 0.146), multiple pregnancy (P = 0.843), or sex (P = 0.570) between the two groups. Additionally, we did not find a statistically significant difference between the two groups in terms of the amount of feeding, oxygen saturation or FiO₂, and oxygen inhalation mode on the day of the blood draw. A previous history of administration of oxygen, switch of oxygen inhalation mode, sepsis, and blood transfusion were comparable between the ROP and control groups. Similarly, no significant difference was found in the level of C-reactive protein and number of platelets between the two groups. Noteworthy, none of the participants had received an insulin infusion. 

Table 1.

 

Characteristics of Participants

 

Classification of Metabolites

We first attempted to categorize the total 742 metabolites pulled down from metabolomics analysis. The pie chart in Figure 1 shows the fraction of eight main classes of metabolites, most of which were lipids (47.04%), followed by amino acids (29.11%). However, a relatively smaller fraction of metabolites was related to nucleotides (5.12%), cofactors and vitamins (4.99%), peptides (4.58%), and carbohydrates (4.04%). To compare the changing trend of these metabolite classes between control and ROP groups, we constructed a heatmap (Fig. 2). The levels of most metabolites in lipids, nucleotides, and carbohydrates were elevated in the ROP group, whereas those in amino acids, energy, cofactors, and vitamins fluctuated. It is noteworthy that most xenobiotics accounted for a small proportion of biochemicals that might indicate less external disturbance. Supplementary Table S2 provides a summary of quantitative differences between the groups. A total of 189 metabolites significantly differed, including elevated levels of 89 metabolites in the ROP group and decreased levels of 100 metabolites in the control group. In addition, an approaching significant difference was detected in the levels of other 67 metabolites. 

Figure 1.

 

Pie chart showing fractions of main classes of metabolites. Categories of metabolites were ordered from the highest to lowest fraction as follows: lipids (47.04%), amino acids (29.11%), nucleotides (5.12%), cofactors and vitamins (4.99%), peptides (4.58%), carbohydrates (4.04%), xenobiotics (3.77%), and energy (1.35%).

Figure 2.

 

Heatmap illustrating categories of metabolites in control and ROP groups. The color scales range from bright blue (low ratio) to bright red (high ratio) and represent the relative ratio of Y = log(Y)(intensity) between the two groups. In general, the levels of metabolites in the ROP group were slightly higher than those in the control group, especially in lipids.

Metabolic Pathway Enrichment Analysis

Construction and Validation of the OPLS-DA Model

We further investigated the OPLS-DA score plot (Fig. 3), which revealed a clear and separate clustering between premature infants with and without ROP. Moreover, the OPLS-DA model achieved a R²Ycum of 0.908 and a Q²cum of 0.523, indicating an excellent goodness of fit and reliable predictive capacity. With statistically validation of the corresponding OPLS-DA model by permutation testing (200 iterations) (Supplementary Fig. S1A), we obtained all permuted R² values below or around 0.9 and most permuted Q² below 0. Furthermore, all R² and Q² values were lower than the original values on the right. Notably, the Q² regression line showed a negative intercept at (0, –0.542). These results collectively suggested a valid model that was not likely to be built by chance. To identify strong and moderate outliers, we performed Hotelling's T-squared test (Supplementary Fig. S1B) and the DModX test (Supplementary Fig. S1C), respectively. No strong outlier in the samples could be identified, whereas infant 45 (in the ROP group) showed an evident deviation in the DModX test. 

Figure 3.

 

OPLS-DA scatterplot of samples from ROP and control groups. Samples from the ROP and control groups clearly separated clusters in the OPLS-DA analysis (R²Ycum = 0.908, Q²cum = 0.523). Green dots in the figure represent ROP (n = 57); blue dots represent non-ROP (n = 57).

Contribution Analysis of all Metabolites in ROP

Based on the OPLS-DA model, we drew a contribution plot, which ranked metabolites according to their contribution to the model (Supplementary Fig. S2). For clarity, only the top 78 metabolites are shown in the figure. The top three contributors were isocitric lactone, arabinose, and mannose. To expand our pathway selection, we considered 649 metabolites (VIP > 0.5) for pathway enrichment. 

Pathway Analysis

We used 649 metabolites with VIP values > 0.5 to carry out metabolite set enrichment analysis (Fig. 4). Pathways at the top right achieved a high impact value and a small P value (green box). We primarily enriched nine pathways with a statistically significant difference (impact > 0.3, P < 0.05) (Supplementary Table S3), including arginine biosynthesis (–logP = 11.584, impact = 0.538); histidine metabolism (–logP = 9.828, impact = 0.533); glycine, serine, and threonine metabolism (–logP = 7.671, impact = 0.729); phenylalanine, tyrosine, and tryptophan biosynthesis (–logP = 7.099, impact = 1.000); alanine, aspartate, and glutamate metabolism (–logP = 6.832, impact = 0.741); phenylalanine metabolism (–logP = 5.936, impact = 0.619); beta-alanine metabolism (–logP = 5.373, impact = 0.672); taurine and hypotaurine metabolism (–logP = 5.306, impact = 1.000); and arginine and proline metabolism (–logP = 3.955, impact = 0.376). 

Figure 4.

 

Results of metabolite set enrichment analysis. Pathway impact was determined by topological analysis (x-axis), and the enrichment logP value was adjusted by the original P value (y-axis). Node color is based on the P value, and the node size is based on pathway impact values. The top nine most significantly affected metabolic pathways (impact > 0.3 and P < 0.05) are inside the green box.

Feature Selection and Random Forest Analysis

To narrow down and identify the optimal biomarkers that are useful in differentiating ROP from non-ROP infants, we implemented 742 metabolites and then used the RF-based χ² feature selection method (Fig. 5). We obtained an unbiased estimate of the predictive performance of the RF model, which typically outputs the top 30 metabolites, and the top six metabolites included isocitric lactone (variable importance, 0.0567), creatinine (0.0416), methylsuccinate (0.0413), arabinose (0.0405), arginine (0.0359), and ribitol (0.0348). However, only differences in creatinine, ribitol, orotidine, cis-4-decenoate (10:1n6), 3-(4-hydroxyphenyl) lactate, undecanoate (11:0), kynurenine, and glutamic acid gamma-methyl ester were statistically significant (P < 0.05). Also, ornithine, 10-undecenoate (11:1n1), and picolinoylglycine showed relatively statistically significant differences (0.05 < P < 0.10). 

Figure 5.

 

Illustration of a random forest graph. The vertical axis shows the top 30 metabolites, and the horizontal axis represents the corresponding importance of the selected feature. Different colors indicate different classes of metabolites. Isocitric lactone had the highest feature importance.

From the RF and the pathway enrichment analysis results, we identified six pathways that corresponded to nine metabolites as the top 30 most significantly affected metabolites in the ROP group (Fig. 6). In the alanine, aspartate, and glutamate metabolism, the average concentration of glutamic acid gamma-methyl ester was elevated in the ROP group compared with that in the control group. This indicated a shift in the redox homeostasis and fatty acid oxidation (FAO), suggesting a more antioxidative state. However, in arginine biosynthesis and arginine and proline metabolism, the average levels of ornithine and dimethylguanidino valeric acid were reduced in the ROP group, whereas the average level of arginine was slightly elevated. In addition, in phenylalanine, tyrosine, and tryptophan biosynthesis, indole-3-carboxylate increased but 3-(4-hydroxyphenyl) lactate and kynurenine showed a decreasing trend in the ROP group. 

Figure 6.

 

Correlation heatmap of the top 30 metabolites and pathways selected from random forest analysis. The green # symbols on the left denote metabolites corresponding to enriched metabolic pathways (on the right with green colors 1, 2, 4, 5, 8, and 9). Statistically significant differences are indicated by purple triangles on the left and right sides (adjusted P < 0.05). Additionally, results before and after adjustment for arginine biosynthesis, FAO, and redox homeostasis pathway enrichment analysis are shown.

After adjusting for statistically significant differences (delivery mode, feeding strategy, and leukocyte), the following metabolites of the first 30 metabolites in the RF remained statistically significant. Creatinine, ribitol, and glutamic acid gamma-methyl ester are found as independent risk factors for ROP (Table 2). The higher the level of creatinine, the lower the risk of ROP, indicating that it is a strong protective factor (P = 0.038; Wald = 4.308; odds ratio [OR] = 8.1465E-7; confidence interval [CI] = 1.4484E-12–0.458). Ribitol was noted as a weak protective factor for ROP (P = 0.016; Wald = 5.770; OR = 0.018; CI = 0.001–0.476). The higher the level of glutamic acid gamma-methyl ester, the higher the risk of ROP, suggesting that it is a strong risk factor (P = 0.015; Wald = 5.908; OR = 1.6491E+05; CI = 10.238–2.6563E + 09). Notable metabolites, which were close to statistical significance (0.05 < P < 0.1), showed a strong statistical significance after adjustment, including ornithine, 10-undecenoate (11:1n1), and picolinoylglycine. Therefore, we demonstrated that creatinine, ribitol, glutamic acid gamma-methyl ester, ornithine, 10-undecenoate (11:1n1), and picolinoylglycine are potential biomarkers of ROP (Fig. 7). Taken together, as shown in Figure 6, pathways related to glycometabolism, redox homeostasis, lipid metabolism, and arginine, as well as tyrosine and tryptophan metabolic pathways, were likely to be closely associated with ROP. 

Table 2.

 

Adjusted Metabolites With Statistically Significant Differences

 

Figure 7.

 

Box plot of potential biomarkers of ROP. Creatinine, ribitol, and glutamic acid gamma-methyl ester showed statistically significant differences (P < 0.05); ornithine, 10-undecenoate (11:1n1), and picolinoylglycine exhibited relatively significant differences (0.05 < P < 0.1). Light red represents elevated levels of the metabolites, and light blue indicates reduced levels of the metabolites. Circles represent outliers, and the plus signs represent the means.

Potential Biomarkers for ROP

We further assessed the receiver operating characteristic curve to determine how good the chosen metabolic biomarkers are at predicting ROP (Supplementary Fig. S3). AUC values for the potential biomarkers were as follows: glutamic acid gamma-methyl ester (AUC = 0.642, CI = 0.541–0.743), ornithine (AUC = 0.671, CI = 0.5719–0.7694), creatinine (AUC = 0.572, CI = 0.467–0.677), ribitol (AUC = 0.766, CI = 0.6783–0.8536), 10-undecenoate (11:1n1) (AUC = 0.602, CI = 0.497–0.707), picolinoylglycine (AUC = 0.618, CI = 0.514– 0.722), all greater than 0.50. Notably, the AUC of ribitol was even greater than 0.766, indicating its reliable predictive ability. According to the hypothesis test (H₀ = 0.5), and after adjusting for confounding factors, creatinine, ribitol, glutamic acid gamma-methyl ester, ornithine, 10-undecenoate (11:1n1), and picolinoylglycine were confirmed to be true positive metabolites and were not random (Table 2). In particular, the accuracy of ribitol-based diagnosis was the highest (AUC = 0.766, sensitivity = 71.9%, specificity = 71.9%) among all potential discriminant metabolites. 

To investigate correlations between levels of potential biomarkers with severity of ROP, we employed the proportional odds model with ordinal logistic regression. The results indicated a valid proportional odds assumption (χ² = 12.892, P = 1.000). In goodness-of-fit tests, the deviance test indicated that the model fits well (χ² = 229.829, P = 1.000), which is better than the model with only a constant (χ² =40.073, P < 0.001). We observed that the levels of creatinine and ribitol decreased with the increase of ROP severity (P = 0.018, OR = 0.027, 95% CI = 0.001–0.544; P = 0.011, OR = 0.448, 95% CI = 0.241–0.834, respectively). However, we did not identify a clear association between glutamic acid gamma-methyl ester, ornithine, 10-undecenoate (11:1n1), or picolinoylglycine and ROP severity after adjusting for delivery mode, feeding strategy, and the number of leukocytes. Our study indicated a potential association between reduced creatinine/ribitol levels and the severity of ROP. Nonetheless, the LOWESS curves of creatinine (first rises and then falls) did not perfectly predict the severity of ROP in the proportional odds model, which could be attributed to coarse grading systems, extreme outliers, or a small number of severe ROP infants (two infants with AP-ROP and one infant with stage 4 ROP) (Supplementary Fig. S4). Thus, further experimental validation is essential, and it is important to plan a confirmatory study. 

Furthermore, we also found that more severe ROP among premature infants was related to feeding with a combination of breast milk and parenteral nutrition (5.9-fold, 95% CI = 2.564–13.699, χ² = 17.315, P < 0 .001), a higher spontaneous labor rate (2.3-fold, 95% CI = 1.034–5.115, χ² = 4.168, P = 0 .041), and a higher count of leukocytes (with every 1 × 10⁹/L increase, the risk of ROP increased by 1.1%; 95% CI = 1.001–1.295, χ² = 3.928, P = 0 .048). These findings are consistent with our previous study on targeted analysis of blood spot samples from ROP infants.⁴² The reason for a higher fraction of ROP infants born via spontaneous labor may be related to those who endured prolonged labor and were more likely exposed to lengthened anoxia and exaggerated fluctuation of oxygen during delivery. 

Beyond this, none of the above six potential biomarkers—glutamic acid gamma-methyl ester, ornithine, creatinine, ribitol, 10-undecenoate (11:1n1), and picolinoylglycine—was linearly correlated with the changes in PMA, BW, and GA based on simple linear regression analysis. For the ROP group, glutamic acid gamma-methyl ester R = 0.025, F = 1.43, P = 0.237; ornithine R = 0.003, F = 0.16, P = 0.689; creatinine R = 0.012, F = 0.68, P = 0.413; ribitol R = 0.019, F = 1.04, P = 0.312; 10-undecenoate (11:1n1) R = 0.041, F = 2.38, P = 0.129; and picolinoylglycine R = 0.011, F = 0.63, P = 0.430. For the non-ROP group, glutamic acid gamma-methyl ester R = 0.005, F = 0.29, P = 0.592; ornithine R = 0.007, F = 0.38, P = 0.539; creatinine R = 0.016, F = 0.89, P = 0.350; ribitol R = 0.006, F = 0.32, P = 0.576; 10-undecenoate (11:1n1) R < 0.001, F = 0.01, P = 0.920; and picolinoylglycine R = 0.025; F = 1.38, P = 0.245. We plotted the second-order polynomial fits to serve as a correlation between biomarkers and PMA (Fig. 8). Our study indicated a potential association between glutamic acid gamma-methyl ester, ornithine, creatinine, ribitol, 10-undecenoate (11:1n1), or picolinoylglycine and the risk and severity of ROP, rather than changes in PMA, BW, and GA (Supplementary Fig. S4). 

Figure 8.

 

Nonlinear regression analysis was used to assess correlation between PMA (days) and biomarkers. The curve was obtained by fitting to a two-parameter model in the nonlinear regression analysis. The potential biomarkers—glutamic acid gamma-methyl ester, ornithine, creatinine, ribitol, 10-undecenoate (11:1n1), and picolinoylglycine—would presumably vary with PMA; however, the analysis did not verify the predicted linear relationship. In the figure, red lines represent the ROP group, and black lines indicate the control group.

Discussion

Although there were statistically significant differences in feeding strategy, we concluded that the groups were comparable at baseline for the following reasons. In the ROP group, blood samples were drawn on the treatment day to reduce interference in preterm infants. To prevent aspiration, ROP infants fasted for 6 and 4 hours before and after the treatment. In order to keep normal energy intake, during fasting, an alternative feeding strategy was supplemented by total parenteral nutrition (TPN). Therefore, a certain amount of ROP infants were mixed fed (breast milk and TPN) on the day of sampling. These procedures may partially account for differences in feeding strategy between the two groups; however, the two groups were strictly matched with the BW, GA, and PMA in the same NICU. According to the nutritional guidelines,⁴³ the feeding volume of infant (similar PMA and weight) feeding strategy and feeding volume. In the present study, the amounts of daily feeding and weights on the day of the blood draw were comparable between the ROP group and the non-ROP group (P = 0.524 and P = 0.140, respectively) (Table 1). Therefore, the amounts of total energy intake were the same or close, and for nutritional formulations the amounts of protein, carbohydrates, and lipids were consistent with similar formulation compositions, whatever their feeding strategy. 

Importantly, for premature infants with a suitable weight (more than 2 kg), NICUs generally do not record specific daily energy intake, which is of no practical significance. Thus, the main concern with such infants is with the total intake through feeding and TPN to maintain daily metabolic requirements. The quantities of nutrients (protein, carbohydrates, and lipids) in the total feeding are relatively fixed, so there is no need to measure the specific calories in breast milk. There were statistically significant differences in feeding methods (children with ROP had lower PMA the day of the blood draw, so the proportion of breast-feeding strategies was small). However, there were no statistically significant differences in total intake, total feeding energy, or weight on the day of the blood draw between the two groups because the two groups shared equivalent total energy intake, protein, carbohydrates, and lipids each day. In addition, we adjusted the feeding strategy to balance its influence. 

We should confirm that different feeding strategies are one of the limitations of the present multicenter study. We will design the experiments more rigorously in the next study. In this strictly matched study (the same NICU, GA, and BW), neither PMA nor the weight on the day of the blood draw showed a statistically significant difference between the two groups; therefore, factors influencing nutritional composition were fairly controlled to some extent. This justifies our intention to design a matched-pair, case-control study. 

In the present study, we compared the plasma metabolomic profiles between strictly matched control and ROP premature infants. The results of analysis revealed few altered metabolic pathways, and six major potential biomarkers of ROP were identified, which were reliable in predicting the development and severity of ROP. We will further discuss these potential biomarkers and relevant pathways to illustrate their ability to predict ROP. We constructed a comprehensive map of potential pathways to indicate the development of ROP (Fig. 9). 

Figure 9.

 

Schematics of potential pathways associated with ROP development. Red words indicate increased metabolites in the ROP group; green words indicate decreased metabolites in the ROP group. Red arrows indicate pathways that had more active reactions, and green arrows represent a less-active reactions. The purple background denotes key enzymes of relevant pathways. CPT1, carnitine palmitoyltransferase 1; PEP, phosphoenolpyruvic acid; ASNS, asparagine synthase; GLS1, glutaminase 1; GPD, glucose-6-phosphate dehydrogenase; GSH, glutathione; GGT, gamma-glutamyl transferase; SAM, S-adenosylmethionine; SAH, S-adenosylhomocysteine; ASS, argininosuccinate synthase; ASL, argininosuccinate lyase; NOS, nitric oxide synthase; IDO, indoleamine-2,3-dioxygenases.

We will also further discuss the most distinguished metabolites and pathways after adjustment, because ROP is characterized by its vascular abnormalities, mainly including neovascularization, and we will concentrate on the implications of the identified metabolite changes in the major vascular components (i.e., ultimate effects on ECs). Meanwhile, other retinal cells, such as photoreceptors, RGCs, and Müller glia, play a significant role in retinal neovascularization, and the metabolic associations require further studies. 

Glycometabolism and Redox Homeostasis

ROP is characterized by its vascular abnormalities, such as neovascularization. Increased endothelial sprouting and proliferation are major cellular events causing pathological proliferative retinopathies; therefore, deciphering the molecular mechanisms underlying these early cellular events is significant to further developing novel therapeutic approaches for the prevention or treatment of ROP. The retina is one of the most energy-demanding tissues of the body. Under hypoxic conditions, there is still a considerable fraction (80%–96%) of glucose to be converted to lactate.^44,45 Similar to the retina, ECs produce adenosine triphosphate (ATP) by aerobic glycolysis, regardless of whether oxygen is present, a process that is referred to as the Warburg effect.⁴⁶ Aerobic glycolysis is rapidly accumulated to supplement cellular energy needs (within seconds to minutes), whereas oxidative phosphorylation (OXPHOS) is 100 times slower than aerobic glycolysis.⁴⁷ This prompts the metabolism to adapt to booming energy demands, as long as the ECs have received stimulation of VEGF. Yizhak et al.⁴⁸ confirmed that production of ATP by aerobic glycolysis, rather than OXPHOS, supports migration of ECs. 

In the current study, we found that the levels of glycolytic intermediates, pyruvic acid, pyruvate, and lactate, were higher in plasma of ROP infants, and the levels of tricarboxylic acid (TCA) cycle metabolites, such as citrate, aconitate, and succinylcarnitine (C4DC), were lower. Such differences are indicative of perturbation in aerobic oxidation. In other words, aerobic glycolysis is overactive in infants with ROP. Importantly, these observations are in agreement with the exuberant proliferation of ECs involved in neovascularization in ROP. In addition, a small amount of mitochondrial pyruvate carriers is responsible for the transportation of pyruvate from glycolysis to mitochondria for OXPHOS, the reduction of which would also lead to retinal disease.⁴⁹ 

Increased levels of the phosphate pentose pathway (PPP) and hexosamine biosynthesis pathway provide another line of evidence supporting the suggestion that aerobic glycolysis is more active in ROP. Glucose-6-phosphate enters the PPP under the action of glucose-6-phosphate dehydrogenase (GPD) to promote the synthesis of ribose-5-phosphate, which is necessary for biosynthesis of nucleotides.⁵⁰ PPP is composed of oxidative (oxPPP) branches and non-oxidative (non-oxPPP) branches; the oxPPP branches play an influential role in the activity and migration of ECs.^51,52 OxPPP branches are facilitated by GPD, which is controlled by VEGF.^53,54 In redox homeostasis, elevated levels of cysteine, alpha-ketobutyrate, 2-aminobutyrate, glutathione, and glutamic acid gamma-methyl ester were found in the ROP group, which could lead to the decreased production of relevant metabolites, such as creatinine. These alterations could also be related to the increase of nicotinamide adenine dinucleotide phosphate (NADPH) from the oxPPP branch, maintaining the redox state of GSH,⁵⁴ which is important in inhibiting the production of ROS. In previous animal (mice) studies on OIR, it was shown that a noticeable amount of ROS was produced by NADPH, which was clearly associated with retinal neovascularization.^55-58 Consistent with the increased GSH levels, we found an accumulation of cysteinylglycine in the ROP group compared with that in the control group, suggesting alterations in gamma-glutamyl transferase activity. However, no significant change was identified in the levels of gamma-glutamyl amino acids (except for gamma-glutamyltyrosine). In addition, we detected downregulation in the levels of several amino sugars, such as N-acetylneuraminate and erythronate. These alterations can be converted from metabolites in the PPP, including ribitol, arabitol/xylitol, and sedoheptulose, which were decreased in the ROP group. This may also account for the fact that oxPPP and non-oxPPP branches are interconvertible in the PPP⁵¹ (Fig. 9). 

In summary, these results suggest a disruption in glycometabolism and a redox imbalance in the plasma of ROP infants leading to an accumulation of correlated metabolites, such as NADPH, GSH, and hexosamine, as well as a decrease in TCA cycle-related molecules. Our prediction model also suggests that creatinine, ribitol, and glutamic acid gamma-methyl ester are likely to be potential biomarkers for ROP. 

Lipid Metabolism Signaling Pathway

Compared with the control group, we found that the FAO process in the ROP group was highly active, which was confirmed by the increased levels of most fatty acids; for example, malonylcarnitine (C3DC) could lead to impaired angiogenesis.^59,60 In addition, in the present study we found that the levels of aspartic acid (nucleotide precursor) and glutamate were also elevated in the ROP group. These two amino acids are known to complement each other to produce α-ketoglutarate and to maintain the protein and/or nucleotide biosynthesis required for proliferation of vascular ECs.⁶¹ In accordance with the results of the present study, our previous targeted metabolomics study reported a similar change in C3DC and glutamate in the ROP group.⁴² This may, to some extent, reveal overactivity of aerobic glycolysis and lipid metabolism in the development of diseases. Metabolic disturbances in the state of ROP infants may then alter photoreceptor metabolism, which regulates pathologic angiogenesis,¹⁵ consistent with our results. Meanwhile, the overactivity of aerobic glycolysis results in a decrease in the product supplied to TCA, which in turn leads to low α-ketoglutarate levels (TCA-cycle intermediates). Low α-ketoglutarate levels (as well as low oxygen levels) stabilize HIF-1α protein and increase the production of VEGF-A, inducing pathological vessel formation (Fig. 9).⁶² 

Corresponding to this finding, in omega oxidation we also detected elevated levels of dicarboxylic fatty acids, eicosanedioate (C20DC), and docosadioate (C22DC). Although beta oxidation represents the primary route of fatty acid metabolism, if this process is overwhelmed or impaired then dicarboxylate-generating omega oxidation (i.e., oxidation of the terminal carbon of fatty acid in peroxisomes and endoplasmic reticulum) may be utilized to meet energy demands. A growing body of evidence suggests that suppressed carnitine palmitoyltransferase 1A could induce retinal vascular deficiency and reduce pathological angiogenesis in a model of ROP.⁶⁰ 

ECs rely heavily on aerobic glycolysis,⁵⁵ whereas FAO contributes only 5% of total ATP production in ECs.⁶³ Schoors et al.⁶⁰ demonstrated that FAO would ensure that ECs have enough de novo synthesis of deoxynucleotides during angiogenic sprouting. It also facilitates selective regulation of the proliferation, rather than migration, of ECs via FAO. This is consistent with the fact that endothelial migration depends mainly on glycolysis,⁴⁸ suggesting that different metabolic pathways regulate different functions of ECs in the process of vascular sprouting. Furthermore, Egnatchik et al.⁶⁴ pointed out that systemic FAO blockade with a chemical inhibitor could alleviate excessive angiogenesis in a mouse model of retinopathy of prematurity. Regarding the FAO process, our exploratory analysis indicated that 10-undecenoate (11:1n1) and picolinoylglycine are potential biomarkers for ROP. 

It is also noteworthy that the retina is very rich in lipids, and lipid oxidation is the main source of the energy for the retina.^65,15 Neuronal metabolic needs control both normal retinal vascular development and pathological aberrant vascular growth. Particularly, photoreceptors, with the highest density of mitochondria in the body, regulate retinal vascular development by modulating angiogenic and inflammatory factors. Dyslipidemia contributes to the development and progression of retinal dysfunction in several eye diseases.¹⁰ Fuel demand from retinal neurons, such as photoreceptors, is very high, and photoreceptor metabolic alterations can affect pathological angiogenesis.¹⁵ 

In the present study, we found an aberrant activation of lipid metabolic pathways, such as fatty acid omega oxidation and fatty acid beta oxidation, which may lead to lactate accumulation and ketone body production. Morphological changes in the outer segment of photoreceptor ultimately result in an absence of retinal angiogenesis physiologically and pathologically. Previous studies^10,66 have also found that metabolic dysfunction and dyslipidemia produce deleterious effects on the eye. Dyslipidemia can be characterized by the increased levels of serum total cholesterol, low-density lipoprotein cholesterol, and triglycerides, or a decrease in serum high-density lipoprotein cholesterol concentrations. In premature infants, high triglycerides are associated with increased severity of ROP.⁶⁷ The level of long-chain polyunsaturated fatty acid is also significantly lower in severe ROP in premature infants at postmenstrual age of 32 weeks.⁶⁸ 

Retinal vasculature is critical for the preservation of normal retinal function.⁶² This can be related to the fact that most nutrients (oxygen, glucose, lipids, ions) delivered to neurons of the inner retina originate from the retinal vasculature, with the majority of photoreceptor cells, Müller glia, and retinal ganglion cells being provided by the choriocapillaris, which lies below the retinal pigment epithelium (RPE). Therefore, altered metabolites in the circulation can influence retinal neuronal function, particularly abnormal vascularization. Metabolic profiling in this study provided a comprehensive snapshot of human metabolism in ROP and clarified specific metabolic mechanisms and pathways of ROP. However, the metabolic associations require further validation.^69,70 

Arginine Metabolism Pathway

In the present study, we found a significantly lower level of ornithine in the ROP group. This is likely to be related to the reduced arginase activity and/or the higher ornithine decarboxylase and/or ornithine–oxo-acid transaminase activity. It has been shown that the upregulation of arginase activity is associated with inflammation, oxidative stress, and peripheral vascular dysfunction.^71,72 This finding is in agreement with the early untargeted metabolomics study performed on OIR by Lu et al.,⁷³ and a study on plasma metabolites of ROP infants conducted by Zhou et al.⁷⁴ However, their baseline data were not comprehensive, and the lack of a feeding strategy and administration of oxygen should be noted, making their results unreliable. Besides, their results were inconsistent without adjusting confounding factors, and their method was defective (not an untargeted, globally accepted standard regimen). Additionally, the sample sizes for their studies were extremely small, which could lead to instability. In the present study, we utilized a rigorous matching strategy using detailed demographic and clinical data (Table 1), whereas other studies have used self-reported data or less stringent matching. To make our results more reliable and scientific, we employed a platform developed by Metabolon and multiple approaches, such as machine learning and multivariate data analysis, after adjustment. Therefore, our speculations regarding several potential pathways are reasonable. 

Furthermore, arginase has been recognized as a therapeutic target for the CNS disorders and cardiovascular diseases, due to its role in mediating damage in the retinal neovascularization, as well as targeting the early phase of neovessel growth.^75,76 In mammalian cells, nitric oxide is known to regulate angiogenesis, neurotransmissions, immune responses, and oxidative stress and to activate endothelial progenitor cells.^77,76 Arginine is a substrate of nitric oxide synthase (NOS) and arginases, and these two enzymes counterbalance each other. Therefore, NOS may be activated when arginase activity is inhibited. In addition, during the phase of ischemic retinopathy in OIR, when arginase activity is downregulated, an increase in inducible nitric oxide synthase expression, physiologic angiogenesis, vitreoretinal neovascularization, and impaired retinal function has been reported.⁷⁸ Also, Wang et al.⁷⁹ found that decreased ROS levels led to angiogenesis and migration of ECs, which could be associated with neovascularization in ROP. 

Additionally, the increased ornithine decarboxylase activity could lead to excessive production of proline and polyamines, causing collagen deposition, pathological neovascularization, and vascular fibrosis.⁸⁰ Moreover, proline is indispensable for the synthesis and maturation of collagen, and it is necessary for the maintenance of blood vessel integrity and angiogenesis; additionally, polyamines are vital factors for cell proliferation, ion channel function, and neuroprotection.⁸¹ Therefore, upregulated ornithine decarboxylase activity may be a potential contributor for fibrovascular proliferation and retinal detachment at the late stage of ROP. 

Retinal Cell Interactions in Retinal Neovascularization

Except for retinal ECs, other retinal cells (photoreceptors, RGCs, and Müller glia) play a significant role in retinal neovascularization. A previous study⁸² showed an alteration in the morphology of the photoreceptor outer segment and a loss of neurons and their synapses in the inner nuclear and plexiform layers of the central retina in the OIR model. Specifically, there was an increase in the number of Müller glia and microglia in mice with OIR, and the elevation in the number of these cells correlated with the absence of deep plexus. This indicates that the activity of both macro- and microglia is altered in regions where the deep plexus blood supply is deficient. In addition, it has been postulated that oxygen demands of the photoreceptors cause the development of vascular abnormalities that are characteristics of ROP.⁸³ In photoreceptor-degenerating mice exposed to oxygen in OIR, there are fewer preretinal vascular endothelial cell nuclei and reduced retinal VEGF expression levels versus wild-type OIR controls.⁸⁴ 

Regarding RGCs, they have been found to play an important role in maintaining the normal structure and function of retinal blood vessels. These findings indicate that retinal vascular networks fail to develop in mice due to the lack of RGCs.¹⁶ RGCs also contribute to retinal blood vessels in experimental models of retinal regeneration induced by ischemia reperfusion.⁸⁵ These results strongly suggest that RGCs play a substantial functional role in both physiological and pathological retinal angiogenesis.⁸⁶ Additionally, Müller glia are the major macroglia and retinal-supporting cells that participate in retinal metabolism, function, maintenance, and protection.⁸⁷ An emerging body of evidence has shown that Müller cells are the predominant source of VEGF.^14,88 

Based on the needs of retinal cells, Fu et al.⁶⁹ demonstrated that RPE cells may uptake glucose from choroid and transport it to photoreceptors, where glucose is converted to lactate through glycolysis and then used as a metabolic fuel. Photoreceptors may spare lactate to Müller glia. Müller glia uptake glutamate released from photoreceptors and convert it to glutamine (Fig. 9). However, the metabolic associations among these retinal cells need to further verification. 

Limitations of the Study

Metabolomics aims at the holistic measurement of a noticeable number of metabolites from a biological sample (biological cell, tissue, organ or organism), which can assist in obtaining a better understanding of the underlying biochemical processes that are altered in disease states, toxicological progression, recovery, etc. However, there are some limitations in this study. First, the size of our cohort was not very large (57 ROP cases and 57 non-ROP cases), and participants were limited to China. For a bigger picture of the ROP conditions worldwide, a multinational study involving a larger population should be carried out. We employed a non-targeted UPLC-MS/MS method to detect as many biochemical species as possible; however, the present study concentrated on only annotated compounds (i.e., with identification available on databases, such as Metabolon). Moreover, un-annotated metabolites might be associated with disease, and their omission can introduce bias to the reported results. Furthermore, to understand the causal roles of metabolic biomarkers and disease development and/or progression, it is essential to perform prospective longitudinal studies that allow the evaluation of future disease risks on the basis of multiomics information or other tissues, such as urine and feces. Additionally, in this preliminary study, the reconciliation between distinguishing metabolites and the occurrence of ROP requires confirmation in a replication cohort. 

There are also limitations related to the modeling. To exclude unfit variables from the list of metabolites, OPLS-DA was applied. OPLS-DA scored plots of ROP and non-ROP showed a clear distinction; validation for the OPLS-DA model was stable, and model overfitting was examined in a 200-time permutation test. However, we selected a certain number of VIP metabolites with a higher, closer value. Another model was developed using a random forest machine-learning algorithm, which achieved greater accuracy (83%) for 30 metabolites. Nevertheless, it is difficult to identify a single biomarker. We also presented the panel of overlap between the results reported by each of the two models and adjusted the confounding factors. Finally, we identified six biomarkers that still could not play a determinant role in the classification of ROP and non-ROP when considered singly. Therefore, in future research, multiple biomarker panels may be more informative than a single biomarker, and all biomarkers should be validated. Furthermore, overlapping two models could potentially increase the possibility of bias and lead to some substances being overlooked. Moreover, because this was a multicenter clinical study, the processing environments could not be completely controlled, which could have some impact on the experiment results, although all blood samples of the different centers were run in the same batch. 

Conclusions

To our knowledge, this is the first multicenter, retrospective, case-control, exploratory study on plasma metabolic profiling using non-targeted UPLC-MS/MS to decipher potential biomarkers for ROP. We demonstrated that creatinine, ribitol, ornithine, 10-undecenoate (11:1n1), picolinoylglycine, and glutamic acid gamma-methyl ester were potential biomarkers for ROP. Furthermore, creatinine and ribitol were strongly correlated with the severity of ROP. As these metabolomic changes could be detected between the asymptomatic phase and disease phase, they might provide clues for the discovery of meaningful metabolic pathways, but their clinical significance requires further validation. Overall, the results provided a relatively unbiased and comprehensive map of ROP and non-ROP infants. The overactivity of aerobic glycolysis and lipid metabolism was noteworthy, and a number of potential biomarkers were screened out by a series of methods. These findings provide a basis for additional studies to further investigate the mechanism of ROP and develop a reliable clinical rationale. 

Along with identification of high-risk individuals for prophylactic treatment, these risk signatures have potential value for clinical screening of new interventional measures. Selecting such individuals for participation has the potential to increase the robustness and benefits of clinical screening, in addition to saving medical resources, thereby increasing therapeutic efficacy. 

Importantly, before confirming clinical application of a metabolomic signature, further studies should be conducted; however, blood metabolomic profiles can provide biological information from an ongoing disease. In summary, the present study provides new insights into experimental studies and our understanding of ROP development, as well as a theoretical rationale for screening and therapeutic strategies for ROP. 

Acknowledgments

The authors thank Marcus Fruttiger from Institute of Ophthalmology at University College London for his great advice and assistance. In addition, we acknowledge the support of the Shenzhen ROP Screening Cooperative Group, which includes, but is not limited to, the NICUs of the following hospitals: Shenzhen People's Hospital, Peking University Shenzhen Hospital, Shenzhen Children's Hospital, Shenzhen Maternal and Child Health Hospital, Shenzhen Luohu Maternal and Child Health Hospital, Shenzhen Nanshan Maternal and Child Health Hospital, Baoan Maternity & Child Healthcare Hospital of Shenzhen, Baoan People's Hospital of Shajing Shenzhen, The People's Hospital of Shajing Shenzhen, Huizhou Central People's Hospital, Guangdong Provincial Maternal and Child Health Hospital, Huizhou Sixth People's Hospital, Huizhou Third People's Hospital, Chaonan District Mingsheng Hospital of Shantou, Shantou Second People's Hospital, The Second Affiliated Hospital of Shantou University, Jieyang Huilai People's Hospital, Puning Maternal and Child Health Hospital, Puning People's Hospital, Zhanjiang Central People's Hospital, Affiliated Hospital of Guilin Medical College, and Guilin Maternal and Child Health Hospital. All of the study participants are thanked for their significant participation. We highly appreciate the efforts dedicated by Haibo Peng, PhD; Chaohui Lian, PhD; and Sisi Luo, and a special thank goes to research assistant Panpan Sun. 

Supported by a grant from the Discipline Layout Project for Basic Research of Shenzhen Science and Technology Innovation Committee, China (JCYJ20170817112542555), and by the Medical and Health Projects of Sanming. 

All data that support the findings of this study are available within the article and Mendeley dataset. The supplementary information, Python code, and metabolomic data are available at https://doi.org/10.17632/7r7xt6pf6x. A demographic data file is available from the corresponding author upon reasonable request. 

Disclosure: Y. Yang, None; Q. Yang, None; S. Luo, None; Y. Zhang, None; C. Lian, None; H. He, None; J. Zeng, None; G. Zhang, None