This retrospective study utilized plasma samples collected before 10am from both female and male individuals provided by the UK ME/CFS Biobank (UKMEB) [17, 18]. The female cohort comprised 56 individuals, categorized into HC (n = 12; 21%), ME/CFSmm (n = 20; 36%), and ME/CFSsa patients (n = 24; 43%). The male cohort included 31 individuals, distributed into HC (n = 17; 55%), ME/CFSmm (n = 8; 26%), and ME/CFSsa patients (n = 6; 19%) (Table 1).

Individuals diagnosed with ME/CFS were evaluated by a clinician based on the Canadian Consensus [19] and/or CDC-1994 (‘Fukuda’) [20] criteria. The diagnosis was validated through responses on the Symptoms Assessment form, ensuring compliance with the case definition and study eligibility. Participants in the study completed a set of questionnaires aimed at evaluating disability levels, including the Fatigue Severity Scale, measuring the severity of fatigue symptoms to provide insights into their impact. Additionally, the Pain and Fatigue Analog Scale assessed the subjective experience of pain and fatigue on a visual analog scale, allowing participants to express the intensity of their symptoms. The Medical Outcomes Survey Short Form (SF-36v2) was employed as a comprehensive multidimensional instrument to evaluate various aspects of participants' physical and mental well-being, offering a holistic view of health-related quality of life. These instruments collectively contributed to a thorough assessment of participants' health status and functional components [21, 22].

Exclusion criteria for participants included individuals who, within the preceding three months, (1) had used drugs known to modify immune function (e.g., azathioprine, cyclosporine, methotrexate, steroids) or had taken antiviral medications; (2) had received any vaccinations; (3) had a history of acute or chronic infectious diseases such as hepatitis B and C, tuberculosis, or HIV (excluding infections by herpes virus or other retroviruses); (4) had another severe medical condition such as cancer, coronary heart disease, or uncontrolled diabetes; (5) had a severe mood disorder; (6) had been pregnant or breastfeeding in the prior 12 months; or (7) presented with morbid obesity (BMI ≥ 40). Home visits were conducted to recruit patients with mobility restrictions (severely affected), while healthy subjects and mild/moderate patients were invited to a recruiting center for clinical assessment and blood sampling [21, 22].

Certified reference standards of corticosterone solution (C-117-1ML; CAS: 50-22-6), cortisol solution (C-106-1ML; CAS: 50-23-7), cortisone solution (C-130-1ML; CAS: 53-06-5), and 11-deoxycortisol solution (D-061-1ML; CAS: 152-58-9) were purchased from Merck (KgaA, Darmstadt, Germany). Stable isotope-labeled standards, namely aldosterone-D7 (CAS: 1261254-31-2), androstene-3,17-dione-13C3 (CAS: 327048-86-2), corticosterone-D4 (CAS: 2243253-91-8), cortisol-D4 (CAS: 73565-87-4), cortisone-13C3 (CAS: 2350278-95-2), 11-deoxycortisol-D5 (CAS: 1258063-56-7), 17α-hydroxyprogesterone-D8 (CAS: 850023-80-2), progesterone-D9 (CAS: 15775-74-3), and testosterone-13C3 (CAS: 327048-83-9) were also obtained from Merck (KgaA, Darmstadt, Germany). Purified water was prepared in-house using a Milli-Q water system from Millipore (Bedford, MA, USA). HPLC-grade methanol (MeOH; 34885-1L-M), ethyl acetate (EtOAc; 1.00868), methyl-tert-butyl ether (MTBE; 34875-1L-M), LC–MS-grade acetonitrile (ACN; 1.00029. 1000), 2-propanol (IPA; 1.02781.1000), formic acid (FA; 1002531000), and ammonium fluoride (338869-25G) were purchased from Merck (KGaA, Darmstadt, Germany).

To prepare the stock solution of the internal standard (IS), the stable isotope-labeled compounds were combined in ACN to achieve final concentrations of 2000 ng/mL for aldosterone-D7, cortisone-13C3 and cortisol-D4, 400 ng/mL for corticosterone-D4 and 11-deoxycortisol-D5, 200 ng/mL androstene-3,17-dione-13C3, testosterone-13C3, 17α -hydroxyprogesterone-D8, and 100 ng/mL for progesterone-D9. The working solution (WSL) was obtained by diluting the stock 1:10 in MeOH. The standard stock solution mixture in MeOH with a concentration of 1000 ng/mL was prepared using the nine individual steroid hormones. A 500 μL aliquot of the stock solution mixture was transferred to a 5 mL volumetric flask and brought up to the highest calibration point with a 50:50 (v:v) mixture of MeOH and H2O. Nineteen additional dilutions were made using a 50:50 (v:v) MeOH:H2O mixture, and the calibration curve samples were prepared by adding 10 μL of the IS WSL solution to 100 μL of the mixed standard solution to cover a calibration range from 0.00019 ng/mL to 100 ng/mL. Stock solutions were stored at –20°C and allowed to reach room temperature before use.

Supported liquid extraction (SLE) was employed to mitigate the impact of the sample matrix. The SLEs were purchased from Agilent Technologies (5610-2005). Briefly, 100 µL of plasma were combined with 100 µL of H2O and 10 µL of IS WSL. The resultant mixture was transferred to the SLE tube, placed onto the sorbent bed with gentle pressure (2–3 psi), and allowed to equilibrate for 5 min. Subsequently, 400 µL of a 1:1 mixture of methyl-tert-butyl ether (MTBE) and ethyl acetate (EtOAc) were introduced into each tube, followed by elution at a rate of 1 drop per second (2 psi). This elution process was iterated three times, and a final application of 6 psi was employed to desiccate the sorbent. The entire eluent was dried with nitrogen (N2) flow at 40°C in a TurboVap water bath (Biotage, Sweden) and then reconstituted with 100 µL of MeOH. Samples were stored at 4 °C overnight, subjected to centrifugation at 2500×g for 5 min, and subsequently loaded into the autosampler.

