Longitudinal soluble marker profiles reveal strong association between cytokine storms resulting from macrophage … – Nature.com
Cohort study
The BEAT-COVID cohort was collected at the LUMC, Leiden, The Netherlands. Individuals diagnosed with PCR-positive COVID-19, who required hospital admission between April 2020 and March 2021 were invited to participate17,18,19,20. The majority of infections was with Wuhan-like viruses; the Netherlands experienced<1% circulating alpha variant until January 2021 nationally (in Dutch: https://www.rivm.nl/corona/actueel/virusvarianten).
Clinical care was provided according to local and national guidelines and not influenced by study participation. Initially (April-Aug 2020), care was supportive only (wave-1, March 2020August 2020), however, from August 2020 onwards patients received dexamethasone (wave-2, August 2020March 2021). Since our study was initiated shortly after the onset of the pandemic in The Netherlands, criteria for hospital admission were also different between the waves. In wave-1 hospital admission was limited to the more severe cases, and generally admission occurred later after onset of symptoms compared to wave-2 as a result of capacity and national guidelines. In addition, SARS-CoV2 PCR testing was only available for hospital admitted patients whereas it was mandatory for all individuals with symptoms in the 2nd wave; hence, people were aware of the infection earlier and may thus have behaved differently.
Patients were included both on the regular ward as well as on the ICU. Inclusion criteria were: admission at LUMC, SARS-CoV-2 PCR positive, and minimum 18 years old. Exclusion criteria were: no informed consent from the patient or a representative. All participants were unvaccinated against SARS-CoV-2. Upon signing informed consent, blood samples were collected 3 times (MondayWednesdayFriday) a week during hospitalization and at an outpatient visit 612 weeks post discharge. Statistical sample size calculation was not performed, the sample size was determined based on availability. Anti-IL6R antibodies were only used in a few patients during our study.
In addition to the patients with COVID-19 disease, a control group of healthy individuals was collected in July 2020 for comparative analyses. None of these healthy volunteers had experienced the infection yet, as was confirmed by serological analyses Semi-quantitative detection of SARS-CoV-2 anti-nucleocapsid (N) protein IgG antibodies was performed on the Abbott Architect platform21,22. In this antibody chemiluminescent microparticle immunoassay (CMIA) test, the SARS-CoV-2 antigen coated paramagnetic microparticles bind to the IgG, respectively, IgM antibodies that attach to the viral nucleocapsid protein in human serum samples. The Sample/Calibrator index values of chemiluminescence in relative light units (RLU) of 1.40 (IgG assay) respectively 1.00 (IgM assay) and above were considered as positive per the manufacturers instructions. Samples were collected from 8 males and 4 females (as the ratio of hospital admissions at that time was 2:1), with a median age of 64 (range 6072). Healthy volunteers were sampled 5 times, also at 2 day intervals and samples were processed identically to those from the individuals included in the BEAT-COVID cohort.
Ethical approval for the study protocol was provided by LUMC Medical Ethics committee (protocol NL73740.058.20). The study was registered at the International Clinical Trials Registry Platform no. NL8589 (https://trialsearch.who.int/Trial2.aspx?TrialID=NL8589). The study complied with the latest version of the Declaration of Helsinki. The principal investigator had access to information to identify individual patients.
The Severity of COronavirus Disease Assessment (SCODA) disease severity-score was developed to track day-to-day disease severity in hospital admitted COVID-19 patients23. The score was based on the 4C mortality score, but constant parameters were substituted with parameters associated with breathing and oxygenation to better reflect the daily condition and its changes23. The SCODA score was developed such that there was continuous scoring during admission on the ward and ICU to permit longitudinal assessment of severity during both types of admission. The SCODA score consisted of the following parameters: respiratory rate, peripheral oxygen saturation on room air (ward only), P/F Ratio (ICU only), oxygen flow (ward only), FiO2 (ICU only), Glasgow Coma Scale (GCS), blood urea level and C-reactive protein (CRP)23. The most severely ill patients admitted to the ICU could reach a maximum severity-score of 17.
We defined the recovery time point, per definition after the highest severity-score for that individual patient, as a SCODA score of a daily severity-score of7, with no subsequent increases. SCODA scores were included in the analysis when matched serum samples and measurements on the soluble analytes were available.
Blood samples were collected by venipuncture including a 8.5ml SST Vacutainer tube (Becton Dickinson, VWR, Amsterdam, The Netherlands). Upon clotting samples were centrifuged (10 min, 3000 rpm) within one hour by the central clinical chemistry laboratory and aliquoted for research purposes. Samples were stored at 80C to guarantee optimal sample quality.
The following cytokine and chemokine reagent kits were selected to analyse the sera of the study cohort; the Bio-Plex Pro human Chemokine panel (40-plex, #171AK99MR2), Bio-Plex Pro human Cytokine Screening panel (48-plex, #12,007,283), Bio-Plex Pro Human Inflammation panel (37-plex, #171AL001M) and a custom made panel of IL-17F, IL-21, IL-23, IL-25, IL-31, IL-33 of the Bio-Plex Pro human Th17 cytokine panel (all Bio-Rad, Veenendaal, the Netherlands). Assays were performed according to manufacturers instructions with manufacturers specific standards and QC control samples. All samples were thawed and diluted 1:4 in sample Diluent HB and run as single measurement with the streptavidin PE (1:200, Becton Dickinson, Erembodegem, Belgium) detection label. Samples were acquired on a Bio-Plex 200 system and analysed with Bio-Plex manager software v6.2. In total 131 analytes were measured, in case measurements were out of range or zero for more than 75% of the samples, analytes were removed from the analysis. When analytes were present in more than one assay panel, the analyte with the largest dynamic working range was selected for downstream analysis to ensure maximum detection window, this resulted in data analysis on 74 unique markers.
Samples were coded, stored directly at 80C upon collection and not thawed before the first measurement, to minimize the loss of detection due to low level or instable analytes. Measurements were performed continuously when 76 samples were available from the study. Thus, longitudinal samples from an individual were not necessarily measured at the same time but throughout the study period. To ensure technical or inter-assay variation would not impact the analysis, a reference control sample, generated by pooling sera of 4 severely ill COVID-19 patients was included in all assays. The reference control was aliquoted to avoid repeated freezing and thawing cycles, maintaining the quality of the reference sample. After each acquisition all standard curves were checked for performance, outliers were deleted and curve fitting optimized (Regression Type: Logistic5PL with standard recovery between 70 and 130%). Datapoints out of range below the LLOQ were set to 0 pg/ml and above the ULOQ were set to maximum (200.000 pg/ml).
Evaluation of the reference control sample showed limited variation in the calculated concentrations for most analytes measured between assays (analytes: n=74), assays: 11 plates; median %CV=56.7% (Supplementary Fig.6A). Analytes with large %CV were mostly produced in high levels with standard curve characteristics of a steep slope and less sensitivity in the lower ranges, impacting the calculated concentration even with minimal fluorescence differences.
Normalized soluble marker levels were+1 (to avoid zero as possible outcome), log2-transformed. Data was visualised and modelled using R (R Core Team, 2023), RStudio (Posit Team, 2023) and packages (R Core Team, 2023), ggplot2 (Wickham H, 2016), ggdendro (De Vries A, Ripley BD, 2022), mixOmics (Rohart F, Gautier B, Singh A, Le Cao K-A, 2017), caret (Kuhn M, 2008), glmnet (Friedman J, Tibshirani R, Hastie T, 2010), psych (Revelle W, 2023), pROC (Robin X, Turck N et al., 2011), ggpubr (Kassambara A, 2023) and GLMMadaptive (Rizopoulos D, 2023). Statistical testing comparing medians of two groups was performed using MannWhitney U. In addition, for analysis of differential soluble analyte levels between disease outcome groups and healthy controls or between waves MannWhitney U was performed and p values were adjusted for multiple testing by BenjaminiHochberg correction. The heatmap is clustered row-wise based on complete-linkage hierarchical clustering. Supervised dimensionality reduction was carried out using partial least squares-discriminant analysis (PLS-DA), separating disease outcome groups and correlation analyses were performed using Spearmans Rank and BenjaminiHochberg correction was performed on p-values.
For identification of the best-classifying soluble analyte signatures-profiles to model outcome (fatal outcome yes or no), ICU admission, high maximum disease severity and time to disease recovery at the earliest available time point (closest to hospital admission), logistic regression with lasso regularisation was performed. Leave-one-out cross-validation and train-test split (training set=70%, test set=30% of dataset) were used to assess the performance of the trained regression models. The classifying performance was further assessed by evaluating sensitivity, specificity, receiver operating characteristic (ROC) curve and area under the ROC curve (AUC).
The progression over time of each soluble analyte was modelled using linear mixed effects models for the total cohort. Linear mixed models are an extension of the simple linear regression models when multiple measurements on the same subject are collected and use random effects to capture the serial correlation. Wave of sample collection was included as confounder into the model. To enhance comparability of the patients and stages of the disease process the number of days since disease onset was used for the timing of the samples within each patient. The non-linear progression over time was captured using natural cubic splines. To compare the progression between the different outcome groups the main effect and the interaction term between the days since disease onset and the outcome groups was added in the model. Similarly to explore association of each soluble analyte with the disease severity the main effect and the interaction term between the days since disease onset and disease severity was added in the model. For soluble analytes with values outside the limits of detection, the truncated normal distribution was used. Finally to capture the serial correlation random intercept and random slope terms were used.
Pathway analyses for exploratory purpose were performed in Ingenuity Pathway Analysis (Qiagen, Germany) with input of 40 significantly different analytes in a core analysis with the Ingenuity Knowledge database as reference set. The 40 significant analytes were selected based on the association between circulating levels and daily SCODA severity-score tested over time to identify different mean progressions during follow up.
The rest is here: