Herbs of JFAD (Coicis Semen, Salviae Miltiorrhizae Radix et Rhizoma, Pinelliae Rhizoma, Fritillariae Thunbergii Bulbus, Persicae Semen, Cremastrae Pseudobulbus Pleiones Pseudobulbus, Phragmitis Rhizoma, Gecko, Pseudostellariae Radix, Arisaematis Rhizoma.) were tested and approved by the First Affiliated Hospital of Guangzhou University of Chinese Medicine.

The herbs were soaked in water for 30 min before starting the decoction process. Initially, Arisaematis Rhizoma and Pinelliae Rhizoma were boiled for 1 h before adding the other herbs. The mixture was then strained and boiled again for an additional hour. The obtained extracts were subsequently filtered and combined. The herbal blend was concentrated to 11.0 g/mL. This decoction process was conducted by a qualified technician at the traditional Chinese medicine pharmacy of the First Affiliated Hospital of Guangzhou University of Chinese Medicine.

For a 70 kg adult, the recommended clinical daily dosage of JFAD is 2.36 g/kg, totaling 165 g. In an initial experiment, JFAD was mixed with double-distilled water to achieve concentrations of 1.375 g/mL, 5.5 g/mL, and 11.0 g/mL of the crude drug (Peng et al., 2024). These concentrations correspond to one times (low dose), four times (medium dose), and eight times (high dose) the adult dosage, and served as the daily dosages for mice. The resulting extracts were collected in centrifuge tubes and stored at −20 °C for future use. Research confirmed that the TCM decoction remained stable and maintained its quality under these conditions throughout the study (Yu et al., 2022; Liang et al., 2023). The decoction was analyzed using HPLC fingerprinting (Peng et al., 2024).

SPF-grade BALB/c nude mice (5 weeks, 15–17 g) were utilized for the study. The experimental animals were obtained from Zhejiang Viton Lihua Laboratory Animal Technology Co., Ltd., under License No. SCXK (ZHE) 2019-0001. The NSCLC cell line A549 was from the Cell Bank of the Chinese Academy of Sciences (Beijing, China), the NSCLC cell line A549-luc was from the OBiO Technology Company (Shanghai, China).

A549-luc cell suspension (approximately 1 × 10⁶ tumor cells) was prepared by mixing with Corning Matrigel (354248) in a 1:1 volume ratio, resulting in a protein concentration of 18 to 22 mg/mL. The suspensions were stored at 4 °C. Twenty-four BALB/c nude mice were randomly selected, and 0.2 mL of the cell suspension (approximately 1 × 10⁶ cells) was injected into the subcutaneous tissue of the right axilla following sterilization. Twenty-four hours post-inoculation, the tumor-bearing mice were randomly divided into four groups of six: the untreated group, low-dose JFAD group (JFAD-L), medium-dose JFAD group (JFAD-M), and high-dose JFAD group (JFAD-H). Additionally, three mice not inoculated with tumor cells comprised the normal group (NG). JFAD was administered daily via gavage at a volume of 0.3 mL for 21 days, at the same time each day. Mice in the untreated group received an equivalent volume of sodium chloride solution in the same manner. Group details are provided in Table 1.

At the conclusion of the experiment, the mice were euthanized, and blood and stool samples were collected for analysis. Approximately 1 mL of whole blood was obtained from each nude mouse. Blood samples were equilibrated at room temperature for 2 h. After centrifugation, the serum was collected, filtered, and stored at −80 °C. Fecal samples were collected and stored at −80 °C.

Bacterial analysis was performed on fecal specimens stored at −80 °C. Genomic DNA was extracted using the SDS method, assessed for concentration and purity via 1% agarose gel, then diluted to 1 ng/μL. The 16S rRNA V4 regions were amplified with barcoded primers using reactions containing Phusion® Master Mix (15 μL), 0.2 μM primers, and 10 ng DNA. The primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used (Walters et al., 2016; Snijders et al., 2016). PCR conditions included an initial denaturation at 98 °C for 1 min, followed by 30 cycles at 98 °C for 10 s, 50 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. PCR products were visualized on 2% agarose gels stained with SYBR Green, pooled equally, and purified using a Qiagen Gel Extraction Kit. Libraries were constructed using the NEBNext Ultra DNA Library Prep Kit, quality-verified by an Agilent Bioanalyzer, and sequenced on an Illumina NovaSeq platform with 250 bp paired-end reads.

