Clinically distinct metabotypes of pediatric MASLD identified through unsupervised clustering of NASH CRN data

Domain: chemistry
Published: 2026-02-24
Editorial outcome: major_revision

Abstract
Metabolic dysfunction-associated steatotic liver disease (MASLD) is the most common liver disease worldwide, yet treatment remains “one size fits all,” despite phenotypic heterogeneity. We analyzed clinical and metabolomics data from 514 children (ages 5-18, 73% male) with biopsy-proven MASLD across three NASH Clinical Research Network studies. Unsupervised clustering of clinical data identified three distinct metabotypes: early-mild (49.4%, youngest, lowest lipids, liver enzymes, insulin resistance), cardiometabolic (36.8%, highest waist circumference, lipids, uric acid, SBP), and inflammatory-fibrotic (13.8%, highest liver enzymes, steatohepatitis, advanced fibrosis). Integrative network and pathway enrichment analyses revealed alterations in tryptophan metabolism within the inflammatory-fibrotic group, including elevated kynurenine pathway metabolites, which were significantly correlated with fibrosis stage. Branched-chain amino acid degradation, butanoate, and purine metabolism demonstrated greater enrichment in the cardiometabolic group. Here, we show that pediatric MASLD subtypes differ in clinical and metabolic features, providing a framework for targeted interventions, with validation needed in independent cohorts.

Introduction
Metabolic dysfunction-associated steatotic liver disease (MASLD), previously known as nonalcoholic fatty liver disease (NAFLD), is the most prevalent liver disease in children and is a leading contributor to liver-related morbidity and mortality1,2. Once considered rare, the prevalence of pediatric MASLD has significantly increased and now affects an estimated 16.5% of U.S. adolescents3, including 26% of those with obesity4. MASLD is characterized by the presence of hepatic steatosis ≥5%, along with at least one out of five specified cardiometabolic risk factors, in the absence of heavy alcohol use and other chronic liver diseases1. In children and adolescents, these cardiometabolic risk factors include increased body mass index (BMI), fasting glucose, blood pressure, triglycerides (TG), and low high-density lipoprotein cholesterol (HDL). MASLD severity exists on a spectrum that ranges from simple steatosis to metabolic dysfunction-associated steatohepatitis (MASH), formerly known as nonalcoholic steatohepatitis (NASH), which can progress to fibrosis and cirrhosis, predisposing one to hepatocellular carcinoma1. The current and future health burden of pediatric MASLD is immense, both for affected individuals and from a public health perspective. The disease presents a significant clinical challenge as it is associated with additional cardiovascular, metabolic, and psychological disorders5,6. Patient quality of life is also adversely affected, and there is an enormous economic burden due to healthcare utilization, which is likely to increase in adulthood7. Despite this burden, there are currently no FDA-approved medications or supplements for treating pediatric MASLD, and lifestyle modifications have demonstrated limited sustained success in reducing hepatic steatosis and serum alanine aminotransferase (ALT) levels, a key biomarker of liver injury and disease progression in MASLD8–12.

Recent studies in adults have highlighted that MASLD is a heterogeneous condition with a broad spectrum of clinical presentations, pathogenesis, and disease severity13–17. This heterogeneity stems from various factors, including age, sex, hormonal status, genetics, gut microbiota composition, comorbidities, and certain exposures, including diet and physical activity14. Moreover, there is substantial inter-individual variation in response to current therapies for MASLD, which complicates treatment approaches13,18,19. Recent research in adults has identified sub-phenotypes of MASLD using clinical, genetic, and omics data20–27. These studies have uncovered differences in MASLD development and progression linked to genetic variants, lipid partitioning, body composition, and inflammation28. Given the heterogeneity observed in adults, it is plausible that sub-phenotypes also exist in children with MASLD. In line with the global shift towards precision medicine and public health, identifying MASLD phenotypes in youth could guide personalized interventions to prevent the onset and progression of hepatic steatosis. Unsupervised clustering analyses have been successfully used to identify subtypes in various diseases29–31. In this work, we identify three clinically distinct metabolic phenotypes (metabotypes) of pediatric MASLD, each characterized by unique metabolic and histological profiles, using unsupervised clustering and high-resolution metabolomics (HRM).

Results
NASH Clinical Research Network (CRN) cohort characteristics
Among the 514 children with biopsy-confirmed MASLD, ages ranged from 5 to 18 years, with a median age of 12 years. The cohort was predominantly male (73%) and of Hispanic ethnicity (74%), with 69% diagnosed with borderline or definite MASH. Notable differences were observed across the three NASH CRN study cohorts. The proportion of male participants included in our analysis was highest from the Treatment of NAFLD in Children (TONIC) study (85%) compared to the NAFLD Pediatric Database (DB1) (57%) and the NAFLD Pediatric Database 2 (DB2) (p=0.006). The percentage of participants of Hispanic ethnicity also varied, with 77% in DB2, 65% in DB1, and 62% in TONIC (p=0.009). Significant differences in liver enzyme levels were observed, with ALT and aspartate aminotransferase (AST) being highest in the TONIC and DB2 cohorts (p<0.001). Glycated hemoglobin (HbA1c), BMI z-score, and NAFLD activity scores (NAS) also varied significantly across cohorts (p<0.05). Please see Supplementary Data 1 for additional details.

Clinically distinct metabotypes of pediatric MASLD
Our unsupervised clustering analysis revealed three clinically distinct metabotype clusters of pediatric MASLD (Table 1; Fig. 1). Based on the underlying metabolic phenotype, the clusters were categorized as follows: 1) early-mild (49.4% of participants), 2) cardiometabolic (36.8%), and 3) inflammatory-fibrotic (13.8%). The early-mild metabotype was significantly younger with lower levels of lipids (very-low-density lipoprotein cholesterol [VLDL], TG, low-density lipoprotein cholesterol [LDL]), liver enzymes (ALT, AST), and insulin resistance (Homeostasis Model Assessment of Insulin Resistance [HOMA2-IR]) (all p<0.001). The cardiometabolic metabotype presented with the highest waist circumference (WC), VLDL, TG, and uric acid levels, with significantly elevated systolic blood pressure (SBP) (p<0.001). Lastly, the inflammatory-fibrotic metabotype exhibited significantly higher levels of liver enzymes and had a greater percentage of participants with steatohepatitis and advanced fibrosis (p<0.001). Interestingly, the early-mild and inflammatory-fibrotic metabotypes had higher percentages of Hispanic participants (81% and 82%, respectively) than the cardiometabolic metabotype (61%; p<0.001). Furthermore, the early-mild metabotype had a higher percentage of participants with advanced fibrosis than cardiometabolic (17% vs. 8%, respectively; p=0.009). Fig. S1 shows that the silhouette widths for the observed data consistently exceed those from the permuted datasets, indicating that the observed clusters are stronger and more meaningful than would be expected by chance. Additionally, as illustrated in Fig. S1, cluster differences were not primarily driven by sex, as male and female participants were represented across all clusters with no clear separation along the first two principal components.

Metabolomic profiling of pediatric MASLD metabotypes
As described above, the primary goal of our HRM analysis was to examine differences in the metabolome across the three metabotypes. The hydrophilic interaction liquid chromatography positive mode (HILIC+) and C18 negative mode (C18-) principal component analysis (PCA) plots showed minimal separation among the early-mild, inflammatory-fibrotic, and cardiometabolic groups. However, the partial least squares discriminant analysis (PLS-DA) plots revealed clearer separation, with cardiometabolic displaying more distinct clustering, while early-mild and inflammatory-fibrotic still exhibited some overlap, suggesting the presence of unique metabolic profiles among the three metabotypes (Fig. S3 and S4). One-way analysis of variance (ANOVA) identified 248 significant features from HILIC+ mode and 178 from C18- mode (false discovery rate [FDR] p<0.05: Supplementary Data 2 & 3). Among these, the most significant HILIC+ annotations included uric acid, which was highest in the cardiometabolic metabotype, as well as a phosphatidylcholine species (PC36:2) and a sphingolipid species (SM34:1), both elevated in the early-mild metabotype (all p<0.05). In C18- mode, the top annotated metabolites included uric acid and related compounds involved in uric acid metabolism, such as 5-hydroxyisourate and urate, along with leucine/isoleucine. Uric acid, 5-hydroxyisourate, and urate levels were highest in the cardiometabolic metabotype (p<0.001). Leucine/isoleucine was also significantly higher in the cardiometabolic group compared to the early-mild group (p<0.0001) and showed a trend toward statistical significance compared to the inflammatory-fibrotic group (p=0.065) (Fig. S3C and S4C).

Pathway enrichment analysis of the significant HILIC+ and C18- features identified by one-way ANOVA (Fig. 2; and Supplementary Data 4) revealed six metabolic pathways that were significantly altered across the three metabotypes. These pathways included tryptophan metabolism (p<0.001), branched-chain amino acid (BCAA) degradation (p=0.001), butanoate metabolism (p=0.004), propanoate metabolism (p=0.008), pantothenate and CoA biosynthesis (p=0.02), and purine metabolism (p=0.02). Group-wise distributions of the enriched features in each pathway are illustrated in Fig. S5-S10 using heatmaps and violin plots. Notably, several metabolites within the tryptophan metabolism pathway, including tryptophan, serotonin, kynurenine, and indole-3-acetaldehyde, as well as pantothenate and CoA biosynthesis (aspartate and serine), were elevated in the inflammatory-fibrotic metabotype compared to cardiometabolic and early-mild (all p<0.05; Fig. S5 and S6). In contrast, purine metabolism, BCAA degradation, and butanoate metabolism showed higher levels of metabolites in the cardiometabolic metabotype relative to the other two groups, while several metabolites of propanoate metabolism were higher in the early-mild metabotype (p<0.05; Fig. S7–S10).

