Both standards of metabolites and stable isotope-labeled internal standards (IS) were obtained from Sigma-Aldrich (St. Louis, MO, United States), Cayman Chemical (Ann Arbor, MI, United States), Bidepharm (Shanghai, China), Steraloids (Newport, RI, United States), Cambridge Isotope Laboratories (Cambridge, MA, United States), and Cayman Chemical or Steroids (Supplementary Table S1). Detailed information was given in our previous study (Hu et al., 2020a). Organic solutions (e.g., acetonitrile, isopropyl alcohol, and methanol) of HPLC grade were purchased from Fisher Scientific (Pittsburgh, PA, United States). The modifier of the mobile phase—formic acid—was obtained from Co., Inc. (Fairfield, OH, United States). Ultrapure Millipore water was prepared by a purification system.

Animal experiments were carried out according to the Guidelines for the Care and Use of Laboratory Animals. Rats (n = 6, 6–8 weeks old, 180–220 g) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China) and raised in controlled environment (25 ± 2°C, 40–70% humidity, and 12-h light on/off cycle). Rats were fed with free drinking water and standard feed. All rats were placed in a single rat IVIVC cage. Animals were accommodated for 1 week prior to the experiment.

For the longitudinal PM of remdesivir, the rats were intravenously administered remdesivir (5 mg/kg) dissolved with 12% sulfobutylether-β-cyclodextrin in water. Blood samples were collected from ophthalmic veins by sterile capillary into (NaF/K-Ox) tubes at before (0 h) and after administration (5, 15, and 30 min and 1, 2, 4, 8, 12, 24, and 48 h) and then directly centrifuged to obtain plasma (3,500 rpm, 10 min, 4°C). All plasma samples were retained at -80°C for further analysis.

For the transversal PM, raw metabolomic data were separated into two parts: before (pre-dose) and after (post-dose) administration. The metabolomic profiling and trajectory effect of remdesivir were investigated and analyzed using multivariate statistical analysis methods.

High-performance liquid chromatography-tandem mass spectrometry system (HPLC-MS/MS, Spark Holland; API 5500, SCIEX, Canada) was adopted for targeted metabolomic analysis. The chromatography columns (Waters BEH, HSS T3) and elution solvent (gradient elution) were all evaluated and used according to our previous study (Hu et al., 2020a). The column temperature was set at 20°C with injection volume of 5 μL. Water phase and organic phase (acetonitrile:isopropyl alcohol = 7:2, v/v) contained 0.1% formic acid, and gradient elution was achieved within 27 min.

All analytes were detected via both negative and positive modes with the help of rapid polarity switching and the advanced MRM algorithm. The MS electrospray voltage was 4500 and 5500 V for negative or positive modes, respectively. The optimized MRM parameters are shown in Supplementary Table S1. Detailed parameters of HPLC-MS/MS are shown in our previous study (Hu et al., 2020a).

The one-step protein precipitation method was adopted for this metabolomic analysis. Briefly, an aliquot of 50 μL plasma was spiked with 10 μL IS mixture (eight ISs, 400 ng/ml) and 140 μL precipitation solution (-20°C methanol). Afterward, the mixture was vortexed for 2 min and centrifuged at 13,500 rpm for 10 min at 4°C. The solutions mentioned earlier were injected into the HPLC-MS/MS system for analysis.

For the purpose of ensuring reliable quantitation of all analytes and better comparability in routine analysis, quality control (QC) samples were prepared by pooling equal volumes of unknown plasma. Briefly, six aliquots of pooled QC samples were constructed as real samples and the analytical sequence was interpolated to check the status of sample injection and the HPLC-MS/MS system.

Pharmacokinetic analysis was performed as reported in our previous study (Du et al., 2021b). Briefly, chromatography separation (LC-20ADXR, Shimadzu, Japan) was accomplished on a Waters XBrige C₁₈ column (50 × 2.1 mm, 3.5 μm) using gradient elution. The temperatures of the autosampler and column were set at room temperature and 40°C, respectively. The flow rate was kept at 0.4 ml/min under the gradient elution mode (Du et al., 2021b).

The mass spectrometry parameters (QTRAP 5500, SCIEX, Canada) were used, and the protonated molecule [M + H]⁺ ion was used for all analytes. The quantitative MRM was set at m/z 292.2→163.2 for Nuc and 237.1→194.1 for the IS (carbamazepine). The calibration curve was linear in the range of 2–1,000 ng/ml (Nuc, the active metabolite of remdesivir). A simple and high-throughput protein precipitation method was used for preparing plasma samples (Du et al., 2021b). Method validation of selectivity, sensitivity, accuracy, precision, recovery, matrix effect, stability, and incurred sample reanalysis met the criteria of method validation guidelines. Detailed results were illustrated in our previous study (Du et al., 2021b).