Barcodes are used to separate and obtain the original data for each sample, followed by the removal of barcodes and primers. R1 and R2 sequence data are concatenated using FLASH software. Quality control is performed on concatenated tags to obtain clean tags, followed by chimeric filtering to yield effective tags for analysis. To ensure accurate and reliable outcomes, raw data undergo concatenation and filtering to produce clean data. Based on clean data, denoising is performed using DADA2 (Callahan et al., 2016), and sequences with an abundance of less than five are filtered out to obtain the final ASVs. Species annotation for each ASV is conducted using the classify-sklearn algorithm in QIIME2 (Bokulich et al., 2018; Bolyen et al., 2019) with a pre-trained Naive Bayes classifier (SILVA 138.1).

To analyze the diversity, richness, and uniformity of the communities in the sample, alpha diversity was calculated using three indices in QIIME2: Chao1, Shannon, and Simpson. Beta diversity, based on weighted UniFrac distances in QIIME2, was calculated to evaluate the complexity of community composition and compare differences between groups. Two-dimensional PCoA results were visualized using the ade4 and ggplot2 packages in R software (Version 2.15.3). MetaStat and T-test analyses were performed using R software (Version 3.5.3) to identify significantly different species at each taxonomic level (Phylum, Class). LEfSe software (Version 1.0) was employed to perform LEfSe analysis (LDA score threshold: 4) for identifying biomarkers. The functional differences among different groups were tested using the T-test method.

A 100 μL aliquot of each sample was transferred to a sterile and enzyme-free EP tube, with 400 μL of extraction solution (acetonitrile:methanol = 1:1, containing isotopically labeled internal standards) added. The mixture was vortexed for 30 s, sonicated in an ice-water bath for 10 min, and incubated at −40 °C for 1 h to precipitate proteins. Samples were centrifuged at 12,000 rpm (RCF = 13,800 × g, R = 8.6 cm) at 4 °C for 15 min, and the supernatant was transferred to a fresh glass vial for analysis.

LC–MS/MS analyses were performed using a Thermo Fisher Scientific Vanquish UHPLC system equipped with a UPLC BEH Amide column (2.1 mm × 100 mm, 1.7 μm), coupled to a Q Exactive HFX Orbitrap mass spectrometer. The mobile phase included 25 mmol/L ammonium acetate and 25 mmol/L ammonia hydroxide in water (pH = 9.75) as phase A and acetonitrile as phase B. The autosampler temperature was kept at 4 °C, with an injection volume of 2 μL. The mass spectrometer operated in information-dependent acquisition (IDA) mode, controlled by Xcalibur software, continuously evaluating full-scan MS spectra. ESI source parameters included: sheath gas flow rate at 30 Arb, auxiliary gas flow rate at 25 Arb, capillary temperature at 350 °C, full MS resolution of 60,000, MS/MS resolution of 7,500, collision energy of 10/30/60 in NCE mode, and spray voltage of 3.6 kV (positive ion mode) or −3.2 kV (negative ion mode). Raw data were converted to mzXML format using ProteoWizard and processed with an in-house program developed in R and based on XCMS for peak detection, extraction, alignment, and integration. Metabolite annotation was conducted using an in-house MS2 database, BiotreeDB version 2.1, which is proprietary and not publicly accessible. This comprehensive database consists of a self-constructed library containing 3,000 standard compounds, an integrated public library with 200,000 substances, and a lipid database comprising 190,000 entries. The confidence level of metabolite identification was determined based on score values from the mean table, using a cutoff threshold of 0.3.

The data were subjected to logarithmic (LOG) transformation and UV scaling using SIMCA software (version 16.0.2, Sartorius Stedim Data Analytics AB, Umea, Sweden). OPLS-DA modeling analysis was conducted on the first principal component. A permutation test was performed by randomly rearranging the order of the categorical variable Y multiple times to obtain different random Q values, validating the model’s validity.

Screening for differential metabolites primarily involves three parameters: VIP, FC, and p-value. VIP (Variable Importance in Projection) indicates the importance of a variable in the projection of the first principal component of the PLS-DA model, showing each metabolite’s contribution to group differentiation. FC (Fold Change) is the ratio of the average quantitative values of each metabolite in the comparison group. Significance thresholds are set at VIP > 1.0 and p-value < 0.05.

KEGG analysis: Pathway enrichment analysis was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database¹ to interpret the biological functions of identified differential metabolites. All pathways corresponding to differential metabolites in Mus musculus (mouse) were identified, and these metabolites were marked on the KEGG pathway map.

K-Means analysis: To study trends in metabolite content changes across different groups, the relative content of all differentially expressed metabolites was z-score standardized according to the screening criteria across all group comparisons, followed by K-Means clustering analysis.

Spearman correlation analysis calculated the correlation coefficient between changes in gut microbiota and metabolites. Inter-group differences in intestinal microbiota were identified based on LEfSe analysis results. Differential metabolites were identified by combining the metabolomic analysis results of JFAD, which can regulate the health of mice, with those of lung cancer-bearing mice. Spearman correlation analysis evaluated potential associations between inter-group differences in intestinal microbiota and differential metabolites, assessing the existence and strength of the correlations. This correlation is considered the most likely target for JFAD’s anti-tumor effect. Spearman’s correlation coefficient, ranging from −1 to 1, described the relationship between these parameters.

We incorporated a supervised multi-omic factor analysis using the DIABLO implementation of block sPLS-DA from the mixOmics package. This model jointly extracts latent components from both the microbiome and metabolome datasets, allowing us to identify features that co-vary across omics layers and contribute directly to the separation of the NG, MG, JFAD-L, JFAD-M, and JFAD-H groups. To further illustrate cross-omic interactions, we added a DIABLO-based Circos plot, which visualizes all significant correlations between the selected microbial features and metabolites identified by the factor model. This provides an intuitive global view of microbe-metabolite connectivity and highlights several key cross-omic pairs that appear to be modulated by JFAD treatment. In addition, we constructed Clustered Image Maps (CIM) based on DIABLO latent components to map block-to-block correlation structures. We also complemented the factor-based analysis with a high-resolution correlation matrix focusing on the top 30 microbes and top 30 metabolites showing the strongest cross-omic associations. This targeted heatmap provides a refined view of microbe–metabolite pairs most relevant to phenotype discrimination and supports the interaction patterns identified by DIABLO. Additionally, to deepen the mechanistic interpretation, we built an integrated correlation network summarizing all high-confidence microbe-metabolite associations. This network highlights several taxa and metabolites that function as central nodes, providing further evidence for coordinated functional pathways potentially involved in the therapeutic effects of JFAD.

Significance levels were determined using appropriate statistical analyses according to data distribution and the number of test groups. GraphPad Prism version 10.1.2 was used for analysis if the data followed a normal distribution. Pairwise comparisons were conducted using one-way ANOVA when the variance was homogeneous.

In this research (Peng et al., 2024), tumor weight decreased with JFAD treatment, with tumor inhibition rates of 15.24% (low), 35.65% (medium) and 7.57% (high). The medium dose produced the greatest and statistically significant tumor inhibition (p < 0.05), indicating that medium-dose JFAD effectively suppresses lung adenocarcinoma in this model. The bioluminescence imaging data was also presented. On a molecular level, JFAD was found to attenuate phosphorylation within the PI3K/Akt pathway, indicating one possible mechanism for its antitumor effects. This earlier work suggests that JFAD, particularly at medium doses, may exert antitumor activity by modulating crucial signaling pathways, warranting further exploration into its therapeutic potential (Figure 1).

To evaluate the regulatory effects of JFAD on the intestinal microbiome, an alpha diversity assessment was initially conducted, revealing no significant differences in the alpha diversity index between the MG group and other groups (Figures 2A–C). Subsequently, beta diversity analysis was performed to examine community structure variations, with PCoA analysis indicating that the NG group was more distant from the MG group, suggesting a more pronounced disparity in species composition. This suggests that the intestinal microbiome of tumor-bearing mice underwent notable changes. Following JFAD intervention, each group (JFAD-L, JFAD-M, JFAD-H) was positioned between the MG and NG groups, with intestinal microbiota structures closely resembling that of the NG group (Figure 2D). This implies that JFAD might assist in restoring the disrupted intestinal flora in lung cancer-afflicted mice to a healthier state. Similarly, NMDS analysis revealed that the species composition of the NG and JFAD groups was more comparable to the MG group, consistent with PCoA results (Figure 2E). Additionally, the Venn diagram shows the number of shared and unique ASVs among the groups, with 318 shared ASVs: 360 unique to the NG group, 647 to the MG group, 952 to the JFAD-L group, 392 to the JFAD-M group, and 531 to the JFAD-H group (Figure 2F). These results suggest that the sequencing depth for each classification unit is sufficient.

The taxonomic composition at the phylum and class levels was analyzed. At the phylum level, Firmicutes and Bacteroidota were significantly dominant across all groups (Figure 3A). Compared to the NG group, the MG group showed a higher abundance of Firmicutes and a lower abundance of Bacteroidota. The JFAD intervention resulted in a reduction of Firmicutes and an increase in Bacteroidota. The Firmicutes to Bacteroidota (F/B) ratio was 1.732 (SD: 0.9933) in the NG group, 10.73 (SD: 13.85) in the MG group, 5.78 (SD: 11.84) in the JFAD-L group, 0.8822 (SD: 0.2861) in the JFAD-M group, and 2.23 (SD: 1.599) in the JFAD-H group. At the class level, Clostridia and Bacteroidia were clearly dominant across all groups (Figure 3B). Compared to the NG group, the MG group showed a higher abundance of Clostridia and a lower abundance of Bacteroidia. The JFAD intervention resulted in a reduction of Clostridia and an increase in Bacteroidia. Therefore, JFAD may enhance the abundance of Bacteroidia while reducing Firmicutes and Clostridia, thus altering the composition of intestinal bacteria.

To identify key phylotypes and biomarkers of gut microbiota across different groups, a comparative analysis using LEfSe was performed. This analysis revealed significant differences in bacterial community dominance between groups. In the LEfSe analysis (Figures 3C,D), the bar chart length represents the LDA score, indicating the extent to which the different species affect group differences. LDA scores greater than 4 were considered statistically significant. Bacteroidaceae and Clostridia exhibited higher LDA scores, suggesting these bacterial taxa can serve as biomarkers, with notable differences observed between groups.

By comparing the MG with NG group, we revealed the significant impact of herbal JFAD on the metabolic functions of microbial communities. Compared to the NG group, the MG group showed abnormal expression in multiple metabolic pathways (Figure 4A). Specifically, the arginine metabolism pathway marker K07052 and the amino sugar metabolism enzyme K00604 were upregulated in the MG group, while pathways involved in core nitrogen metabolism, including K00278, K00334, K12267, K03340, and K00335, were downregulated. Following JFAD-L intervention, the elevated levels of K07052 and K00604 in the MG group were suppressed. In contrast, the nitrogen metabolism module K03340 and phosphonate metabolism-related markers K01546, K01547, and K01548 were significantly upregulated (Figure 4B). JFAD-M intervention further enhanced nitrogen and phosphonate metabolism functions, including the upregulation of K00278, K00334, and K12267, as well as stress-adaptive transport and regulatory factors K18331 and K18332. This suggests that medium-dose intervention may have a more profound impact on improving microbiota function (Figure 4C). JFAD-H, in the above metabolic pathways, showed no significant differences in regulation (Figure 4D).

The OPLS-DA model was used to assess the extent of separation among various groups. The OPLS-DA score plot illustrates that the serum metabolome of the NG group is clearly distinguishable from that of the MG group (Figure 5A). This indicates that there are differences in the metabolic profiles between the two groups. The scattered points are relatively concentrated, indicating strong repeatability within each group, with similar metabolite profiles, thus demonstrating the reliability of each group’s model. Permutation test results for the OPLS-DA model, indicated by R²Y and Q² values, suggest that the model does not suffer from overfitting, as the blue regression line crosses the left ordinate below the zero mark (Figure 5B). Furthermore, an analysis of metabolic pathways involving differing metabolites between the NG and MG groups showed that these metabolites were mainly associated with the metabolism of nicotinate and nicotinamide, alpha-Linolenic acid metabolism, as well as the synthesis of phenylalanine, tyrosine, and tryptophan (Figure 5C). The combined analysis of the two ion patterns yielded results that were nearly identical to the current results.