Among the significantly enriched pathways, several contained at least one Schymanski Level Confidence (SLC) Level 1 identification from our internal reference library, supporting the confidence of these results32. Three pathways, propanoate metabolism, BCAA degradation, and pantothenate and CoA biosynthesis, contained features that could not be independently confirmed by secondary annotation methods and were therefore reported as SLC Level 5. These results should be interpreted cautiously but may represent true biological signals and provide important hypotheses for future validation.

Network analysis of clinical & metabolomic features across pediatric MASLD metabotypes
Figure 3 shows the integrative networks of the clinical and metabolic features for the three metabotypes. VLDL and TG had the highest centrality scores across all networks, underscoring their prominent roles within each group (Supplementary Data 5). The inflammatory-fibrotic network also showed high centrality for clinical markers associated with inflammation and metabolic dysfunction, such as ALT, AST, fibrosis stage, and HOMA2-IR, reflecting their central roles within this network. Similarly, the early-mild network included clinical variables such as ALT, AST, fibrosis stage, and HOMA2-IR within the same community, suggesting related pathways or functions among these variables in both networks. In contrast, the cardiometabolic network was predominantly lipid-focused, with LDL and WC positioned alongside VLDL and TG, highlighting a central orientation around lipid and adiposity markers. These orientations and centrality differences across the networks underscore the distinct metabolic heterogeneity characterizing each of the three metabotypes.

Comparing the early-mild and cardiometabolic metabotype networks, we identified three significantly enriched pathways based on metabolic features differing in centrality: tryptophan metabolism, pantothenate and CoA biosynthesis, and BCAA degradation (p<0.05; Fig. 3; and Supplementary Data 6). Pairwise comparisons showed that the majority of these metabolites, including 3-methyl-2-oxobutanoic acid, valine, and 4-methyl-2-oxopentanoate, were significantly elevated in the cardiometabolic metabotype versus early-mild (p<0.05; Fig. S11). Similar pathway alterations were observed in the comparison between the early-mild and inflammatory-fibrotic metabotype networks; however, glycine, serine, and threonine metabolism was also significantly enriched, with features in these pathways such as tryptophan, kynurenine, aspartate, and hydroxyproline, notably higher in the inflammatory-fibrotic metabotype (p<0.05; Fig. 3 and S12). Lastly, when comparing the inflammatory-fibrotic and cardiometabolic networks, tryptophan metabolism, BCAA degradation, pantothenate, and CoA biosynthesis, as well as glycine, serine, and threonine metabolism, were significantly enriched. Metabolic features from tryptophan and glycine, serine, and threonine metabolism were significantly higher in the inflammatory-fibrotic metabotype (p<0.05), while those from BCAA degradation, including leucine/isoleucine, were elevated in the cardiometabolic metabotype (Fig. 3 and S13). Given the consistent alteration of metabolites in the tryptophan pathway across metabotypes, we further examined their associations with clinical variables in the network analysis. Notably, tryptophan, indole-3-acetaldehyde, serotonin, 5-methoxyindoleacetate, 5-phenyl-1,3-oxazinane-2,4-dione, and kynurenine showed moderate to strong positive correlations with ALT, AST, and fibrosis stage (Supplementary Data 7).

Correlations between metabolic features and fibrosis stage
Among the 3,758 mass-to-charge (m/z) features from HILIC+ mode and 3520 from C18- mode, 351 and 249 features, respectively, showed significant associations with the fibrosis stage (p<0.05; Supplementary Data 8). As shown in Fig. 4, the top metabolites correlated with fibrosis stage were 4-hydroxy-5-methyl-3(2H)-thiophenone, a sulfur-containing thiophene derivative potentially involved in oxidative stress pathways33,34 (r=0.23, p=1.42 × 10⁻⁷; HILIC+ mode), and an unknown feature (437.1615_66.9; r=−0.23, p=3.03 × 10⁻⁷; C18- mode). Additionally, several tryptophan-related metabolites, such as indole, tryptophan, and kynurenine, were positively associated with fibrosis, suggesting a link between tryptophan metabolism and fibrosis severity. Lipid-related metabolites, including phosphatidylethanolamine (PE), phosphatidylcholine (PC), and phosphatidic acid (PA) species, also showed positive correlations with fibrosis stage. Furthermore, polycyclic aromatic hydrocarbons (PAHs), such as pyrene and 2-aminonaphthalene, as well as nucleotide metabolism-related features, including deoxyinosine and xanthine, and bile acid-related metabolites such as chenodeoxycholic acid derivatives, were also positively correlated with fibrosis.

Methods
All participants included in this study provided written informed consent and assent for their samples and data to be used for future research. Written informed consent was obtained from a parent or legal guardian for all participants, regardless of age. Participants aged 13-17 years provided written adolescent assent, and those under 13 years provided written child assent. Research was conducted in accordance with both the Declarations of Helsinki and Istanbul, and the study protocol was approved by the institutional review boards at Emory University and Children’s Healthcare of Atlanta (Atlanta, GA; STUDY00001715).

This retrospective cross-sectional study employed a two-step framework (Fig. 6). First, a data-driven clustering analysis was performed using anthropometric and clinical data from youth enrolled in NASHCRN studies to identify metabotypes of pediatric MASLD. Second, an exploratory high-resolution metabolomics analysis of participant serum samples was conducted to examine differences in the metabolome among the identified metabotypes.

NASH CRN cohort
The sample population included a combined cohort of children and adolescents aged 5–18 years with liver biopsy-proven MASLD enrolled in one of three National Institute for Diabetes and Digestive and Kidney Diseases (NIDDK)-funded NASH CRN studies: the Treatment of NAFLD in Children (TONIC) trial (n=93) [NCT00063635], the NAFLD Pediatric Database (DB1) study (n=23), and the NAFLD Pediatric Database 2 (DB2) study (n=398) [NCT01061684].

The NAFLD Pediatric Database studies (DB1 and DB2)64 were observational analyses of patients two years and older with MASLD, and TONIC was a phase III, masked, randomized, placebo-controlled trial of metformin or vitamin E in children ages 8–17 years with MASLD65. Inclusion and exclusion criteria were similar among the three studies64,66. Briefly, eligible participants were ≤18 years of age with MASLD confirmed by liver biopsy. Exclusion criteria included any clinical or histologic evidence of alcoholic liver disease, evidence of other causes of chronic liver disease, greater than one month of total parenteral nutrition six months before the liver biopsy, short bowel syndrome, biliopancreatic diversion, bariatric surgery, and known HIV infection. Participants in the TONIC study were also required to have a baseline ALT value of ≥60 U/L, and additional exclusion criteria were the diagnosis of diabetes, metabolic acidosis or renal dysfunction, and the use of anti-diabetic or anti-MASLD drugs.

For the present study, children were included if they met the current diagnostic criteria for MASLD confirmed by liver biopsy (hepatic steatosis ≥5%) along with at least one marker of cardiometabolic dysfunction1 and had serum samples processed for high-resolution metabolomics at the Emory Clinical Biomarkers Laboratory. The initial eligible sample included 556 children and adolescents with available clinical and metabolomics data. A small percentage of missing data (0.34%) was imputed using the median for clinical variables. Outliers, defined as participants with any clinical variable greater than five standard deviations from the mean, were excluded, resulting in 515 participants. During this process, most participants with an HbA1c ≥6.5% were excluded; however, one participant remained in the dataset. To maintain consistency and avoid bias from a single value, this participant was also excluded, yielding a final sample of 514 children and adolescents for analysis (Fig. S14).

Clinical and laboratory assessments
Demographic data for the three studies were collected through structured interviews and questionnaires. Height, weight, and waist circumference measurements were measured in duplicate, with participants standing in light clothing without shoes. Body mass index (BMI) was calculated by dividing weight (kg) by height (m) squared. BMI percentiles were determined based on age and sex using the Centers for Disease Control and Prevention (CDC) criteria67.

Fasting blood samples were collected from all participants after an overnight fast (minimum of 8 h or more) and processed for plasma and serum within 2 h64,65. Laboratory assays included in the present study were ALT, AST, alkaline phosphatase, gamma-glutamyl transferase (GGT), platelet count, fasting glucose, insulin, HbA1c, bilirubin, and uric acid levels. Additionally, a lipid panel including fasting TG, total cholesterol, LDL, and HDL was included. VLDL was calculated by dividing triglyceride levels by five, consistent with established methods68. The updated homeostasis model assessment of insulin resistance (HOMA2-IR) was calculated using the online HOMA calculator version 2.2.369. Systolic and diastolic blood pressure were assessed using standard procedures70. Additional fasting samples from all study participants were stored at −80 °C at the NIDDK Central repository for future use. Please see our Supplementary Materials for a table summarizing participant baseline clinical and demographic characteristics stratified by study (Supplementary Data 1).

