Licorice roots (Glycyrrhiza uralensis, Inner Mongolia of China) and TR fruits (Yueyang, China) were purchased from Zhejiang University of Chinese Medicine Traditional Chinese Medicine Sliced Medicine Co., Ltd. (Hangzhou, China). All these plant materials were authenticated by Prof. Shuili Zhang (School of Pharmaceutical Sciences, Zhejiang Chinese Medical University). The voucher specimens of the licorice roots and TR fruits were deposited in the School of Pharmacy, Zhejiang Chinese Medical University (Nos. 20180903-gc and 20180819s, respectively). Reference standards were purchased from Chengdu Mansite Biotechnology Co., Ltd (Chengdu, China).

Raw and processed TR were used in the study. The processed TR and their extracts were produced based on the previous method (Shan Q et al., 2020). In brief, 6 g of licorice pieces were used for processing 100 g of TR. Licorice pieces were immersed in water and boiled for 30 min. The extract was collected and fresh water was added to the licorice pieces and boiled again for the second extraction. The extractions were combined, concentrated, and filtered. The condensed licorice aqueous extract was added to 100 g of raw TR fruit for full absorption, which were then stir-fried until dry (No.20180819z).

The dried fruit of raw and processed TR (30 kg) were extracted with an 8-fold volume of 70% ethanol for 1 h, twice, under refluxing. The solvent was evaporated by a vacuum rotary system (IKA, Germany) and reduced to 3 g/ml (high dose). The extracts were diluted by distilled water to low (0.15 g/ml) and moderate (1.5 g/ml) doses before use.

The components of raw and processed TR were identified by ultra-high-performance liquid chromatography (UPLC)-quadruple time of flight mass spectrometry (QTOF-MS) (Shan Q et al., 2020) as shown in Supplementary Figure S1, and the main components were determined (Supplementary Table S1) by the method we previously established (Hui et al., 2021).

Male and female Sprague-Dawley (SD) rats (180–220 g) were purchased from Shanghai Laboratory Animal Center (Shanghai, China) and acclimatized for 7 days in the housing condition of the specific pathogen-free environment (22 ± 1°C temperature, 50 ± 5% humidity, 12 h light/dark cycle, standard pellet diet, and water ad libitum) before the experiment. Animal care and use were conducted in accordance with the principles and guidelines of National Institutes of Health (National Research Council, 2011), and the experimental protocol was approved by the Animal Care and Use Committee of Zhejiang Chinese Medical University (No. 20180716-11).