The first principal component (PC1) accounts for the largest variance, providing the most explanation for differences in the original data. In comparison, the second and third principal components are less significant. In positive ion mode, the scores of the first principal component are relatively similar, resulting in samples from each group clustering together. Closer sample points indicate more similar types and contents of metabolites. In negative ion mode, the detection results reveal that post-JFAD intervention, the distribution points align closely with those of the NG group, clustered around −20 on the PC1 axis. Conversely, the MG group is positioned nearer to 0 on the PC1 axis. This suggests that cancer may disrupt the metabolic pattern of mice, while JFAD intervention appears to regulate the metabolic pattern in lung cancer-bearing mice (Figure 6A).

In positive ion mode, identified metabolites were categorized by their chemical classifications, showing that a substantial portion consisted of lipids and lipid-like molecules (28.111%), organoheterocyclic compounds (23.272%), and organic acids and derivatives (20.507%) (Figure 6B). In negative ion mode, the metabolites were similarly categorized, with lipids and lipid-like molecules (36.889%), organic acids and derivatives (19.556%), and organoheterocyclic compounds (14.222%) being predominant. In both modes, the top three categories with the highest proportions of identified metabolites were lipids and lipid-like molecules, organic acids and derivatives, and organoheterocyclic compounds, with slight differences in proportions. Positive and negative ion patterns were combined for analysis. Using KEGG pathway annotation, comparisons were made between the MG group and the NG group. It was observed that metabolites related to metabolic processes such as pyrimidine metabolism (17.24%), nicotinate and nicotinamide metabolism (13.79%), glycine, serine, and threonine metabolism (6.9%), and alpha-linolenic acid metabolism (6.9%) had relatively high occurrence rates. Following JFAD intervention, the proportions of metabolites associated with these metabolic pathways decreased (Figure 6C).

Differential metabolites associated with lung cancer were identified through a screening process. Using data from the NG and MG groups, metabolites with a VIP (Variable Importance in Projection) score above 1 and a t-test p-value below 0.05 were selected. Sixty-seven differential metabolites were identified in POS mode, with significant differences noted after the JFAD intervention (Table 2). Metabolites related to nicotinate and nicotinamide metabolism included 1-methylnicotinamide, niacinamide, and N1-methyl-4-pyridone-3-carboxamide. Compared to the NG group, the relative levels of these metabolites were significantly lower in the MG group. Significant increases in the levels of 1-methylnicotinamide and niacinamide were observed in the JFAD-M and JFAD-H groups. Additionally, differential metabolites associated with glycine, serine, and threonine metabolism included betaine aldehyde and guanidoacetic acid. The levels of these metabolites were significantly higher in the MG group compared to the NG group, but JFAD treatment reduced their levels across all doses. Using the same screening technique, 34 distinct metabolites were identified in NEG mode in the NG and MG groups (Table 3). Notably, succinate, a metabolite associated with pyrimidine metabolism, increased in concentration in the MG group relative to the NG group, but JFAD treatment successfully reduced its levels. The Euclidean distance matrix for the quantitative values of differential metabolites across each group was computed, and the complete linkage approach was used to cluster these metabolites.

K-Means analysis was applied to all significant differential metabolites. In positive ion mode, a total of 301 significant differential metabolites were divided into 9 clusters, spanning from the NG group to the MG group and each JFAD group (Figure 7A). Notably, compared to the MG group, metabolite concentrations in clusters 1 and 8 showed a reversal following each JFAD intervention, suggesting these 81 metabolites may relate to the anti-tumor effect mediated by JFAD. In negative ion mode, 160 significant differential metabolites were divided into 9 clusters, spanning from the NG group to the MG group and each JFAD group (Figure 7B). Compared to the MG group, metabolite concentrations in clusters 8 and 9 showed a reversal after each JFAD intervention, indicating that these 47 metabolites (clusters 8 + 9) may also relate to the anti-tumor effect mediated by JFAD. Specific data are provided in Supplementary Tables S1, S2.