Histological evaluation
For each study, biopsies were centrally assessed by the NASH CRN Pathology Committee as a group for consensus scoring according to the system by Kleiner et al71. Assessments included the NAFLD activity score (NAS) on a scale of 0-8, with higher scores indicating more severe disease; the components of this measure include steatosis (0-3), lobular inflammation (0–3), and hepatocellular ballooning (0-2). Biopsies were also evaluated for the presence or absence of MASH and the stage of fibrosis, which was classified as stage 0, stage 1a (mild perisinusoidal), stage 1b (moderate perisinusoidal), stage 1c (portal/periportal fibrosis only), stage 2 (zone 3 and periportal), stage 3 (bridging fibrosis), and stage 4 (cirrhosis). For analysis, participants with fibrosis stages 1a, 1b, and 1c were grouped into a single category, and fibrosis stage ≥3 was considered advanced fibrosis. The diagnosis of MASH was divided into four categories: no MASH, borderline zone 1, borderline zone 3, and definite MASH72. The borderline and definite MASH categories were also combined into a single group for analysis based on studies in adults indicating a similar risk of fibrosis progression in borderline MASH patients compared to those with definite MASH73.

Unsupervised clustering analysis using clinical data
The clinical variables used for unsupervised clustering were selected based on their association with MASLD, ease of availability in clinical practice, and prior use in adult MASLD phenotyping studies. These ten variables included age (years), ALT (U/L), AST (U/L), WC (cm), VLDL (mg/dL), LDL (mg/dL), TG (mg/dL), HOMA2-IR, uric acid (mg/dL), and SBP (mm Hg), which informed the cluster composition and helped characterize the underlying phenotypes. All variables were scaled to a mean of 0 and a standard deviation of 1 before performing k-means clustering to ensure comparability across different scales. The optimal number of clusters (k) was determined using several approaches, including the silhouette width, and the NbClust package in RStudio version 4.2.3 (R Foundation for Statistical Computing, Vienna, Austria)74, which provides a comprehensive evaluation of 26 different clustering indices to recommend the optimal number of clusters75. Based on these methods, the majority indicated three clusters as optimal (Fig. S15 and S16). K-means clustering was then performed with a k value of three, using the k-means function from the stats package76 in RStudio (maximum iterations: 1000, nstart: 25), and results were visualized with the factoextra R package77. The k-means algorithm partitions the data by minimizing within-cluster variance using a centroid-based distance metric for efficient unsupervised clustering78. Permutation testing was also conducted to assess cluster stability by generating 100 randomized datasets, each created by independently permuting the values within each variable (Fig. S1). K-means clustering was applied to each permuted dataset, and the average silhouette width was calculated for cluster numbers ranging from k=2 to 10. Since sex is a categorical variable and could not be included in k-means clustering, principal component analysis (PCA) was used to reduce the data to the first two principal components, enabling cluster visualizations to highlight differences by sex (Fig. S1).

Cluster profiling
All variables were evaluated for normality, and descriptive statistics were calculated to summarize the characteristics of the three clusters. Categorical variables were reported as counts and frequencies, while continuous variables were presented as medians with interquartile ranges (IQR). Pairwise comparisons between clusters were performed using the Mann-Whitney U test for continuous variables and either the chi-squared test or Fisher’s exact test, as appropriate, for categorical variables. Key results were illustrated using boxplots and barcharts. The clusters were labeled based on the clinical and histological variables that showed statistically significant differences and best characterized their underlying phenotypes. Statistical analyses were conducted in RStudio version 4.2.3, with statistical significance set at p<0.05. OpenAI’s ChatGPT version 4 was used to refine R code and assist with minor language improvements, including grammar and readability, during manuscript preparation. The final manuscript, including all scientific content and conclusions, was thoroughly reviewed and validated by the authors.

High-resolution metabolomics (HRM)
Participant samples from each study were safely transferred from the NIDDK Central Repository to Emory University, where they were stored at -80 °C. While serum samples from the DB2 cohort were processed in 2020 and those from the TONIC and DB1 cohorts in 2022, all analyses were performed at the Emory Clinical Biomarkers Laboratory using established liquid chromatography-mass spectrometry (LC-MS) methods, including identical instrumentation, protocols, and internal quality control procedures to minimize potential batch-related variability32,79. Additionally, serum samples were randomized prior to extraction and LC-MS acquisition to minimize potential bias from instrument drift and batch position effects. All samples were extracted with acetonitrile (ACN) plus an internal standard and analyzed in batches of 40 using a dual chromatography platform. The platform consisted of hydrophilic interaction liquid chromatography (HILIC) with positive electrospray ionization (ESI) and C18 reversed-phase chromatography with negative ESI. Serum samples from the TONIC and DB1 studies were analyzed using an Orbitrap Fusion mass spectrometer (Thermo Fisher Scientific, San Jose, CA), while samples from the DB2 study were analyzed on an HFQE mass spectrometer (Thermo Fisher Scientific, San Jose, CA), both operated at 120,000 resolution in full-scan mode. Each batch of samples included pooled human plasma (QStd-3) at the beginning, middle, and end, along with NIST 1950-calibrated reference pooled human plasma preceding and following each batch. This dual chromatography platform provides broad coverage of the metabolome, including ~500 endogenous, dietary, and environmental chemicals80, with identities confirmed by co-elution and MS/MS fragmentation spectra matching authentic standards32. Raw data files were extracted to generate feature tables using apLCMS81 and xMSanalyzer82, which facilitate adaptive alignment of untargeted metabolomics data and merging of features across multiple extraction parameters. Spectral feature intensities were corrected for batch effects using ComBat83 and filtered for the coefficient of variations less than 75%. Triplicate analyses were tested for repeatability using Pearson correlation (r>0.7), and triplicates were median summarized prior to further analysis. The resulting data file from the TONIC and DB1 studies consisted of 18,826 HILIC/+ESI and 18,024 C18/-ESI chemical features defined by an accurate mass to charge (m/z), retention time (RT), and ion abundance, while the resulting HILIC/+ESI and C18/-ESI data files from the DB2 study consisted of 12,614 and 8468 m/z features, respectively. The two separate HILIC+ and C18- datasets were matched and merged based on an m/z cutoff of 5 ppm and an RT cutoff of 30 seconds using RStudio software version 4.2.3. The resulting 3758 m/z-matched features from HILIC+ and 3520 from C18- were standardized across the studies, quantile normalized, and log2-transformed using the online software MetaboAnalyst version 6.0 prior to downstream analysis. A summary of our analytical workflow is illustrated in Fig. 6.

Metabolomics data analysis and feature annotation
Principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were used to examine overall patterns and group separation among the three metabotypes. PCA was initially applied to identify potential outliers based on the Hotelling’s T² statistic with a 99% confidence interval across the first five principal components84. Outlier removal, performed in RStudio version 4.2.3, led to the exclusion of 24 samples (4.7% of the dataset), possibly due to pre-analytical variations among samples, such as participant location or sample freeze times, for which specific data were not available. One-way ANOVA was performed in MetaboAnalyst version 6.0, both with and without the identified outliers, to confirm that their removal did not alter the overall results. Significant features identified by ANOVA, following FDR correction, were further analyzed using Tukey’s Honest Significant Difference (HSD) post hoc test for pairwise comparisons. The top 25 differentiating HILIC+ and C18− features were first annotated using an internal reference library and confirmed by matching ion dissociation patterns and retention times with authentic standards, which were reported as Schymanski Level 185. The remaining features were computationally annotated using xMSannotator, which performs accurate mass matching to common positive and negative mode adducts in the Human Metabolome Database (HMDB) with an m/z tolerance of ±5 parts per million (ppm) and a retention time tolerance of 10 seconds86. Annotation confidence was reported using the Schymanski Level Confidence (SLC) framework, which ranges from Level 1 (confirmed structure with reference standard) to Level 5 (exact mass of interest only)32.

Pathway enrichment was conducted with Mummichog version 2.087 in MetaboAnalyst 6.088. Significant features identified through one-way ANOVA with FDR-adjusted p values <0.05 were used as input. Enrichment was evaluated via permutation testing using all quality-filtered m/z features as the background. Pathways were considered enriched if the permutation-derived p value was <0.05 and at least three significant features overlapped with known metabolites in the pathway. Pathway membership was defined according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) library. Analyses were performed with the following parameters: mass tolerance 5 ppm; RT in seconds; primary ions enforced; and inclusion of pathways/metabolite sets containing ≥3 entries. Significant features identified through pathway analysis were also annotated as described above. The pathway results include ratios that indicate the number of significant features mapped to each pathway (numerator) out of the total number of pathway-relevant features (denominator). Group-wise distributions of the significant features were examined to interpret pathway differences between groups and visualized using the R packages pheatmap89 for heatmaps and ggplot290 for violin plots. Pathway enrichment was performed as an exploratory approach to investigate potential biological differences between clinically defined metabotypes, recognizing that further targeted studies will be necessary to confirm these findings. The datasets used as input for pathway analysis via MetaboAnalyst are provided in Supplementary Data 9 and 10.

