Between May 2022 and January 2023, 43 participants (30 paediatric patients and 13 healthy volunteers) were recruited from the Paediatric University Hospital La Paz (Madrid, Spain). A discovery cohort of 29 participants included eight healthy paediatric volunteers (HVs) and 21 children with acute hepatitis. Of these, 17 patients met the ECDC/WHO criteria for PAHUA, which considers hepatitis virus serology, transaminase levels, and age ( Figure 1 ). Four children were excluded due to confirmed alternative diagnoses of autoimmune hepatitis (AIH). The AIH group served as an internal control for conditions with similar clinical presentations. A validation cohort was established with 7 PAHUA patients, 2 AIH patients, and 5 HVs to confirm the predictive model. All blood samples were collected from the participants prior to treatment.

Fresh blood from all participants was collected in EDTA-anticoagulated tubes (Fisher Scientific, Massachusetts, United States, #Cat: 10081844). Whole white blood cells were enriched following red blood cell lysis. Briefly, 2 mL of EDTA-anticoagulated venous blood were treated with 20 mL of 1X Pharm Lyse Buffer (BD Biosciences, New Jersey, United States, #Cat: 555899) for 15 minutes at room temperature (RT) on a rocker. After two PBS (Phosphate Buffered Saline, Sigma, Massachusetts, United States, #Cat: P4417-100TAB) washes, cells were stained with LIVE/DEAD Fixable Blue Dead Cell Stain Kit (Invitrogen, Massachusetts, United States, #Cat: L23105) to exclude dead cells. Prior to labelling, True-Stain Monocyte Blocker (BioLegend, California, United States, #Cat:426103) was added to prevent nonspecific binding of some fluorochromes to monocytes. Cells were then labelled with fluorochrome-conjugated monoclonal antibodies for 25 minutes at RT in the dark. Two separate cytometry panels ( Supplementary Tables S1 , S2 ) were used to analyse immune cell populations and immune checkpoint expression, respectively. Labelled cells were acquired using a Cytek Aurora Spectral Cytometer (Cytek Biosciences, California, United States). Two additional cytometry panels ( Supplementary Tables S3 , S4 ) were used to assess CD4/CD8 proliferation and memory profiles and monocyte phenotyping. Data were analysed using FlowJo software (v10.6.2, TreeStar).

Peripheral blood mononuclear cells (PBMCs) were extracted from the same EDTA-anticoagulated venous blood samples used for immunophenotyping in both paediatric patients and healthy paediatric volunteers. The isolation was performed using Ficoll-Plus (Cytiva, Massachusetts, United States, #Cat: 17144003) in accordance with the manufacturer’s instructions. The PBMCs were washed twice with phosphate-buffered saline (PBS), stained with Trypan blue, and counted. Freshly isolated PBMCs were cultured in RPMI 1640 medium (ThermoFisher Scientific, Massachusetts, United States, #Cat: 72400112) supplemented with 10% foetal bovine serum (FBS, ThermoFisher Scientific, Massachusetts, United States, #Cat: A5669701), 25 mM HEPES, 2 mM L-glutamine, and 1% Penicillin-Streptomycin (ThermoFisher Scientific, Massachusetts, United States, #Cat: 15140122) prior to stimulation to induce cellular activation or proliferation. The cultures were maintained at 37°C in a humidified incubator with 5% CO₂.

Plasma samples were collected from EDTA-anticoagulated venous blood using a Ficoll-Plus (Cytiva, Massachusetts, United States, #Cat: 17144003) solution, following standard density gradient centrifugation. The plasma samples were aliquoted and stored at -80°C until further use.

Plasma samples were thawed, and centrifuged at 1000 RCF for 30 minutes to remove particulates before use. The titres of SARS-CoV-2-specific IgG antibodies against the Spike (S1 and receptor binding domain [RBD]) and Nucleocapsid (N) proteins were measured using the LEGENDplex SARS-CoV-2 Serological IgG Panel (Biolegend, California, United States, #Cat: 741132), following the manufacturer’s instructions. The samples were analysed using a FACSCalibur flow cytometer (BD Biosciences, New Jersey, United States), and data were processed using LEGENDplex software (v.8, Biolegend).

PBMCs were labelled with Carboxyfluorescein succinimidyl ester (CFSE, ThermoFisher Scientific, Massachusetts, United States, #Cat: C34554) according to the manufacturer’s protocol and seeded into 96-well round-bottom plates (2x10⁵ cells per well) in RPMI 1640 medium (ThermoFisher Scientific, Massachusetts, United States, #Cat: 72400112) supplemented with 10% foetal bovine serum (FBS, ThermoFisher Scientific, Massachusetts, United States, #Cat: A5669701), 25 mM HEPES, 2 mM L-glutamine, and 1% Penicillin-Streptomycin (ThermoFisher Scientific, Massachusetts, United States, #Cat: 15140122). CFSE-labelled PBMCs were stimulated with either pokeweed mitogen (PWD, 2.5 μg/Ml, Sigma, Massachusetts, United States, #Cat: L8777-5MG) or Dynabeads Human T-activator CD3/CD28 (0.5 μL/well, ThermoFisher Scientific, Massachusetts, United States, #Cat: 11131D) for three days. Supernatants were collected, aliquoted, and stored at -80°C until use.

PBMCs were seeded into 96-well flat-bottom plates (based on 2x10⁴ monocytes per well) in RPMI 1640 medium supplemented (ThermoFisher Scientific, Massachusetts, United States, #Cat: 72400112) supplemented with 10% foetal bovine serum (FBS, ThermoFisher Scientific, Massachusetts, United States, #Cat: A5669701), 25 mM HEPES, 2 mM L-glutamine, and 1% Penicillin-Streptomycin (ThermoFisher Scientific, Massachusetts, United States, #Cat: 15140122). Cells were stimulated with lipopolysaccharide (LPS, Sigma, Massachusetts, United States, #Cat: L2630-10MG) to simulate Gram-negative bacterial infection, R848 (Resiquimod, InvivoGen, California, United States, #Cat: tlrl-r848-1) as a TLR7/8 agonist to simulate viral infection, or beta-glucan (β-glucan, InvivoGen, California, United States, #Cat: tlrl-wgp) to simulate fungal infection for 24 hours (early response) and seven days (late response). After stimulation, supernatants were collected, aliquoted, and stored at -80°C.

Soluble immune checkpoints and cytokine concentrations in plasma samples and supernatants were measured using the LEGENDplex Human Immune Checkpoint Panel 1 (12-plex: sCD25, 4-1BB, sCD27, sCD86, sCTLA-4, sPD-L1, sPD-L2, sPD-1, sTim-3, sLAG-3, Gal-9 and TGF-β; Biolegend, California, United States, #Cat: 740867) and LEGENDplex Human Essential Immune Response Panel (13-plex: IL-1β, IL-2, IL-4, IFNγ, TNF-α, MCP-1, CXCL10, IL-6, IL-8, IL-10, IL-12p70, IL-17A, and Free Active TGF-β1; Biolegend, California, United States, #Cat: 740930), according to the manufacturer’s instructions. Samples were acquired on a FACSCalibur flow cytometer (BD Biosciences, New Jersey, United States), and data were processed using LEGENDplex software (v.8, Biolegend).

Raw data from immune factors and COVID-19 history were normalised using the Z-score method. Hierarchical clustering was performed using R (v.4.0.2) with the heatmap, geom_tile, and ggplot2 packages (v.1.16.0). The overall accuracy was calculated by dividing the number of correctly classified samples in the two-group analysis by the total number of samples. Principal Component Analysis (PCA) was conducted using the prcomp, princomp and PCA packages with k-nearest neighbours (kNN) and Euclidean clustering.

Group differences were assessed using Student’s t-test and Levene’s test, with statistical significance set at p<0.05. Receiver Operating Characteristic (ROC) curve analysis identified potential predictors, with optimal cut-offs determined by the Youden Index.

Due to the nature of the SPSS software, the binary logistic regression model was performed on a subset of cases that presented all the available variables, comprising 10 PAHUA patients and 11 non-PAHUA cases. The binary logistic regression model predicted PAHUA using nine variables with AUC/ROC>0.9 and statistically significant in univariate regression applying Wald forward stepwise regression (13, 14). The final model combined Gal-9 and sTim-3, expressed as: Score = 0.00178 x Gal-9 levels + 0.03321 x sTim3 levels.

To further evaluate the predictive validity and generalisability of the binary logistic regression model, internal validation was conducted using R (version 4.3.3). A 10-fold cross-validation was carried out with the caret package, assessing predictive performance via the area under the curve (AUC), which yielded a mean AUC of 1. Sensitivity and specificity were both estimated at 0.9.

Prior to validation, variable selection was performed using the MASS package, applying a forward stepwise logistic regression approach based on the Wald statistic. In addition, a non-parametric bootstrap procedure comprising 10,000 iterations was implemented with the boot package to estimate the stability and 95% confidence intervals of the AUC, resulting in a mean AUC of 0.999.

Receiver operating characteristic (ROC) curves and AUC values were calculated using the pROC package. The dplyr and readxl packages were employed for data pre-processing and importation, respectively.

Statistical analyses were performed with SPSS v.23 and GraphPad Prism 8.0 and R software.

Between May 2022 and January 2023, 43 children were recruited from the Paediatric University Hospital La Paz. A discovery cohort of 29 participants included 21 children with acute hepatitis and eight healthy paediatric volunteers (HVs). Of these, 17 met PAHUA criteria per ECDC/WHO guidelines, while 4 were diagnosed with autoimmune hepatitis (AIH), which served as an internal control ( Figure 1 ).

Most PAHUA patients presented with gastrointestinal symptoms such as vomiting, diarrhoea, and abdominal discomfort, occasionally accompanied by hepatomegaly or splenomegaly. These symptoms progressed to jaundice, with fever notably absent ( Supplementary Table S5 ). Hepatitis viruses A–E were excluded, and no shared environmental exposures were identified.

The cohort’s median age was 7.2 years (range: 1–15), with gender distribution balanced (52% male) ( Table 1 ). PAHUA patients exhibited significantly elevated bilirubin and liver enzyme levels (AST, ALT, and GGT) compared to HVs ( Table 1 ). Anaemia, coagulopathy, abnormal INR and prothrombin activity were common among PAHUA cases. Ten PAHUA and one AIH patient developed paediatric acute liver failure (pALF); one required liver transplantation, and two succumbed to the illness. Treatments included plasmapheresis (24%), steroids (12%), intravenous immunoglobulin (IVIG) (6%), and N-acetylcysteine (6%).

Adenovirus was detected in 17.6% of PAHUA patients, with two cases co-infected with influenza A ( Supplementary Table S6 ). Additionally, 41% tested positive for other infections, including norovirus, cytomegalovirus (CMV), and herpesvirus. SARS-CoV-2 IgG antibodies against the spike protein were found in 88.2% of PAHUA patients and 100% of AIH and HVs, indicating prior infection rather than vaccination ( Supplementary Table S6 ). Nucleocapsid protein IgG analysis revealed that 11.8% of PAHUA patients lacked evidence of prior SARS-CoV-2 infection, supporting hypotheses of immune-mediated liver inflammation potentially linked to the spike protein acting as a superantigen (11).

To investigate immune profiles in PAHUA and non-PAHUA (AIH and HV) groups, immunophenotyping of whole blood was performed using spectral flow cytometry, alongside ex vivo immune response assessments, T cell proliferative capacity, and quantification of cytokine and soluble immune checkpoints (ICs) levels ( Supplementary Figure S1 ).

Combining all immunological data, principal component analysis (PCA) revealed distinct clustering patterns among PAHUA (red: 1-17), AIH (green: 18-21), and HV (blue: a-h) groups ( Figure 2A ). AIH cases were most distinct, separating undoubtedly from PAHUA and HV groups. Two-group heatmap analysis using Euclidean clustering and correlation grouping achieved 93.1% overall accuracy, segregating non-PAHUA cases (top cluster) from PAHUA patients (bottom cluster). Two HVs (a and c cases) grouped with PAHUA patients, likely due to prior infection with SARS-CoV-2 (a) and HIV exposure during pregnancy (c) ( Figure 2B ).

Analyses focusing on specific datasets showed varying accuracies. Whole blood immunophenotyping, analysing 529 immune cell populations, achieved 86.2% accuracy in the two-group heatmap analysis, though group separation was less distinct in the PCA ( Supplementary Figure S2 ). In this analysis, the top cluster included most of the non-PAHUA patients, with two PAHUA patients (7 and 13). As observed in the previous heatmap ( Figure 2B ), which incorporated all immunological data, two HVs (a and c) were again grouped with the PAHUA patients in the bottom cluster. In contrast, PCA and heatmap analyses based on T cell proliferation and monocyte phenotyping following various stimulations yielded worst results, achieving an overall accuracy of 65.5% ( Supplementary Figure S3 ). The two-group heatmap revealed that most non-PAHUA cases clustered in the top group, which also included seven PAHUA patients (1–6 and 15). The bottom cluster composed of the remaining PAHUA patients and three HVs (a, b and c).

Given the potential of soluble molecules as biomarkers and therapeutic targets (15–17), we assessed their capacity to classify the three groups. ICs achieved 79.3% accuracy, outperforming cytokines (75.9%), but neither matched the full dataset’s accuracy ( Supplementary Figures S4, S5 ). Nevertheless, the showed promise in distinguishing PAHUA from other groups and warrant further investigation as biomarkers.

A volcano plot analysis identified differentially expressed variables between PAHUA, AIH, and HV groups ( Supplementary Figure S6 ). Using a p-value threshold of 0.01 and a fold-change of 2, we found 15 upregulated and eight downregulated variables in PAHUA vs. HV; 28 upregulated and 23 downregulated in PAHUA vs. AIH; and 17 upregulated and 18 downregulated in AIH vs. HV ( Supplementary Figure S6A ).

Immunophenotyping variables were most frequently highlighted and were effective in clustering the groups ( Supplementary Figure S7 ). Notably, Gal-9 and sTim-3 were commonly upregulated in PAHUA compared to both HV and AIH, marking them as key biomarkers ( Supplementary Figure S6B and Supplementary Table S7 ).

To develop a reliable prognostic tool for PAHUA, we assessed 790 variables, including COVID-19 history (3), immune cell populations (529) and ex vivo data (258), using AUC/ROC and univariate logistic regression. Variables with AUCs>0.90 and statistically significant results by univariate logistic regression (p<0.05) were incorporated into a binary logistic regression model ( Figure 3A ).

Nine variables were identified: %CM (central memory) CXCR3⁺CD95⁺ Treg (T regulatory), %EM (effector memory) CXCR3⁺CD95⁺ Treg, %EM CD27^-CD28^- CD4^-CD8^- T cells, sCD25, Gal-9, sTim-3, %CM CXCR3^-CD95⁺ Treg, %Unswitched on B cells, and %CD108 on CD4⁺ T cells. Each variable demonstrated strong discrete predictive capacity, with AUCs ranging from 0.9136 (%CD108 on CD4⁺ T cells) to 0.9818 (sCD25) ( Figure 3B ). Optimal cut-off values were determined using the Youden Index ( Supplementary Figure S8 ).

The nine variables were incorporated into a forward stepwise regression model (Wald test). After three steps, the final optimised model included the combination of two ICs, Gal-9 and sTim-3, both of which are easily quantifiable in clinical settings. This resulted in a highly effective predictive tool, achieving an AUC of 1.000 (95% CI 1.000–1.000) in ROC analysis with a statistically significant p-value of 0.0001 ( Figure 4A ). The optimal cut-off score, determined using the Youden Index, provided maximum sensitivity (1.000; 95% CI 0.723–1.000) and specificity (1.000; 95% CI 0.741–1.000), effectively distinguishing PAHUA cases from non-PAHUA cases (AIH and HV) ( Figure 4B ). It is important to note that into the non-PAHUA individuals, the AIH patients (green dots) were separated totally from the PAHUA patients.

To validate the model, Gal-9 and sTim-3 levels were quantified in the validation cohort (PAHUA, n=7 and non-PAHUA, n=7). The Chi-squared test confirmed significant differences between high and low scoring groups (χ²=14.00; z=3.742; p=0.0002), with perfect sensitivity and specificity in distinguishing PAHUA cases ( Figure 4C ). This finding underscores the discriminatory power of the predictive model in distinguishing between patient groups.

Furthermore, the predictive score demonstrated perfect specificity (1.000; 95% CI 0.646–1.000) and sensitivity (1.000; 95% CI 0.646–1.000) in differentiating PAHUA from non-PAHUA cases ( Figure 4D ).