Upon acquiring data for various groups of differentially expressed metabolites, the KEGG database specific to Mus musculus (mouse) was used to perform pathway searches and examine regulatory interaction networks, as shown in Figure 8A. This analysis identified 15 pathways, 5 modules, 54 enzymes, 111 reactions, and 35 compounds. It reveals intersections between the metabolic pathways of the MG and NG groups, highlighting potential targeted enzymes and metabolites. To elucidate the potential anti-tumor effects of JFAD through regulation of metabolites, the focus was directed toward the central carbon metabolism pathway (mmu05230) in cancer, as shown in the complex metabolic network in Figure 8A. The KEGG pathway diagram indicates a notable increase in succinate levels and a significant decrease in arginine levels in the MG group compared to the NG group (Figures 8B,C). Compared to the MG group, all doses of JFAD treatment (low, middle, and high) led to down-regulation of succinate. Additionally, essential amino acids were significantly down-regulated in the MG group relative to the NG group. After treatment with JFAD-L and JFAD-H, no notable changes were observed in the differential metabolites. A significant increase in essential amino acids was observed in the JFAD-M group. This observation may be associated with the PI3K/Akt and mTOR signaling pathways.

Spearman correlation analysis was performed to explore relationships between six significantly different families and three differentially altered phyla with differential metabolites in both positive and negative ion modes (Figure 9). Differential metabolites in positive and negative ion modes are shown in Tables 2, 3. The results showed that 33 differential metabolites exhibited a strong correlation with these nine bacterial groups in positive ion mode. In negative ion mode, 15 differential metabolites exhibited a strong correlation with the same nine bacterial groups. These findings suggest that in vivo metabolite levels can reflect the gut microbiota structure. Additionally, betaine aldehyde and guanidinoacetic acid, related to glycine, serine, and threonine metabolism pathways, were negatively correlated with Tannerellaceae and Campylobacterota. Succinic acid, a differential metabolite in the central carbon metabolism pathway of cancer, was also negatively correlated with the gut microbiota families Tannerellaceae and Campylobacterota.

To further clarify the coordinated changes between gut microbiota and serum metabolites, we performed an integrated sPLS-DA analysis. The score plots showed a clear separation of the NG, MG, and JFAD groups in both the microbe and metabolite blocks (Figure 10A). MG samples clustered apart from NG, indicating a marked disruption of microbe-metabolite structure in the disease state. JFAD-treated mice shifted closer to NG along the first component, suggesting a partial restoration of the joint microbial-metabolic profile after treatment. The microbe-metabolite clustered heatmap revealed several consistent correlation patterns (Figure 10B). Taxa such as Prevotella, Anaerotruncus, Blautia, and Oscillibacter showed strong positive associations with metabolites including 2-methyl-5-propylpyrazine, and formyl-methionyl-serine. In contrast, Methanosaeta, Methanomassiliicoccus, and Staphylococcus were mainly associated with lower levels of these metabolites. These opposing clusters matched the differences observed between NG and MG mice and aligned with the reversal trend seen in the JFAD groups.

The loading plots further identified variables contributing most to the group separation (Figures 11A–D). Staphylococcus, Thalassospira, Methanolinea, and IS1 showed the highest loading weights in the microbe block on component 1, whereas 2-pyridylacetic acid, 3-methyl-L-tyrosine, inosine, and bilobalide were the main contributors in the metabolite block. Component 2 highlighted additional taxa and metabolites that distinguished the different JFAD doses. These findings indicate that specific combinations of bacteria and metabolites, rather than isolated changes, are responsible for the metabolic imbalance in MG mice and the modulatory effects of JFAD.

The circos plot revealed a dense network of strong correlations (|r| ≥ 0.7) between gut microbes and serum metabolites. Although most connections were negative, these relationships still represented close associations between specific bacterial taxa and metabolite profiles. A large cluster of strong negative correlations was observed linking taxa such as Methanobrevibacter, Acinetobacter, Staphylococcus, Erysipelatoclostridiaceae, Prevotella, and Trichococcus with multiple differential metabolites involved in amino acid metabolism, steroid derivatives, and lipid-related pathways. These patterns indicate that changes in microbial composition were tightly mirrored by coordinated shifts in metabolite abundance across groups. Positive correlations were less frequent, but one notable strong positive association was identified between Methanosaeta and estriol, suggesting a distinct microbe-metabolite pair that varied in the same direction. The circos network illustrates that gut microbiota and metabolites form a highly interconnected system, in which both positive and negative associations contribute to the integrated metabolic response underlying the group differences (Figure 12).