Integrative network analysis and pathway enrichment analysis
The 284 significant metabolic features from HILIC+ and 178 from C18-, identified through one-way ANOVA, were integrated with the clinical dataset using the data integration software xMWAS (v.0.552)91. To reduce skewness, the following clinical variables were log-transformed prior to data integration: ALT, AST, TG, VLDL, LDL, and HOMA2-IR. The fibrosis stage was also included in the clinical dataset due to its significant variation across the three clusters, and was square root transformed to address zero values. In xMWAS, partial least squares (PLS) canonical mode was used to conduct pairwise association analysis between the metabolomics and clinical datasets for each metabotype. This approach identifies bidirectional associations between clinical variables and metabolic features92. The clinical dataset was entered as the reference, and associations were selected based on association scores (|r|>0.3) and significance (p<0.05 by Student’s t-test). Network visualization and community detection were performed for each metabotype based on the resulting association matrix using a multilevel algorithm designed to identify hierarchical community structures consisting of tightly connected nodes representing metabolic features and clinical variables93. The eigenvector centrality method was utilized to identify nodes exhibiting changes in their topological characteristics between networks. Nodes with higher centrality scores have more connections and are considered more influential within the network94,95. An absolute difference in centrality score (Δ centrality) of 0.025 was selected to identify metabolic features with meaningful changes in centrality across networks. This threshold was selected to equilibrate loading of selected metabolites for pathway analysis across all comparisons while limiting inclusion of minimal centrality shifts. The features with a Δ centrality of ≥0.025 between network pairs were entered into pathway enrichment analysis using Mummichog version 2.087 via MetaboAnalyst 6.088. Pathways were considered significantly enriched based on p<0.05 in permutation-based testing and an overlap size ≥3 m/z features.