The overall experimental design is shown in Supplementary Figure S2. The rats were orally administrated with TR extracts for the subacute toxicity study. One hundred SD rats were randomly divided into 10 groups (10 rats per group; 5 females and 5 males). The low, moderate, and high doses were selected as 1.5, 15, and 30 g/kg, respectively, which were approximately 3 times of the recommended human dose (the recommended human dose is 2–5 g per day, the conversion factor from human to rat is 6.2 according to the FDA Guidance for Industry https://www.fda.gov/media/72309/download), 30 fold (about one-tenth of LD₅, raw TR: 145.6 g/kg; processed TR: 160 g/kg), and 60 fold (one-eighth of LD₅₀, raw TR: 222.0 g/kg; processed TR: 230.4 g/kg), respectively (Shan Q et al., 2020). Among the 10 groups of rats, 6 groups of rats were treated with low-, moderate-, and high-dose raw or processed TR extracts by oral gavage for 28 consecutive days. Two recovery groups were treated with raw or processed TR (4-weeks high-dose administration), followed by 2 weeks of recovery. The control group and control-R group (the negative control for the recovery group) were treated with water.

The body weight, food, and water consumption of all rats were recorded weekly before and during the treatment using the methods reported in previous studies (Yang et al., 2019). Food and water consumption were calculated accordingly. At the end of the oral administration experiment, the rats were anesthetized with sodium pentobarbital at 0.06 g/kg (i.p.) and blood were collected from heart, and carbon dioxide was applied to sacrifice. Samples and organs were collected for hematological, biochemical, histological, molecular, and metabolomics analyses as described in the following sections.

Approximately 2 ml of blood sample per rat was used in the hematological analysis (Upadhyay et al., 2019). Serum samples were used in the biochemical analysis of aspartate aminotransferase (AST), glucose (GLU), triglycerides (TG), creatine kinase (CK), and lactate dehydrogenase (LDH) based on reported methods (Upadhyay et al., 2019). Heart, liver, lung, kidney, and brain samples were fixed in 10% neutral buffered formalin for 24 h at room temperature and then subjected to standard procedures for hematoxylin-eosin (H&E) microscopic examination (Loha et al., 2019).

Metabolomics analysis was conducted using serum samples to further study the dose-effect and time-course factors of the oral administration of raw and processed TR in rats (Supplementary Figure S2). For the dose-response analysis, serum samples were harvested from male rats in the subacute toxicity study as described above. For the analysis of time-course factor, another set of male SD rats were used. In brief, 90 male SD rats were randomly divided into nine groups (10 rats/group). One group served as the negative control, and six groups were treated with high-dose raw or processed TR (30 g/kg) for 1, 2 or 4 weeks. Two groups of rats were used as the recovery groups, which were treated with high-dose raw or processed TR for 4 weeks, followed by 14 days of clearance. At the end of the experiment, rats in all groups were anesthetized with sodium pentobarbital (0.06 g/kg, i.p.) for blood sample collection.

Blood samples were allowed to stand still for 30 min for coagulation and then centrifuged (3000 rpm, 10 min). Supernatants were collected and stored at −80°C. Next, samples were prepared for the metabolomics study. Serum sample (200 µl) was mixed with 400 µl of cold methanol and vortexed for 30 s to precipitate proteins. After centrifugation at 12,000 rpm for 15 min, 100 µl of the supernatant was transferred to sample vials with inserts for injection. A pooled quality control (QC) sample was generated by combining 10 µl of supernatant from each tested sample. The QC sample was used to equilibrate the analytical platform at the beginning and during the process (one injection per every 10 samples) (Dunn et al., 2012). All samples were kept at 4°C before injection.

The serum samples were analyzed on a UPLC system coupled with a SYNAPT G2 QTOF high-definition MS system (Waters, Milford, MA, United States) and separated on a Acquity UPLC^® BEH C₁₈ column (2.1 × 100 mm, 1.7 µm). The mobile phase was composed of solvent A (acetonitrile) and solvent B (deionized water containing 0.1% v/v formic acid). The detailed chromatography and MS conditions were set based on our previous work (Shan Q et al., 2020).

Raw data were acquired and preprocessed by MassLynx, UNIFI 1.7, and Progenesis QI (Waters, Milford, MA, United States). Data was collected in both ion modes and subjected to signal drift correction by StatTarget (https://www.bioconductor.org/packages/release/bioc/html/statTarget.html) using the suggested parameters (Luan et al., 2018). The converted data were then exported to MetaboAnlayst 4.0 for principal components analysis (PCA), partial least squares discriminant analysis (PLS-DA), orthogonal partial least-squares discrimination analysis (OPLS-DA), volcano plot analysis, and pathway analysis (Chong et al., 2018). A heatmap generation was produced by R (R Core Team, https://www.r-project.org). After normalization, potential biomarkers were selected according to the PLS-DA and volcano plot results with minor parameter adjustment [variable importance in projection (VIP) > 1, fold change (FC) > 2 or < 0.5, and p value after false discovery rate (FDR) < 0.05], and receiver operation characteristic (ROC) curve analysis was performed (Dunn et al., 2012; Chen et al., 2016; Di Guida et al., 2016). In the ROC analysis of individual biomarkers, the optimal threshold was determined by max (sensitivity + specificity), followed by the web-based software instruction. The potential biomarkers were compared with the chemical standards, filtered by interquantile range, normalized by median, log transformed (base 10), and then tested by multivariate ROC curve based exploratory analysis to explore, test, and cross-validate the biomarkers by three machine learning mathematical models [PLS-DA, Random Forest, and Linear Support vector machine (SVM)] using the biomarker analysis module in MetaboAnalyst 4.0 (Chong et al., 2018).

Total RNA was extracted from liver tissue samples and reverse-transcribed to cDNA. QRT-PCR analysis was performed with SYBR Green PCR Master Mix (Roche, Germany) on QuantStudio PCR (Life technologies, United States). The primers synthesized by General Biosystems (Chuzhou, China) are listed in Supplementary Table S2. Gene expression levels were calculated using the 2^−△△Ct method (Schmittgen and Livak, 2008).

Results were expressed as mean ± standard derivation, unless otherwise specified. For group comparison of pharmacological effects, Dunnett T test and Dunnett’s T3 test were applied when the variances among the groups had homogeneity and had no homogeneity, respectively. All statistical analyses were performed with IBM SPSS statistics (IBM Corporation, Armonk, NY) and the graphs were produced by GraphPad PRISM 8.0 (GraphPad Software, Inc., San Diego, CA, United States).

Male and female rats were treated with a broad dose range of raw and processed TR extracts, which were approximately 3, 30, and 60 fold of the recommended human dose, for 28 consecutive days to test the possible subacute toxicity of TR. Two groups were also treated with high-dose raw and processed TR for 4 weeks and placed in recovery for 2 weeks to test the recovery of rats after TR treatment.

Microscopic inspection was conducted to evaluate the changes in the main organs of the male and female rats. As shown in Figure 1A, local hepatocyte necrosis, focal inflammatory cell infiltration, mild bile duct hyperplasia, and partial hepatocyte vacuolation lesions were found in the liver of male and female rats treated with raw and processed TR, but the rats treated with processed TR had some extent of toxicity attenuation in the liver (Supplementary Figure S3A). No other systemic changes in the other main organs of the rats treated with high-dose TR were observed (Supplementary Figure S3B).

The body weight, food consumption, and water intake of rats in different treatment groups were assessed every week for general toxicity evaluation in the subacute toxicity study. Male rats treated with high-dose of raw TR had a remarkable weight gain reduction at days 21 and 28 compared with those in the control group. The female recovery group, after high dose processed TR treatment at day 35, also had weight gain reduction. However, this phenomenon may be caused by individual differences, because no other weight gain reduction was observed before or after day 35. Besides, the female groups treated with high-dose raw and processed TR drank more water and no other significant changes in food and water intake were found (Supplementary Figure S4, Supplementary Table S3). No mortality was observed in any group throughout the experiment.

Organ coefficient is one of the most sensitive toxicity parameters, and remarkable changes may reveal early damage even before morphological changes happen (Piao et al., 2013). As shown in Figure 1B and Supplementary Table S4, the relative liver/weight ratio (liver index) in both males and females administrated with raw TR at high dose were significantly increased when compared with those of control (p < 0.001). Only the liver index of male rats elevated when administered with high-dose of processed TR, but attenuation was observed when compared with the raw TR-treated groups (p < 0.05). Therefore, although male and female rats had increased liver index after the oral administration of high-doses of TR, processed TR had reduced toxicity of the liver in male and female rats compared with raw TR. The relative weights of the kidney, testicle, and epididymis of rats treated with high-dose raw TR significantly increased (p < 0.05, p < 0.05, and p < 0.01, respectively) and processed TR treatment also had reduced toxicity in these organs of male rats when compared with those of raw TR (Supplementary Table S4). In addition, the liver index of rats treated with high-dose TR returned to normal after 2 weeks of recovery (Figure 1, Supplementary Table S4).

Blood samples of male and female rats treated with TR were collected and assessed for hematology and serum biochemical analysis. As shown in Figure 1B, Supplementary Tables S5, S6, high-dose raw TR treatment resulted in the remarkable elevation of platelet (PLT) levels in both female and male groups, whereas those of processed TR did not lead to a PLT increase. Mean corpuscular hemoglobin concentration (MCHC) remarkably decreased in female rats, whereas plateletcrit (PCT) remarkably increased in male and female rats. Another parameter affected by TR was hemoglobin (HGB), which decreased in the female group treated with high-dose -raw TR. Other parameters, such as eosinophil percentage (EOS) and mean corpuscular hemoglobin (MCH), were within normal laboratory ranges, which indicated that no relative histopathological changes occurred (Giknis and Clifford, 2008). Interestingly, although the high dose is 60 times the recommended clinical amount, most of these disturbed parameters returned to normal after the recovery period. Together, although high doses of TR, especially raw TR, could introduce toxicity to the blood system with 28 days of administration, its toxicity can be reversed after the withdrawal of TR treatment.

AST, CK, and LDH levels increased, whereas ALP, GLU, and TG values decreased in the male groups treated with raw TR at different doses when compared with those of control (Figure 1B and Supplementary Table S6). High-dose processed TR did not cause AST elevation in the female and male groups, which indicated that the processing procedure could be a potential solution to reduce TR’s liver toxicity. Similarly, the elevation of CK and LDH in male rats treated with processed TR were lower than those caused by raw TR treatment. A reduction in TG was also observed in the processed TR-treated female and male rat groups, which demonstrated that treating with processed TR has equal or greater TG decreasing effects when compared with those of raw TR. In addition, most of the affected indicators were ameliorated after recovery.

A metabolomics approach was performed to landscape the broad toxic impact of TR treatment to rats under high-dose exposure, and the key enzymes and genes those may regulate the biosynthesis or transformation of metabolites in the most disturbed pathways were validated. The screened potential biomarkers were further evaluated by dose–effect and time-course group data with reference standard confirmation and then tested and cross-validated by training and prediction data sets. The results provided sensitive and specific potential biomarkers for TR toxicity prevention, diagnosis, and treatment.

A total of 14,854 and 11,879 peaks were detected in the serum samples under positive and negative modes, respectively. After whole groups were subjected to signal drift correction (detailed in Methods), 5,595 and 5,643 variables were retained in positive and negative ion modes. All groups of rat serum samples under both positive (Figure 2A) and negative ion mode (Supplementary Figure S5A) were analyzed to identify the metabolic perturbations among different treatment groups. The PCA plot is unsupervised and provides an informative first look at the relationships among multiple groups. The PLS-DA is supervised and can provide the detailed information of discriminations among groups. OPLS-DA is suitable for comparison between two groups and could also provide the most important variables which could explain the difference between groups (Blaise et al., 2021). Figure 2A, Supplementary Figure S5A are the PCA plots of general clustering profiles. The figures indicate the distinct separations of serum samples treated with different doses of TR extracts. The two major groups are clustered separately and the QC groups are clustered tightly in the center. Figure 2B, Supplementary Figure S5B are PLS-DA plots and show the dynamic progress profiles of dose–effect and time–course factors. The figures indicate a clear dose- and time - accumulated developing cycle. Figures 2C,D, Supplementary Figures S5C,D are PCA plots and Figures 2E,F, Supplementary Figures S5E,F are OPLS-DA plots.

Metabolic changes that may be caused by dose, low-, moderate-, and high-dose of raw TR treated to rats were studied. As shown in Figure 2C and Supplementary Figure S5C, clusters of dose response in rats were clearly identified. The low-dose group was closer to the control, whereas the moderate- and high-dose groups were clustered together. The results indicate that the metabolic perturbations caused by TR-induced toxicity are dose dependent.

Based on the dose–effect results, high–dose TR was chosen to conduct a time–course study of TR–induced toxicity. As shown in Figure 2D and Supplementary Figure S5D, the PCA plots show that the groups clustered by treatment time clearly. The recovery group was clustered with the controls, whereas the 1-, 2-, and 4-week administration groups were clustered on the other side. Compared with the controls, the group treated for 4 weeks had the largest perturbations at the metabolic level. Moreover, the clustering result of the recovery group indicated that metabolic disturbance was reversible even in the treatment with 60 fold of the clinically recommended dose. These metabolic observations are in accordance with the serum AST level and liver/body weight alterations (Figure 1B).

Global profiling of serum samples revealed the remarkable metabolic disturbances in the rats treated with raw TR for 28 days (Figure 2 and Supplementary Figure S5). OPLS-DA plots (Figure 2E and Supplementary Figure S5E) and related statistical analysis show the detailed information of TR toxicity. A total of 499 (positive) and 545 (negative) peaks were filtered out in positive and negative modes, respectively, based on the OPLS-DA model and volcano plot analysis of the control group versus the high-dose raw TR group treated for 28 days (VIP value > 1, FC > 2 or < 0.5, adjusted p < 0.05; Table 1). These potential features were further annotated by standard identification procedure, in which the accurate molecular weight, retention time, and MS/MS fragments were compared with the Human Metabolome Database (HMDB), MassBank, Kyoto Encyclopedia of Genes and Genomes (KEGG), and literature records. Sixty-three metabolites were identified under both ion modes (Table 1). Most of them were arachidonic acid and its derivates, bile acid, and cholesterol derivates. Based on the identified metabolites, the pathway analysis by MetaboAnalyst revealed that TR administration resulted in distinct shifts at the metabolic level. Lipid metabolism, amino acid metabolism, metabolism of cofactors and vitamins, and carbohydrate metabolism were perturbed by TR challenging (Figure 2G). Specifically, the perturbed pathways are primary bile acid biosynthesis (P1), steroid biosynthesis (P2), steroid hormone biosynthesis (P3), and tryptophan metabolism pathways (P4). Half of the top 10 perturbed pathways belong to lipid metabolism according to KEGG classification (Figure 2G).

In the subacute toxicity study, we found that the processing method attenuated TR toxicity in histopathological (Figure 1A) and physiological conditions (Figure 1B). We further explored at the metabolic level to find clues for the mode of action of the toxicity attenuation caused by TR processing. We compared the 28-day, high-dose raw, and processed TR treatment groups (Figure 2F) and focused on the metabolites in the lipid metabolism pathways (primary bile acid biosynthesis, steroids biosynthesis pathway, and arachidonic acid metabolism). Interestingly, processed TR treatment alleviated toxicity in the metabolic level (Figure 3). When we focused on the main metabolites in lipid metabolism pathways, we found that the contents of the perturbed metabolites, such as primary bile acid (taurocholic acid), steroids (cholesterol, lathosterol, vitamin D3, calcidiol, and calcitriol), and inflammatory factors (arachidonic acid and leukotriene A4) were remarkably back-regulated when treated with processed TR. The expression of key genes in the pathways were analyzed to test if the toxicity and toxicity attenuation mechanisms were also regulated at the transcription level. As shown in Figure 4, the rate-limiting enzyme of bile acid biosynthesis, Cyp7a1, increased in the transcript level, whereas the biliary excretion of bile acid enzyme, Abcb11 (or bile salt export pump BSEP), decreased considerably. The reduced Acbc11 transcription level, linked with less steatosis and more inflammation after raw TR treatment, was in accordance with previous report (Fuchs et al., 2020). The transcription levels of pro-inflammatory factors, cPla2, Alox5, and Cox2, also increased and these key enzymes are involved in leukotriene biosynthesis and liver inflammation or injury (Enyedi et al., 2016). All these disturbed genes at transcript level were back regulated when treated with processed TR. The results revealed the possible mode of toxicity attenuation via these three pathways.

Based on the identified potential toxicity metabolites, receiver operation characteristic (ROC) curve analysis (the optimal threshold was determined by max (sensitivity + specificity), followed by web-based software instruction, was performed, and valuable features with area under the ROC curve (AUC) values larger than 0.8 were screened. The screened potential toxicity biomarkers were further validated in male groups to observe for histopathological, physiological, and metabolic progress patterns with dose and time dependence. Afterward, only the performance of three biomarkers were fitted with all these factors (Figures 5A–C). As shown in Figures 5A,B, the three metabolites (7α-hydroxycholesterol, calcitriol, and taurocholic acid) can serve as sensitive and specific potential biomarkers for TR-induced toxicity. All these potential markers are involved in lipid metabolism and show dose-dependent changes after TR administration (Figure 5C). The steroid-related metabolites (7α-hydroxycholesterol and calcitriol) remarkably decreased after 7 days of treatment and remained at a low level until withdrawal. The recovery group demonstrated a reverted pattern when compared with the high-dose and TR treatment groups. The third metabolite, taurocholic acid, had an opposite trend during administration compared with the other biomarkers (Figure 5C). The potential toxicity biomarkers were compared with chemical standards, filtered by interquantile range, normalized by median, log transformed (base 10), and applied in multivariate ROC curve-based exploratory analysis to explore, test, and cross-validate the results by three machine learning mathematic models (PLS-DA, Random Forest, and Linear SVM) through the biomarker analysis module in MetaboAnalyst 4.0 (Chong et al., 2018). The whole metabolomics data set size of the control group, high-dose recovery group, and three raw TR treatment groups was 60 (we excluded the processed TR treatment groups, because we could not explain how much the toxicity attenuation effect would be at the metabolic level). We chose the control group and 4-week high-dose TR treatment group to build the prediction models using the three biomarkers with PLS-DA, Random Forest, and Linear SVM modeling methods (with 20 samples). The models were then run on 2/3 of the left-out samples as the test data set (26 samples from the low-, moderate-, and high-dose TR treatments for 1 and 2 weeks and the recovery groups) to test whether we could obtain a clear classification through each prediction model. Cross-validation was then conducted on 1/3 of the left-out sample set (14 samples) according to the instruction of MetaboAnalyst (https://www.metaboanalyst.ca/docs/Faqs.xhtml). The Linear SVM method obtained the highest AUC and accuracy in the test and cross-validation prediction data sets (Figure 5D). Together, the results showed that the validated potential biomarkers were dose- and time-dependent, sensitive, and specific for TR administration even at the early stage, and the TR-induced metabolic perturbations were reversible after recovery. Therefore, the three biomarkers could be applied in TR-induced toxicity prevention, diagnosis, and treatment.