All separations were performed on a 1290 Infinity UHPLC system (Agilent Technologies) equipped with an Agilent ZORBAX RRHD Eclipse Plus C18 column (2.1 × 100 mm; 1.8 µm; 821725-902) and a ZORBAX RRHD C18 guard column (2.1 × 5 mm; 1.8 µm; 821725-901). The mobile phases comprised 0.2 mM ammonium fluoride in H2O as mobile phase A and 0.2 mM ammonium fluoride in MeOH as mobile phase B. The gradient conditions were as follows: 0–3.0 min; 50–60% B; 3.0–7.0 min, 60–86% B; 7.0–7.1 min; 86-100% B, followed by a return to the initial conditions. The total chromatographic run time was 8.5 min. The flow rate was set to 0.4 mL/min, and the column temperature was maintained at 40 °C. The injection volume was 3 µL, and a needle wash with 1:1:1:1 ACN/MeOH/IPA/H2O with 0.2% FA was utilized. Mass detection was carried out in dynamic multiple reaction monitoring (dMRM) mode on an Agilent 6460 triple quadrupole system using positive electrospray ionization (ESI) mode. Specific settings can be found in Tables S1 and S2, while Figure S1 provides a representative UHPLC-dMRM chromatogram.

For comparing the means of clinical variables among the three groups (HC, ME/CFSmm, and ME/CFSsa), one-way ANOVA was employed when the respective data followed a normal distribution. In cases where normal distribution was not met, the non-parametric Kruskal–Wallis test was used. The significance level for these tests was set at 5%. Central tendency and variability in the dataset were estimated using the median and interquartile range (IQR).

To assess differences in the circulating levels of steroid hormones in female and male ME/CFS patients compared to HC, ANOVA, and Kruskal–Wallis tests were utilized to determine statistical significance. Unadjusted p values were computed for each test. To address multiple testing, we controlled the false discovery rate (FDR) at a 5% level using the Benjamini–Hochberg (BH) and Benjamini–Yekutieli (BY) procedures.

We also conducted a predictive analysis to assess the disease status (HC versus ME/CFS) of each participant using the following classifiers and the nine steroid hormones as respective predictors: (i) linear discriminant analysis (LDA), (ii) random forest (RF), and (iii) partial least square discriminant analysis (PLS-DA). This analysis was carried out separately for male and female datasets. In PLS-DA, the analysis was performed with one and two latent components. Using a higher number of components in PLS-DA, while potentially enhancing predictive performance, was prone to overfitting the data. For LDA and PLS-DA, the probability of an individual being an ME/CFS patient was estimated based on a leave-one-out procedure. For RF, the same probability was estimated using Bootstrap with 10,000 simulated decision trees.

After obtaining the classification probability for each participant, we constructed a Receiver Operating Characteristic (ROC) curve, where 1-specificity (x-axis) was plotted against sensitivity (y-axis). The corresponding area under the curve (AUC) and the point on the curve closest to the pair (0,1) (representing perfect classification) were determined. This point served as the optimal cutoff to predict the health status of each individual. Finally, we utilized the predicted health status to estimate accuracy (proportion of individuals with correctly predicted health status), sensitivity (Se, proportion of cases with correctly predicted health status), and specificity (Sp, proportion of controls with correctly predicted health status) associated with each classifier. This analysis was carried out in the R software using the following packages: MASS (for LDA) [23], randomForest (for RF) [24], caret (for PLS-DA) [25], pROC (for ROC curve and AUC calculation) [26], OptimalCutpoints (for optimal cutoff estimation, accuracy, sensitivity, and specificity) [27].

Analyses utilizing Spearman's correlation coefficient (Rsp) were conducted to examine the statistical relationship between data from two steroid hormones. Each steroid hormone was ranked from lowest to highest, and the correlation coefficient along with corresponding p values were computed using GraphPad Prism. Next, to compare Spearman's correlation matrices between HC and individuals with ME/CFS (or its subgroups), with or without the inclusion of data on aldosterone (which contains missing data), a permutation variant of Jennrich's test was implemented. This test was originally designed for comparing two Pearson's correlation matrices [28]. The computational procedure involved the following steps: (i) fixing the values of the group variable; (ii) permutating the original data set; (iii) computing the Jennrich's test statistic in the permutated data set; (iv) repeating steps (ii) and (iii) until obtaining 1000 values of the test statistic; (v) estimating the p value by determining the proportion of times the observed test statistic in the original dataset was higher than the values of the test statistic based on the permutated datasets.

Principal component analysis (PCA) of the steroid hormones data was employed to evaluate the similarity of study participants concerning sex and disease severity. In this analysis, we utilized the first and second principal components, as these components explained over 90% of the data variability. Heatmaps were generated for the entire dataset and separately for the female and male cohorts using the MetaboAnalyst R package [29] to conduct a thorough analysis of the data. Clustering of the data was performed using Euclidean distance, aiming to unveil underlying patterns and relationships. To identify distinct subgroups, the resulting clusters were further segmented into multiple groups. Data visualization was executed using the pheatmap R package [30].