Correlation analysis of metabolic features with fibrosis stage
To identify metabolic features significantly correlated with the fibrosis stage, we used the rcorr function from the Hmisc package96 in RStudio to calculate Spearman correlation coefficients and p values for metabolites in both HILIC+ and C18- modes. The fibrosis stage was included as a continuous variable and square root transformed to address zero values and reduce skewness. The top 20 significant feature correlations (p<0.05) were annotated using an internal reference library or xMSannotator, as described above, and visualized in a Manhattan plot, with p values transformed to negative log10 to highlight significance levels.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The study investigates the Metabolic Dysfunction-Associated Steatotic Liver Disease (MALDS) in 3 pediatric cohorts assembled into one dataset containing demographic data, clinical targeted measurements and untargeted metabolomics features. The authors present the analytical chemistry method for metabolic feature discovery and the data analysis pipeline. In the first step of the data analysis, clinically distinct metabotypes were formed with unsupervised clustering based on the 10 manually selected features. Following that, the comparative analysis of metabolomics data was performed to identify the most distinct metabolic features of each cluster. Some features were annotated and used in pathway enrichment analysis. The authors present important work on pediatric cohorts in a clear and comprehensive way. The reviewer is not a wet lab specialist and therefore cannot comment on the specifics of the LC-MS method. Thus, the comments are focused on the data analysis methods.\n1. The central claim of the paper is the existence of the three distinct metabotypes, which directs the follow up statistical and pathway enrichment analysis. This analysis is performed based on the 10 preselected variables. While the unsupervised methods used for clustering and arriving at the conclusion about the existence of exactly 3 metabotypes are solid, the reviewer has concerns about the possible bias introduced by the choice of variables. In the Discussion section, the authors make connections between the findings of the presented study and the existing literature. However, the 10 variables are also preselected on the basis of “prior use in adult MASLD phenotyping studies”, which makes them more likely to return similar results to the existing studies. The reviewer suspects that the results of the unsupervised clustering might significantly change when all the clinical and demographic variables are included. The results of the clustering that are based on all the variables must be included and compared to the presented metabotypes. Otherwise, the more detailed justification for not including the variables must be described.\n2. Do the same metabotypes also occur when the unsupervised clustering of untargeted metabolomics features is performed? Do some other clusters form if only the metabolomics is measured and clinical variables are not included? How do they relate to the metabotypes formed by the clinical variables and the previous studies?\n3. No participants from DB1 ended up in the inflammatory-fibrotic metabotype. Do the authors suspect that this due to chance or due to the selection bias of the participants in DB1?\n4. To make an analysis on the Figure 4 even more explainable, the reviewer suggests including some potential cofactors (age, clinical variables) into a linear regression alongside with the most correlated metabolic features to predict the levels of fibrosis. This would allow to list out some possible cofactors and establish associations in the datasets.\n5. The reviewer suggests that the Supplementary Tables are provided in the Excel/CSV format for easy access.", "reviewer_2": "In this study the authors used a very well characterized cohort of 514 children with biopsy-proven MASLD to identify clinically distinct metabolic phenotypes of paediatric MASLD using unsupervised clustering and high-resolution metabolomics. This work with the aim to uncover distinct metabolic phenotypes in children with MASLD is very timely and important and a significant step in paediatric liver disease. In the first step, they used clinical and histological parameters and identified 3 unique clusters. This part is done very well. Although validation using external cohort is required . The second part related to metabolomics is more concerning. 3 cohorts were done in two groups possibly in two different labs and 2 years apart, hence the initial number of features are quite different and then no randomization in analysis which is quite important. The methodology is also less suitable for analysing lipids which are important in MASH and important potential biomarkers were reported previously in adult MASH patients. In addition, the results are not well-described. They reported 284 significant features but it is not clear if these are the metabolites that differentiate the 3 identified clusters or did they performed unsupervised analysis independently? Are these fully annotated using standards (tier 1) or these are only features or mixed? Then without discussing specificity of these metabolites, they moved on with performing pathway analysis using Mummichog or MetaboAnalyst which are ok for quick visualization but not reliable for such an important differentiation. I really don’t think this type of analysis is relevant for this work and it offers any meaningful outcome. I suggest to remove it! This study has a lot of potential with a very important cohort. It wold be great to use the same thorough evaluation of the metabolites as it was done for the clinical and histological parameters.\nMinor comment: Please use the recently approved term MASH/MASLD across the manuscript. Currently both MASH and NASH/NAFLD is used across the text.", "reviewer_3": "Thank you for asking me to review this manuscript. This was a very interesting to read, and relevant to the field, piece of work.\nNoteworthy results:\nThe work presented looks to identify distinct phenotypes/metabotypes of paediatric MASLD by analysing clinical and metabolomic data from children with biopsy proven MASLD who had previously participated in 3 NASH CRN studies. (The inclusion of biopsy proven MASLD is a strength of this work).\nWith the use of unsupervised clustering 3 distict metabotypes were identified: the early-mild, the cardiometabolic and the inflammatory-fibrotic, with differences in clinical, laboratory and metabolic parameters between the groups. In addition six metabolic pathways were shown to be significantly altered between the three metabotypes.\nSignificance of work for the field:\nThis work would be of significance in the field, as it is very important to understand the different 'types' of MASLD in children from a disease pathogenesis point of view, from a prognostic point of view and of course from a management/ treatment point of view. Even more so as the prevalence of paediatric mASLD is rapidly increasing world wide. The work compares well to the relevant literature in adults with MASLD and the authors analyse this adequately in the discussion.\nThe work does support the conclusions and claims. The authors do highlight the limitations of the study, mainly the lack of a control group and the predominance of Hispanic males in the participants. This does limit generalisability of the conclusions, but the work still provides valubale insight that would be very helpful in the design of future studies including longitudinal research.\nIn terms of data analysis, interpretation and conclusions, there are a few points that may strengthen the manuscript should the authors wish to address them:\n1. For parameters like waist circumference and systolic blood pressure, I believe it would be more appropriate to use z- scores for age and sex, particularly as there were significant differences in terms of age between clusters and there is an over-representation of male participants.\n2. As mentioned on line 143/144 and seen in Table 1, the early-mild metabotype has a significantly higher percentage of patients with advanced fibrosis when compared to the cardiometabolic metabotype and about 50% of the percentage of advance fibrosis seen in the inflammatory-fibrotic metabotype. Knowing that fibrosis is the driver of outcomes from a liver point of view and knowing that this group was actually younger in age than the other 2 groups, makes it more worrying that this group potentially has more aggressive disease, rather than mild disease. How would the authors comment?\n3. As mentioned on line 193/194 and Figure S5, tryptophan metabolism is enriched mostly in the inflammatory-fibrotic group, then in the cardiometabolic group and then in the early-mild metabotype (albeit smaller differences between the CM and the EM group in comparison to the IF group). Knowing the association of tryptophan metabolism and fibrosis stage how would the authors reconcile that the fibrosis patterns don't quite follow this pattern amongst the groups (i.e. fibrosis was more in the IF>EM>CM)?\n4. How do the authors believe that the exclusion of diabetic children from the study impact on the results?\n5. Would it be possible, to have some suggested way from the authors of cluster prediction based on some of the clinical parameters? In other words, when presented with a patient how can we know which metabotype they potentially fit in?"}, "author_response": "Dear Nature Communications Editorial Team and Reviewers,\nThank you for the careful evaluation of our manuscript and the constructive feedback. We have\nthoroughly addressed all comments and revised the manuscript accordingly.\nAll changes in the revised manuscript are highlighted and indicated in bold. We have also\ncompleted the required checklists and other materials in accordance with Nature\nCommunications policies.\nWe look forward to your further consideration of our revised submission.\nSincerely,\nHelaina Huneault, PhD, RD (on behalf of the research team)\nReviewer #1 (Remarks to the Author):\nThe study investigates the Metabolic Dysfunction-Associated Steatotic Liver Disease (MASLD)\nin 3 pediatric cohorts assembled into one dataset containing demographic data, clinical targeted\nmeasurements, and untargeted metabolomics features. The authors present the analytical\nchemistry method for metabolic feature discovery and the data analysis pipeline. In the first step\nof the data analysis, clinically distinct metabotypes were formed with unsupervised clustering\nbased on the 10 manually selected features. Following that, the comparative analysis of\nmetabolomics data was performed to identify the most distinct metabolic features of each cluster.\nSome features were annotated and used in pathway enrichment analysis.\nThe authors present important work on pediatric cohorts in a clear and comprehensive way. The\nreviewer is not a wet lab specialist and therefore cannot comment on the specifics of the LC-MS\nmethod. Thus, the comments are focused on the data analysis methods.\n1. The central claim of the paper is the existence of the three distinct metabotypes,\nwhich directs the follow-up statistical and pathway enrichment analysis. This analysis is\nperformed based on the 10 preselected variables. While the unsupervised methods used for\nclustering and arriving at the conclusion about the existence of exactly 3 metabotypes are\nsolid, the reviewer has concerns about the possible bias introduced by the choice of\nvariables. In the Discussion section, the authors make connections between the findings of\nthe presented study and the existing literature. However, the 10 variables are also\npreselected on the basis of “prior use in adult MASLD phenotyping studies”, which makes\nthem more likely to return similar results to the existing studies. The reviewer suspects that\nthe results of the unsupervised clustering might significantly change when all the clinical\nand demographic variables are included. The results of the clustering that are based on all\nthe variables must be included and compared to the presented metabotypes. Otherwise, the\nmore detailed justification for not including the variables must be described.\nWe thank the reviewer for these thoughtful comments. We agree that variable selection can\ninfluence clustering results and appreciate the opportunity to clarify our approach. As described\nin our Methods section, we intentionally selected 10 clinically accessible variables based on their\nprior use in adult MASLD phenotyping studies, relevance to MASLD pathophysiology, and\navailability in routine clinical care. (Methods section, Page 16, Lines 433-437) This approach\nensured that the identified metabotypes would be clinically interpretable, reproducible, and\ncomparable to prior MASLD phenotyping studies.\nThe importance of consistent variable selection is well demonstrated in diabetes sub-phenotyping\nresearch. Since the initial study by Ahlqvist et al.1, researchers have consistently used the same\nset of variables to perform clustering across different populations. This approach has led to novel\ninsights; for example, in European populations, clusters tend to be driven by insulin resistance,\nwhereas in South Asian populations, defects in insulin secretion are the predominant factors.\nSuch comparisons reveal distinct underlying mechanisms in South Asians, which would have\nbeen difficult to identify/compare if different variable sets were used. Conversely, including a\nlarge number of additional variables can obscure biologically meaningful signals, introduce\nnoise, and reduce the clarity of clustering outcomes. This concern is especially relevant for k-\nmeans clustering, where the Euclidean distance metric is sensitive to the number of variables; as\nthis number increases, the meaning of distance changes, often leading to fewer and less distinct\nclusters2.\nAs suggested by the reviewer, we conducted an exploratory clustering analysis using all\navailable continuous clinical and demographic variables (n=21), including markers such as GGT,\nHbA1c, creatinine, bilirubin, and additional lipid measures. As shown below, the NbClust R\npackage identified k=3 (a consensus of 26 metrics), and the silhouette method suggested k = 2\nclusters. Despite this variability, the 3-cluster solution aligned most closely with the clinical\nphenotypes observed in our primary 10-variable model. Inclusion of additional variables led to\nminor changes in cluster membership but did not improve interpretability or alter the overall\nclustering structure.\nNbClust (consensus of 26 clustering metrics):\nSilhouette method:\nBased on these findings, we retained the 10-variable clustering model to preserve clinical\nrelevance, improve interpretability, and support comparability with prior MASLD phenotyping\nstudies.\n2. Do the same metabotypes also occur when the unsupervised clustering of untargeted\nmetabolomics features is performed? Do some other clusters form if only the metabolomics is\nmeasured and clinical variables are not included? How do they relate to the metabotypes\nformed by the clinical variables and the previous studies?\nWe thank the reviewer for these important questions. In this study, clustering was intentionally\nperformed using clinical variables to identify clinically meaningful and translatable pediatric\nMASLD phenotypes. Untargeted metabolomics was then used to explore metabolic differences\nand compare these clinically defined metabotypes. This two-step framework is now clarified in\nour updated Figure 6 (shown below). Additionally, we updated our methods section to highlight\nthis: Methods section, Pages 13, Lines 365-369: “This retrospective cross-sectional study\nemployed a two-step framework (Figure 6). First, a data-driven clustering analysis was\nperformed using anthropometric and clinical data from youth enrolled in NASH CRN\nstudies to identify metabotypes of pediatric MASLD. Second, an exploratory high-\nresolution metabolomics analysis of participant serum samples was conducted to examine\ndifferences in the metabolome among the identified metabotypes.” We did not perform\nunsupervised clustering based solely on untargeted metabolomics features, as this was beyond\nthe scope of our analysis.\nFigure 6. Two-step framework for identifying metabotypes and exploring metabolic\ndifferences in pediatric MASLD. (1) Unsupervised clustering of clinical variables from 514\nchildren enrolled in NASH CRN studies with biopsy-confirmed MASLD identified distinct\nmetabotypes using k-means clustering. (2) High-resolution metabolomics was performed on\nstored serum samples, including feature extraction (apLCMS, xMSanalyzer), normalization,\nand\nlog2 transformation. Features that differed significantly across metabotypes (p < 0.05, FDR-\nadjusted) were identified using one-way ANOVA, entered into untargeted pathway analysis\n(Mummichog v2.0), and annotated using an internal reference library (tier 1) and xMSannotator.\nMetabolomics features were integrated with clinical biomarkers and fibrosis stage in xMWAS for\nnetwork and pathway analysis. Features differing in centrality scores by ≥0.025 between\nmetabotype networks were further explored in pathway enrichment analysis.\nWhile we did not cluster on metabolomics data alone, our study did examine how untargeted\nmetabolomics features differed across the three clinically defined metabotypes. This integrative\napproach allowed us to explore underlying biological differences, such as alterations in\ntryptophan metabolism, branched-chain amino acid degradation, and purine metabolism, which\nprovided mechanistic insights into the distinct clinical profiles of the metabotypes (Discussion\nsection, Pages 9-11, Lines 272-324). We agree that future studies could explore clustering\ndirectly on metabolomics data; however, such analyses would benefit from advanced statistical\nmethods and integration with clinical variables to ensure biological and clinical relevance.\nFurthermore, we have discussed how the metabotypes formed in our analysis relate to previously\nreported subtypes in our Discussion section, Pages 8-9, Lines 243-271)\n3. No participants from DB1 ended up in the inflammatory-fibrotic metabotype. Do the\nauthors suspect that this is due to chance or due to the selection bias of the participants in\nDB1?\nWe appreciate this thoughtful observation. We suspect that the absence of DB1 participants in the\ninflammatory-fibrotic metabotype is primarily due to sampling variability, given the small\nnumber of DB1 participants included in the present analysis (n=23 with complete clinical and\nmetabolomics data). With such a limited sample size, it is likely that the distribution across the\nthree identified metabotypes did not fully capture the heterogeneity present in the broader DB1\ncohort. While we cannot entirely rule out the possibility of selection bias, we believe the lack of\nDB1 representation in the inflammatory-fibrotic group is more likely due to sample size\nconstraints than to a systematic bias in participant selection.\n4. To make an analysis on the Figure 4 even more explainable, the reviewer suggests\nincluding some potential cofactors (age, clinical variables) into a linear regression alongside\nwith the most correlated metabolic features to predict the levels of fibrosis. This would\nallow to list out some possible cofactors and establish associations in the datasets.\nWe thank the reviewer for this suggestion. While we agree that incorporating clinical covariates\ninto a regression model alongside selected metabolite features to predict fibrosis severity could\nprovide additional insights, the primary focus of this study was to define clinically relevant\nmetabotypes of pediatric MASLD and characterize their distinct metabolomic profiles. Our goal\nwas not to identify individual metabolites associated with fibrosis stage or to build a predictive\nmodel for fibrosis, as this represents a separate research question beyond the scope of the current\nanalysis. We acknowledge that modeling fibrosis by integrating clinical and metabolomics data is\nan important future direction, and we have added this to our Discussion section as an avenue for\nfuture research in pediatric MASLD. Discussion section, Page 13, Lines 354-356: “Moreover,\nfurther research could explore the relationship between specific metabolite features,\nclinical variables, and fibrosis severity in pediatric MASLD to identify potential\nbiomarkers of fibrosis risk and refine our understanding of disease progression.”\n5. The reviewer suggests that the Supplementary Tables are provided in the Excel/CSV\nformat for easy access.\nWe appreciate the reviewer’s suggestion. We have now provided all Supplementary Tables in\nExcel format to facilitate ease of access.\nReviewer #1 (Remarks on code availability):\nThe code is provided alongside with installation instructions.\nThe reviewer tested one file and it was executing as expected.\nReviewer #2 (Remarks to the Author):\n1. In this study the authors used a very well characterized cohort of 514 children with biopsy-\nproven MASLD to identify clinically distinct metabolic phenotypes of paediatric MASLD using\nunsupervised clustering and high-resolution metabolomics. This work with the aim to uncover\ndistinct metabolic phenotypes in children with MASLD is very timely and important and a\nsignificant step in paediatric liver disease. In the first step, they used clinical and histological\nparameters and identified 3 unique clusters. This part is done very well. Although validation\nusing external cohort is required . The second part related to metabolomics is more concerning. 3\ncohorts were done in two groups possibly in two different labs and 2 years apart, hence the initial\nnumber of features are quite different and then no randomization in analysis which is quite\nimportant. The methodology is also less suitable for analysing lipids which are important in\nMASH and important potential biomarkers were reported previously in adult MASH patients. In\naddition, the results are not well-described. They reported 284 significant features but it is not\nclear if these are the metabolites that differentiate the 3 identified clusters or did they performed\nunsupervised analysis independently? Are these fully annotated using standards (tier 1) or these\nare only features or mixed? Then without discussing specificity of these metabolites, they moved\non with performing pathway analysis using Mummichog or MetaboAnalyst which are ok for\nquick visualization but not reliable for such an important differentiation. I really don’t think this\ntype of analysis is relevant for this work and it offers any meaningful outcome. It would be great to use the same thorough evaluation of the metabolites as it was done for the clinical and histological parameters.\nWe thank the reviewer for recognizing the strengths of our clinical and histological analyses and\nfor their detailed critique of the metabolomics component. We appreciate the opportunity to\nclarify and strengthen this aspect of the study and have made several key revisions and additions\nas detailed below.\nWe agree that validation in an external cohort is an important next step, and we have noted this\nas a limitation in the Discussion section, Page 11-12, Lines 326- 329.\nWe clarify that while samples from the DB2 cohort were processed in 2020 and those from the\nTONIC and DB1 cohorts in 2022, all analyses were conducted at the Emory Clinical Biomarkers\nLaboratory using consistent instrumentation, protocols, and internal quality control procedures.\nThis minimized potential batch-related variability. To clarify this, we have added the following\nsentence to our Methods section (Page 17, Lines 468-474): “While serum samples from the\nDB2 cohort were processed in 2020 and those from the TONIC and DB1 cohorts in 2022, all\nanalyses were performed at the Emory Clinical Biomarkers Laboratory using established\nliquid chromatography-mass spectrometry (LC-MS) methods including identical\ninstrumentation, protocols, and internal quality control procedures to minimize potential\nbatch-related variability.”\nAdditionally, as indicated in the Methods section (Page 18, Lines 493-497: “The two separate\nHILIC+ and C18- datasets were matched and merged based on a m/z cutoff of 5 ppm and\nan RT cutoff of 30 seconds using R Studio software version 4.2.3. The resulting 3758 m/z-\nmatched features from HILIC+ and 3520 from C18- were standardized across the studies,\nquantile normalized, and log2-transformed using the online software MetaboAnalyst\nversion 6.0 prior to downstream analysis.” Furthermore, we clarify that all serum samples\nwere randomized within each batch before extraction and LC-MS acquisition to reduce potential\nbias from instrument drift and batch effects (Methods, Page 17, Lines 473–474). “Additionally,\nserum samples were randomized prior to extraction and LC-MS acquisition to minimize\npotential bias from instrument drift and batch position effects.”\nWe have also updated Figure 6 to outline the two-step analysis framework used in this study (see\nresponse to reviewer #1, comment #2 above). Unsupervised clustering was first performed solely\non clinical variables, followed by metabolomics analysis, which was conducted as a separate,\nexploratory analysis to characterize metabolic differences between these predefined\nmetabotypes.\nMoreover, we have improved transparency regarding metabolite annotation. Among the\nsignificant features identified through one-way ANOVA (284 features from HILIC+ and 178\nfrom C18-), the top 25 have now been annotated as shown in Supplemental Tables S2 and S3.\nThese features were either confirmed using an internal reference library and authentic standards\n(MSI Level 1) or computationally annotated (MSI Level 2 or lower) following the Metabolomics\nStandards Initiative guidelines as described in the updated table descriptions. The significant\nfeatures identified through pathway analysis were also annotated using the same process as\nindicated in our updated Methods section. Methods section (Page 19, Lines 509-520): “The top\n25 differentiating HILIC+ and C18− features were annotated using an internal reference\nlibrary, confirmed by matching ion dissociation patterns and elution times with authentic\nstandards3. These matches were considered level 1 ‘confirmed’ per the Metabolomics\nStandards Initiative (MSI)4. The remaining m/z features were computationally annotated\nusing xMSannotator, which performs accurate mass matching to common positive and\nnegative mode adducts in HMDB with an m/z tolerance of ± 5 parts per million (ppm) and\na retention time tolerance of 10 seconds5. Pathway enrichment was conducted with\nMummichog version 2.06 in MetaboAnalyst 6.07, with features identified through one-way\nANOVA with FDR-adjusted p-values below 0.05 to assess pathway-level differences.\nGroup-wise distributions of significantly enriched features were examined to interpret\npathway differences between groups and visualized using the R packages pheatmap8 for\nheatmaps and ggplot29 for violin plots. Significant features identified through pathway\nanalysis (p<0.05) were also annotated as described above.”\nWe acknowledge that the pathway enrichment analysis conducted using Mummichog and\nMetaboAnalyst is exploratory and intended to provide initial insights into potential biological\ndifferences across metabotypes, rather than definitive evidence of specific pathways or\nbiomarkers. We have added a clarification to the Methods section (Page 19, Lines 521-523) to\nreflect this, stating: “Pathway enrichment analysis was performed as an exploratory\napproach to investigate potential biological differences between the clinically defined\nmetabotypes, recognizing that further targeted studies are necessary to confirm these\nfindings.” In addition, the following statement has been added to the Discussion/Limitations\nsection (Page 13, Lines 357-358): “Finally, our pathway enrichment analysis was\nexploratory and not intended to establish definitive pathways or biomarkers; targeted\nstudies are required to validate these findings.”\nWe also acknowledge the limitations of our untargeted metabolomics platform for analyzing\nlipids. Several lipid-related features, including phosphatidylcholine and sphingolipid species,\nwere differentially abundant across metabotypes. While our dual LC-MS platform (HILIC+ and\nC18−) approach offers broad coverage of the metabolome, it is not optimized for targeted\nlipidomics. To address this, we have added the following sentence to the Discussion (Page 12-\n13, Lines 351-354): “While our untargeted metabolomics approach offers broad coverage\nof the metabolome, it is not optimized for targeted lipidomics. Future studies employing\ndedicated lipidomics platforms are needed to better characterize lipid-driven mechanisms\nunderlying pediatric MASLD heterogeneity.”\nWe believe these substantial revisions, including the updated Figure 6 workflow, the clarified\nanalysis strategy, the annotation of the top 25 features, and the expanded transparency around\nmethods and limitations, address the reviewer’s concerns and significantly strengthen the\nmetabolomics component of our study. We thank the reviewer for their constructive feedback,\nwhich has improved the quality and clarity of our manuscript.\n2. Minor comment: Please use the recently approved term MASH/MASLD across the\nmanuscript. Currently both MASH and NASH/NAFLD is used across the text.\nThank you for this helpful suggestion. We have carefully revised the manuscript to use the\nupdated terminology MASH (metabolic dysfunction-associated steatohepatitis) and MASLD\n(metabolic dysfunction-associated steatotic liver disease) throughout, in accordance with the\nrecently approved nomenclature. References to NASH/NAFLD have been retained only when\nreferring to historical study names (e.g., NASH CRN) or previously published data that used the\nformer terminology. These instances are now clarified in the text to avoid confusion.\nReviewer #3 (Remarks to the Author):\nThank you for asking me to review this manuscript. This was a very interesting to read, and\nrelevant to the field, piece of work.\nNoteworthy results:\nThe work presented looks to identify distinct phenotypes/metabotypes of paediatric MASLD by\nanalysing clinical and metabolomic data from children with biopsy proven MASLD who had\npreviously participated in 3 NASH CRN studies. (The inclusion of biopsy proven MASLD is a\nstrength of this work).\nWith the use of unsupervised clustering, 3 distinct metabotypes were identified: the early-mild,\nthe cardiometabolic, and the inflammatory-fibrotic, with differences in clinical, laboratory, and\nmetabolic parameters between the groups. In addition six metabolic pathways were shown to be\nsignificantly altered between the three metabotypes.\nSignificance of work for the field:\nThis work would be of significance in the field, as it is very important to understand the different\n'types' of MASLD in children from a disease pathogenesis point of view, from a prognostic point\nof view and of course from a management/ treatment point of view. Even more so as the\nprevalence of paediatric MASLD is rapidly increasing world wide. The work compares well to\nthe relevant literature in adults with MASLD and the authors analyse this adequately in the\ndiscussion.\nThe work does support the conclusions and claims. The authors do highlight the limitations of\nthe study, mainly the lack of a control group and the predominance of Hispanic males in the\nparticipants. This does limit the generalizability of the conclusions, but the work still provides\nvaluable insight that would be very helpful in the design of future studies, including longitudinal\nresearch.\nIn terms of data analysis, interpretation and conclusions, there are a few points that may\nstrengthen the manuscript should the authors wish to address them:\n1. For parameters like waist circumference and systolic blood pressure, I believe it would\nbe more appropriate to use z-scores for age and sex, particularly as there were significant\ndifferences in terms of age between clusters and there is an over-representation of male\nparticipants.\nWe thank the reviewer for this helpful suggestion. While age- and sex-adjusted percentiles or z-\nscores are commonly used in pediatric studies, their application in our MASLD cohort is limited.\nThe 2017 American Academy of Pediatrics guidelines provide SBP percentiles based on normal-\nweight children, which may not generalize to youth with overweight or obesity (Flynn et al.,\n2017)10. Similarly, the z-score equations from Wühl et al. (2002)11 were derived from normal-\nweight European children and may not apply to diverse or overweight populations. For WC, no\nwidely accepted z-score references exist across the full pediatric age and BMI spectrum\n(Yamanaka et al., 2021; Fernández et al., 2004)12,13.\nIn our study, age was included as a clustering variable, and all features, including SBP and WC,\nwere standardized prior to clustering to ensure equal weighting, as described in the Methods\n(Page 16, Lines 437–439): “All variables were scaled to a mean of 0 and a standard\ndeviation of 1 before performing k-means clustering to ensure comparability across\ndifferent scales.” This approach maintains clinical interpretability and aligns with prior MASLD\nphenotyping studies. While male participants were overrepresented, sex was not a primary driver\nof cluster separation based on principal component analysis (Figure S2).\n2. As mentioned on line 143/144 and seen in Table 1, the early-mild metabotype has a\nsignificantly higher percentage of patients with advanced fibrosis when compared to the\ncardiometabolic metabotype and about 50% of the percentage of advanced fibrosis seen in\nthe inflammatory-fibrotic metabotype. Knowing that fibrosis is the driver of outcomes from a\nliver point of view and knowing that this group was actually younger in age than the other 2\ngroups, makes it more worrying that this group potentially has more aggressive disease, rather\nthan mild disease. How would the authors comment?\nWe thank the reviewer for raising this important observation. We hypothesize that the higher\nproportion of advanced fibrosis in the early-mild metabotype, despite its younger age profile,\nmay reflect a mixed, earlier-stage phenotype that shares features of both the cardiometabolic and\ninflammatory-fibrotic groups. This group may include individuals who could progress toward\neither of these more defined phenotypes over time (Figure 5).\nAdditionally, while the early-mild group shows fewer cardiometabolic risk factors overall, it is\nimportant to note that fibrosis can develop early in life, even in the absence of overt metabolic\ndysfunction. The presence of advanced fibrosis in this younger group may therefore indicate an\nunderlying vulnerability or disease trajectory that warrants close monitoring.\nWe have added this point to the Discussion to highlight the need for future longitudinal studies to\nbetter understand progression pathways and the clinical significance of the early-mild group.\nDiscussion section, Page 8, Lines 228-230: “As illustrated in Figure 5, we hypothesize that\nthe early-mild metabotype represents a younger, less severe phenotype with mixed features\nof the other two metabotypes, which may progress into one of them over time.” Also, Page\n8, Lines 239-242: “The higher proportion of advanced fibrosis observed in the early-mild\ngroup underscores the potential for fibrosis to develop early in life, even in the absence of\novert metabolic dysfunction, and highlights the need for longitudinal studies to better\nunderstand progression pathways in pediatric MASLD14.”\n3. As mentioned on line 193/194 and Figure S5, tryptophan metabolism is enriched mostly\nin the inflammatory-fibrotic group, then in the cardiometabolic group and then in the\nearly-mild metabotype (albeit smaller differences between the CM and the EM group in\ncomparison to the IF group). Knowing the association of tryptophan metabolism and\nfibrosis stage how would the authors reconcile that the fibrosis patterns don't quite follow\nthis pattern amongst the groups (i.e. fibrosis was more in the IF>EM>CM)?\nWe thank the reviewer for this thoughtful observation. While tryptophan metabolism was most\nenriched in the inflammatory-fibrotic group, followed by the cardiometabolic and early-mild\ngroups, the fibrosis distribution followed a different pattern (IF > EM > CM). We hypothesize\nthat this reflects the complexity of tryptophan metabolism, which is involved in multiple\nbiological processes beyond fibrosis, including obesity-associated low-grade chronic\ninflammation, metabolic dysfunction, and immune regulation. Therefore, while tryptophan\npathway enrichment is highest in the inflammatory-fibrotic group, its presence in the\ncardiometabolic group may reflect metabolic inflammation rather than fibrosis per se.\nWe have added this interpretation to the Discussion to emphasize the need for future targeted\nstudies to clarify the specific roles of tryptophan metabolism in pediatric MASLD progression.\nDiscussion section, Page 10, Lines 284-290: “Tryptophan catabolism is also upregulated in\nobesity and is associated with systemic low-grade chronic inflammation15,16. While\nenrichment of tryptophan metabolism in the inflammatory-fibrotic group likely reflects\nfibrosis-related processes, its elevated presence in the cardiometabolic group may reflect\nunderlying metabolic inflammation. Fibrosis progression in pediatric MASLD likely results\nfrom the combined effects of altered tryptophan metabolism, lipid metabolism, and\ninflammation. Future studies are needed to clarify the specific contributions of tryptophan\nmetabolism to fibrosis development and disease progression in this population.”\n4. How do the authors believe that the exclusion of diabetic children from the study impact\non the results?\nWe appreciate the reviewer’s thoughtful question. As shown in the participant selection\nworkflow (Supplementary Figure S14), participants with diabetes (HbA1c ≥6.5%) were not\nexplicitly excluded during initial data selection. However, because k-means clustering is\nsensitive to extreme values, we applied an objective outlier removal process, excluding\nparticipants with any clinical variable greater than five standard deviations from the mean. This\napproach resulted in the exclusion of all but one participant with an HbA1c ≥6.5%. To maintain\nconsistency and avoid undue influence from a single extreme value, we also excluded this\nremaining participant prior to analysis. This exclusion slightly reduced the sample size but\nimproved the robustness of the clustering by minimizing the influence of extreme outliers.\nWe have clarified this step in the Methods section (Page 14, Lines 390-397): \"The initial\neligible sample included 556 children and adolescents with available clinical and\nmetabolomics data. A small percentage of missing data (0.34%) was imputed using the\nmedian for clinical variables. Outliers, defined as participants with any clinical variable\ngreater than five standard deviations from the mean, were excluded, resulting in 515\nparticipants. During this process, most participants with an HbA1c ≥6.5% were excluded;\nhowever, one participant remained in the dataset. To maintain consistency and avoid bias\nfrom a single extreme value, this participant was also excluded, yielding a final sample of\n514 children and adolescents for analysis (Figure S14).”\n5. Would it be possible, to have some suggested way from the authors of cluster prediction\nbased on some of the clinical parameters? In other words, when presented with a patient\nhow can we know which metabotype they potentially fit in?\nThank you for this insightful suggestion. Developing a clinical tool for predicting metabotype\nmembership is an important focus of our ongoing work. Similar to the approach by Raverdy et\nal. (Nature Medicine, 2024)17, we plan to create a calculator that will allow clinicians to input a\npatient's values for the 10 clinical variables used in clustering and compute the Euclidean\ndistance to each cluster centroid. The patient would then be assigned to the metabotype with the\nclosest centroid. However, we recognize that this approach requires careful validation, and our\nfirst priority is to validate the metabotypes identified in this study in an independent cohort.\nOnce validated, we will refine and test the calculator, which we plan to include in a future publication."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The paper investigates the Metabolic Dysfunction-Associated Steatotic Liver Disease (MASLD)\nin 3 pediatric cohorts. In the review round, the authors addressed the questions about the data analysis, mostly in the domain of potential bias introduced by variable pre-selection. The authors conducted new experiments that include more variables and a different method benchmark, updated relevant text and figures. Overall, the comments are addressed and clarity improved.", "reviewer_2": "Thank you for addressing the previous comments and revising the manuscript. The updated version is significantly improved. However, I remain concerned about the pathway analysis, which plays a central role in the biological interpretation of the metabolomics data in this manuscript.\nWhile pathway analysis can offer useful insights into metabolic alterations, it must be applied rigorously. In this manuscript, the approach appears to rely on publicly available tools that, while convenient, often lack transparency and reproducibility, and may yield misleading results without careful validation.\nFor example, in Figure 2, \"Tryptophan metabolism\" is listed as a top pathway (6/41), and \"Purine metabolism\" as another (4/70). However, what do these ratios signify? Does \"6/41\" mean that only 6 of the 41 known metabolites in the tryptophan pathway were statistically significant? If so, that would be a relatively weak signal and may not justify emphasizing the pathway as enriched. The criteria used for pathway selection — including enrichment metrics, statistical thresholds, and background sets — are not clearly described and need to be explicitly stated.\nSimilarly, in Figure 4A, the metabolite PA (20:4) is highlighted. What is this molecule? Is it phospholipids???! It does not appear among the directly annotated metabolites and seems to be inferred only through the pathway tool. This raises questions about annotation confidence and interpretability.\nThe authors have really nice data. Using appropriate statistical analysis with correctly annotated metabolites, they can assess the biological importance without introducing any bias.", "reviewer_3": "Thank you for the revision submitted.\nI am happy with the response to my comments."}, "author_response": "Manuscript Title: Clinically distinct metabotypes of pediatric MASLD identified through\nunsupervised clustering of NASH CRN data\nManuscript ID: NCOMMS-25-11773A\nCorresponding Author: Helaina Huneault, PhD, RD\nDate: August 28, 2025\nDear Nature Communications Editorial Team and Reviewers,\nThank you for your thoughtful feedback and the opportunity to revise our manuscript. We have\ncarefully addressed all remaining reviewer comments, with particular attention to clarifying our\npathway analysis and annotation methods. Below, we provide a detailed point-by-point response.\nAll changes in the revised manuscript are highlighted and indicated in bold.\nWe appreciate your continued consideration.\nSincerely,\nHelaina Huneault, PhD, RD (on behalf of the research team)\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe paper investigates the Metabolic Dysfunction-Associated Steatotic Liver Disease (MASLD)\nin 3 pediatric cohorts. In the review round, the authors addressed the questions about the data\nanalysis, mostly in the domain of potential bias introduced by variable pre-selection. The authors\nconducted new experiments that include more variables and a different method benchmark,\nupdated relevant text and figures. Overall, the comments are addressed and clarity improved.\nReviewer #1 (Remarks on code availability):\nThe code is provided alongside with installation instructions.\nThe reviewer tested one file and it was executing as expected\nThank you for reviewing our manuscript and code.\nReviewer #2 (Remarks to the Author):\nThank you for addressing the previous comments and revising the manuscript. The updated\nversion is significantly improved. However, I remain concerned about the pathway analysis,\nwhich plays a central role in the biological interpretation of the metabolomics data in this\nmanuscript.\nWhile pathway analysis can offer useful insights into metabolic alterations, it must be applied\nrigorously. In this manuscript, the approach appears to rely on publicly available tools that, while\nconvenient, often lack transparency and reproducibility, and may yield misleading results\nwithout careful validation.\nFor example, in Figure 2, \"Tryptophan metabolism\" is listed as a top pathway (6/41), and \"Purine\nmetabolism\" as another (4/70). However, what do these ratios signify? Does \"6/41\" mean that\nonly 6 of the 41 known metabolites in the tryptophan pathway were statistically significant? If\nso, that would be a relatively weak signal and may not justify emphasizing the pathway as\nenriched. The criteria used for pathway selection — including enrichment metrics, statistical\nthresholds, and background sets — are not clearly described and need to be explicitly stated.\nSimilarly, in Figure 4A, the metabolite PA (20:4) is highlighted. What is this molecule? Is it\nphospholipids???! It does not appear among the directly annotated metabolites and seems to be\ninferred only through the pathway tool. This raises questions about annotation confidence and\ninterpretability.\nThe authors have really nice data. Using appropriate statistical analysis with correctly annotated\nmetabolites, they can assess the biological importance without introducing any bias.\nWe thank the reviewer for their careful assessment and constructive feedback regarding pathway\nenrichment analysis and metabolite annotation. We have taken several steps to improve the\ntransparency of our methods and clarify the interpretation of our results.\nWe agree that the ratios in Figure 2 were not clearly explained in the initial submission. While\nthis information was included in the figure legend, it was missing from the Methods. In the\nrevised manuscript, we now explicitly state that the numbers in parentheses indicate the number\nof statistically significant metabolic features (numerator) out of the number of detected metabolic\nfeatures within each pathway (denominator), as defined by the KEGG pathway library. In line\nwith default parameters in Mummichog and MetaboAnalyst, we required at least three\nsignificant features per pathway to reduce the likelihood of spurious enrichment results.\nWe have also revised the Methods section to clearly describe the enrichment criteria, statistical\nthresholds, and background sets used for pathway analysis. Methods section, Pages 19-20,\nLines 527–542: “Pathway enrichment was conducted with Mummichog version 2.01 in\nMetaboAnalyst 6.02. Significant features identified through one-way ANOVA with FDR-\nadjusted p-values <0.05 were used as input. Enrichment was evaluated via permutation\ntesting using all quality-filtered m/z features as the background. Pathways were considered\nenriched if the permutation-derived p-value was <0.05 and at least three significant\nfeatures overlapped with known metabolites in the pathway. Pathway membership was\ndefined according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) library.\nAnalyses were performed with the following parameters: mass tolerance 5 ppm; RT in\nseconds; primary ions enforced; and inclusion of pathways/metabolite sets containing ≥3\nentries. Significant features identified through pathway analysis were also annotated as\ndescribed above. The pathway results include ratios that indicate the number of significant\n\nfeatures mapped to each pathway (numerator) out of the total number of pathway-relevant\nfeatures (denominator). Group-wise distributions of the significant features were examined\nto interpret pathway differences between groups and visualized using the R packages\npheatmap3 for heatmaps and ggplot24 for violin plots. Pathway enrichment was performed\nas an exploratory approach to investigate potential biological differences between clinically\ndefined metabotypes, recognizing that further targeted studies will be necessary to confirm\nthese findings.”\nAdditionally, we updated Supplementary Table S4 to include the Schymanski Level Confidence\n(SLC) for each annotation (see Liu et al., 2020, Analytical Chemistry). Several of the\nsignificantly enriched pathways contained at least one Level 1 annotation from our internal\nreference library, providing strong confidence in these results. Three pathways, propanoate\nmetabolism, BCAA degradation, and pantothenate and CoA biosynthesis, included features such\nas 3-methyl-2-oxobutanoic acid and 4-methyl-2-oxopentanoate that could not be corroborated by\neither xMSannotator or our internal reference library. These features were therefore\nconservatively reported as SLC Level 5. We have noted these as exploratory, hypothesis-\ngenerating findings that warrant confirmation in future targeted studies. The datasets used as\ninput for pathway analysis are now provided in Supplementary Tables S9 and S10, and our\nMethods and Results sections have been updated accordingly.\nMethods section, Page 19, Lines 542-543: “The datasets used as input for pathway analysis\nvia MetaboAnalyst are provided in Supplementary Tables S9 and S10.”\nResults section, Pages 5-6, Lines 174-180: “Among the significantly enriched pathways,\nseveral contained at least one Schymanski Level Confidence (SLC) Level 1 identification\nfrom our internal reference library, supporting the confidence of these results5. Three\npathways, propanoate metabolism, BCAA degradation, and pantothenate and CoA\nbiosynthesis, contained features that could not be independently confirmed by secondary\nannotation methods and were therefore reported as SLC Level 5. These results should be\ninterpreted cautiously but may represent true biological signals and provide important\nhypotheses for future validation.”\nWe have also updated the explanation of our metabolite annotation procedures for clarity.\nMethods section, Page 19, Lines 518–526: “…features were first annotated using an\ninternal reference library and confirmed by matching ion dissociation patterns and\nretention times with authentic standards, which were reported as Schymanski Level 16. The\nremaining features were computationally annotated using xMSannotator, which performs\naccurate mass matching to common positive and negative mode adducts in the Human\nMetabolome Database (HMDB) with an m/z tolerance of ± 5 parts per million (ppm) and\na retention time tolerance of 10 seconds7. Annotation confidence was reported using the\nSchymanski Level Confidence (SLC) framework, which ranges from Level 1 (confirmed\nstructure with reference standard) to Level 5 (exact mass of interest only)5.\nRegarding Figure 4A, the feature initially described as PA (20:4) has been corrected to PA(O-\n22:4), consistent with the HMDB11156 candidate returned by xMSannotator. PA refers to\nphosphatidic acid, with O- denoting an ether-linked species. Because this annotation is based on\naccurate mass without MS/MS confirmation and without determination of an unequivocal\nmolecular formula, the feature is conservatively reported as SLC Level 5. We clarify that all\nmetabolite features in Figure 4 were annotated using either our internal reference library or\nxMSannotator, as described above. This is now stated in our updated Methods.\nMethods section, Page 21, Lines 571–574: “The top 20 significant feature correlations (p <\n0.05) were annotated using an internal reference library or xMSannotator, as described\nabove, and visualized in a Manhattan plot, with p-values transformed to negative log10 to\nhighlight significance levels.”\nWe have also updated Supplementary Figure S8 to include the compound name, adduct, and\nSLC for each of the top 20 annotated features from HILIC+ and C18– shown in Figure 4, with\nthe corrected annotation listed as PA(O-22:4), SLC Level 5. To improve interpretability, we\nadded an abbreviation key to the legend of Figure 4: “Figure 4. Top 20 HILIC+ (A) and C18–\n(B) metabolite features associated with fibrosis among NASH CRN study participants.\nCorrelations are presented on a color scale from −1 (blue) to 1 (red). The dashed line\nindicates a p-value cutoff of <0.05. Metabolite features were annotated using the internal\nreference library and xMSannotator. Please see Table S8 for correlations, p-values,\nadducts, and annotation confidence. Abbreviations: PA, phosphatidic acid; PC,\nphosphatidylcholine; PE, phosphatidylethanolamine; PI, phosphatidylinositol. O- denotes\nan ether-linked species.”\nLastly, we have updated our Discussion section to emphasize the exploratory nature of the\nmetabolomics analysis and findings. Discussion section, Page 13, Lines 368-372: “Our\nfindings underscore the potential for precision healthcare approaches to tailor\ninterventions to specific metabolic dysfunctions in MASLD subtypes. Given the exploratory\nnature of the metabolomics analysis, further longitudinal and targeted studies are needed\nto validate these findings and explore strategies that may improve outcomes for children\nwith MASLD.”\nWe sincerely appreciate the reviewer’s comments, which helped us significantly strengthen the\nclarity and interpretability of our pathway analysis and feature annotation methods.\nReviewer #3 (Remarks to the Author):\nThank you for the revision submitted.\nI am happy with the response to my comments.\nThank you for reviewing our manuscript."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors acknowledged the exploratory nature of pathways analysis and included the description of the methodology used. I have no further comments."}, "author_response": "Responses to the Reviewers\nClinically distinct metabotypes of pediatric MASLD identified through unsupervised\nclustering of NASH CRN data\nNature Communications manuscript NCOMMS-25-11773B\nDecember 3, 2025\nReviewer #2 (Remarks to the Author):\n1) The authors acknowledged the exploratory nature of pathways analysis and included the\ndescription of the methodology used. I have no further comments.\nWe thank the reviewer for their feedback and are pleased that the clarification of the exploratory\npathway analysis and methodological description addressed their concerns. No further changes\nwere required in response to this comment."}]

Decision letter