Raw data files were processed and checked by MultiQuant 3.0.1 (SCIEX). The concentrations of analytes were calculated according to the calibration curve. Pharmacokinetic parameters, such as C_max (maximum concentration) and AUC (area under the curve) were calculated using Phoenix (Pharsight 8.3, Mountain View, CA) software. Pearson’s correlation was used to investigate the relation between metabolomics data and PK parameters using IBM SPSS 26.0 (Armonk, New York, United States). SIMCA-P software (v14.1, Umetric, Umeå, Sweden) was used to build mathematic models including unsupervised principal component analysis (PCA), supervised orthogonal projection to latent structures-discriminant analysis (OPLS-DA), and partial least squares (PLS). The compounds with values of variable importance in the projections (VIPs) > 1 and statistical significance of p < 0.05 were picked out for further identification and metabolic pathway analysis. Two hundred random permutation tests were used to check overfitting and random effects, which can assess the predictive ability of the model. Pathway analysis was achieved using online MetaboAnalyst 5.0 (http://www.metaboanalyst.ca), while the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was also used for hierarchical cluster analysis (HCA) and the t-test and mechanism analysis. A p value less than 0.05 was considered statistically significant.

The comprehensive metabolomics method developed in our laboratory was utilized for present quantitation in plasma samples (Hu et al., 2020a). As described in Figure 1, a total of 289 metabolites, which contained amino acids, bile acids, and vitamins, were covered in the present metabolomics method. All biologically active metabolites can be quantitated during the 27-min analysis period. The calculation linearity ranged from 0.2 to 5,000 ng/ml, which provides powerful capability for successful quantitation of low-abundance compounds. Furthermore, other parameters were all carefully investigated (Supplementary Table S1).

With respect to the longitudinal metabolic fingerprints after remdesivir treatment, metabolite peak area ratios of all plasma samples were pooled into the statistical data set. For metabolic data quality analysis, a distance-to-model (DModX) plot was used to check the outliers, and all samples were in the limit of 2 (data not shown). Both PCA and OPLS-DA were used to integrate and co-analyze all data to explore the longitudinal metabolic trajectory in all rats. Individuals were utilized as the grouping basis, and dots of the same color represented samples of one rat at different time points (Du et al., 2021a). From the results of Figure 2A, plasma samples of one rat were divided into tight clusters, which indicated that the longitudinal metabolic fingerprint of the same rat was relatively stable after remdesivir treatment. The metabolic fingerprint changes generated by individuals were greater than the metabolic disturbance induced by remdesivir treatment. This model was validated by 200-time permutation, and no overfitting was observed (Figure 2B, R ²= (0.0, 0.743), Q²= (0.0, −0.835)). HCA was calculated based on the Euclidean correlation with the Ward clustering algorithm (Figure 2D). In order to analyze the metabolic trajectory at different treated time points, the metabolic changes between pre-dose and post-dose samples (5, 15, 30 min, 1, 2, 4, 8, 12, 24, and 48 h) were analyzed and are shown in Figure 2C. The metabolic trajectory was almost steady from 5 min to 24 h; however, distinguished metabolic profiles were shown for that at 48 h. Besides, the number of VIP >1.0 was 80, and the highest VIP was phenylacetylglycine (1.97), l-kynurenine (1.87), and cholic acid (1.83), respectively (Supplementary Figure S1).

For the purpose of exploring the metabolic phenotype variation caused by remdesivir treatment, the metabolomics features at baseline (pre-dose) were compared with those at the treated period (post-dose). For the transversal PM of remdesivir, Figure 3A illustrated that all metabolic data were introduced for OPLS-DA. Although limited metabolic samples were used in the present study at baseline time points, the individuals of both groups were discriminated well in the OPLS-DA model. Moreover, the random permutation test with 200 iterations was performed to investigate the validity and predictability of the OPLS-DA model. Figure 3B shows that no overfitting was observed for all introduced data (R ²= (0.0, 0.354), Q²= (0.0, -0.424)). The most significantly changed metabolites between pre- and post-dose were picked up using the independent samples t-test. A total of 65 metabolites were observed to be significantly changed with the cutoff of VIP >1, p < 0.05 (Supplementary Table S2). The results of pathway analysis indicated that the high-impact pathways were linoleic acid metabolism; phenylalanine, tyrosine, and tryptophan biosynthesis; phenylalanine metabolism; alpha-linolenic acid metabolism; arachidonic acid metabolism; glycine, serine, and threonine metabolism; and arginine biosynthesis. As shown in Figure 3C, although the impact of linoleic acid metabolism was the highest, neither the hits nor the p value met the acceptable criteria. After pathway analysis, two pathways, arginine biosynthesis (two hits) and aminoacyl-tRNA biosynthesis (nine hits), were considered the disturbance metabolism pathway, with impact of 0.228 and 0.167 and p value of 0.046 and 1.26 E-06. The arginine biosynthesis pathway consisted of l-citrulline and l-glutamine; the aminoacyl-tRNA biosynthesis pathway comprised nine amino acids (l-phenylalanine, l-glutamine, l-serine, l-valine, l-lysine, l-isoleucine, l-leucine, l-threonine, and L-tryptophan) (Supplementary Table S3).

As shown in Figure 3D, the significantly changed metabolites correlated with each other positively or negatively. Taken together, inherent metabolic phenotype variations had taken place as a result of the treatment of remdesivir.

Additionally, the metabolic intensity of 65 significantly changed metabolites was compared transversally. The differentially regulated metabolites are presented by fold change (FC > 1: upregulated metabolites; FC < 0.5: downregulated metabolites). The upregulated results (three metabolites) are shown in Figure 4A. Besides, the intensity of the remaining 62 metabolites was downregulated, and the five most significantly downregulated metabolites included two fatty acids, glycohyodeoxycholic acid and glycoursodeoxycholic acid compared between pre- and post-dose treatment (Figure 4B). Figure 4C of the volcano map indicates that 41.54% (27/65) metabolites whose FC < 0.5 were thought to be significantly disturbed after remdesivir treatment.

To further explore whether metabolic fingerprint disturbance accompanies the plasma drug exposure, Pearson’s correlation analysis was conducted to determine what metabolites are highly interplayed with this tendency. As shown in Figure 5, a total of 15 metabolites—12 carnitines, one N-acetyl-D-glucosamine, one allantoin, and one corticosterone—were significantly correlated with the concentration of Nuc (r > 0.5, p < 0.05).

The mean concentration-time curve and pharmacokinetic (PK) parameters of Nuc were elucidated in rats after intravenous administration of remdesivir (5 mg/kg) according to our previous study (Du et al., 2021b). The individual blood concentration-time curves of Nuc measured in six rats presented a high degree of interindividualized variation of PK behavior (Supplementary Figure S1). During the four parameters summarized in table, t_1/2 variation (4.22-fold difference) was the most significant due to individual treatment among all rats. The variation fold change of C_max, T_max, and AUC_0-t was 1.64-, 2.00-, and 3.50-fold difference between the maximum and minimum values. It is known that AUC and C_max can be considered as indicators of drug efficacy or toxicity to some extent (Xing et al., 2019). Therefore, AUC and C_max were chosen for further PLS model analysis.

A supervised PLS model was constructed for capably predicting PK parameters and identifying relationship between two groups of variables (Xing et al., 2019). Thus, the 207 endogenous metabolites were described as one group of variables (X, the predictive variables), while AUC or C_max were represented as a group of variables (Y, the response variables), respectively. First, the PCA model was constructed to find out outliers to avoid deviation of prediction. All rats are distributed according to pre-dose metabolic profiles (Supplementary Figure S2).

Second, the intensities of 207 metabolites were correlated with AUC or C_max of Nuc in the initial PLS model to roughly investigate the relationship between the X and Y variables (Figure 6). The two-component PLS model is adopted for AUC and C_max prediction, which indicates a visible positive linear regression (Figure 6A, R ² = 0.9824; Figure 6B, R ² = 0.9797). Figures 6C, D indicate the loading plot of the aforementioned models and the relationship between predictive variable (X, triangle) and the response variable (Y, box). As shown in this loading plot, X variables on the top right or low left corner represent positive correlation to AUC or C_max and negative correlation to pharmacokinetic response variables. Besides, 88 (AUC) and 83 (C_max) VIP >1.0 X variables were identified due to contribution of X variables to the PLS model (red triangles, Figures 6C, D), which were chosen for following prediction of AUC and C_max, respectively.

Considering the complex and difficult situation to predict PK parameters based on 88 and 83 variables, some significant and representative variables were screened for predictive biomarkers. Pearson’s correlation analysis was used for the association between PK parameters and VIP >1.0 variables. Regarding the prediction model, 16 or 12 VIP >1.0 screened variables were significantly correlated with AUC or C_max, respectively (Table 1). Also, five common metabolites including adenosine, spermine, guanosine, sn-glycero-3-phosphocholine, and l-homoserine were found in both predictive models. As illustrated in Figure 7A, a PLS model was built based on the previous 16 variables, which helps explain about 99.5% variation (R²Y) and predict 98.6% variation (Q²) for the AUC. Simultaneously, as shown in Figure 7B, it could explain about 97.6% variation (R²Y) and predict 95.8% variation (Q²) in C_max based on the 12 variables. Permutation tests were performed with 100 iterations in order to avoid overfitting of this prediction model (Figures 7C, D). Overall, the results indicated that the prediction model exerted ability to predict AUC and C_max with no risk of overfitting. Variables listed in Table 1 were considered potential biomarkers for predicting AUC or C_max.

For the purpose of validating this prediction ability of the aforementioned screened potential biomarkers, all rats were divided into high- and low-value groups according to their AUC and C_max. Discrimination between high- and low-value groups based on the screened biomarkers was performed using OPLS-DA models. As described in Figure 8, the selected 16 biomarkers (AUC model) and 12 biomarkers (C_max model) completely distinguish the two groups, which indicate that these biomarkers are capable of discriminating drug exposure/response.

For the sake of obtaining reliable and reproducible results, several approaches and all sources of fluctuations have been identified and taken to minimize undesirable variation/bias, such as sample handling and preparation and HPLC-MS/MS system status. As shown in Supplementary Figure S3, the peak area ratios of pooled QC plasma samples were clustered under OPLS-DA which indicated that the fluctuation of pooled QC samples was small and constant pertaining to each analyte.