Prevalent and persistent new-onset autoantibodies in mild to severe COVID-19 – Nature Communications

Study cohorts

The COMMUNITY study is an ongoing longitudinal study which enrolled 2149 HCW and 118 COVID-19 patients admitted to Danderyd Hospital, Stockholm, Sweden, between April and May/June 2020 (HCW/patients)^29,30,31. The cohort is followed with blood sampling every four months and we assess anti-SARS-CoV-2 IgG using a multiplex bead array of SARS-CoV-2 proteins²⁸. In addition, HCW reported symptoms post-COVID-19 through electronic self-assessment forms at selected visits, including visits 3 to 5. The symptoms for self-assessment were anxiety, brain fatigue, impaired concentration, cough, depressed mood, diarrhea, dyspnea, dizziness, fatigue, fever, hair loss, headache, impaired hearing, impaired memory, ageusia, muscle/joint pain, nausea, numbness, anosmia, palpitations, skin disorders, sleep disturbance, and stomach ache. Symptom severity was graded as mild, moderate, or severe. Neuropsychiatric symptoms were defined as reporting one or several of anxiety, brain fatigue, impaired concentration, depressed mood, and impaired memory. The disease severity of the hospitalized patients ranged from mild to severe. Disease severity and comorbidities are presented in Supplementary Table 1. At the start of the study in April 2020, all patients currently hospitalized with confirmed COVID-19 and all employees were eligible for inclusion, with no exclusion criteria.

In the present study, we retrospectively considered visits 1 to 5, i.e., May 2020 to September 2021. We selected a subgroup of 478 HCW and 47 hospitalized COVID-19 patients based on previous serological results and reported symptoms post-COVID^29,30,31. HCW were selected in two steps. First, we selected HCW who were seronegative for anti-SARS-CoV-2 IgG at the first visit in Apr-May 2020, seroconverted prior to vaccination, and had participated in the visits immediately before and after seroconversion (n = 381). Second, we selected HCW who were seropositive at the first visit, had participated in all of the four subsequent visits, and had reported several symptoms post-COVID-19 (n = 97). Patients were selected based on participation in all of the first four follow-up visits (n = 47). In total, 525 individuals were selected, with an average of 4.8 samples per individual, yielding 2532 samples for autoantibody analysis. Demographics are presented in Supplementary Table 1.

Ethical approval for the COMMUNITY study was obtained from the Swedish Ethical Review Authority (Nos. 2020-01653 and 2021-04113). All healthcare workers left written informed consent for study participation. For the patients, oral informed consent was obtained instead of written informed consent due to the risk of contagion. In the case of incapacity, informed consent was obtained from patients’ next of kin. Oral informed consent was recorded in each patient’s medical record as well as in a separate file held by the responsible researcher. The use of oral consent was approved by the Swedish Ethical Review Authority.

The pre-pandemic healthy control group was university employees and students who had not received psychiatric care during their lifetime and were recruited in the timeframe April 2014–April 2017 as healthy controls for the Uppsala Psychiatric Patient Samples (UPP) cohort. Ethical approval for the collection of CSF was acquired by the Regional Ethical Review Board in Uppsala, Sweden (Nos. 2012/081 and 2014/148). Oral and written consent was obtained from all controls. All available controls with matched CSF and serum were selected for this study. Participants underwent a clinical health examination, including blood pressure and body mass index (BMI), and answered questionnaires on socio-demographics, medical history, heredity, and current medication, as well as an interview to evaluate any psychiatric symptoms. While no controls had ongoing symptoms that required specialist psychiatric care, the frequency of prior or ongoing mild and subclinical states of psychiatric conditions was higher than expected. All samples were acquired pre-pandemic. No other exclusion criteria were applied⁵⁵. Demographics are presented in Supplementary Table 1.

The neuro-COVID population has been described previously^56,57. Patients were prospectively included between April 2020 and June 2021. Inclusion was based on a positive PCR for SARS-CoV-2 in upper and/or lower airway samples, and at least one new-onset neurological symptom and presence of anti-SARS-CoV-2 IgG in serum, or typical COVID-19 symptoms in combination with pulmonary ground-glass opacities and consolidations on computed tomography scan of the thorax. Patients with previous central nervous system insults were excluded. One patient had a previous cerebrovascular insult. Clinical neurological evaluation was performed by an experienced neurologist. All lumbar punctures were performed as part of the clinical routine for neurological investigations. Neurological manifestations at the time of lumbar puncture included cranial nerve affection, central or peripheral paralysis, extrapyramidal, sensory symptoms, altered mental status including confusion, encephalopathy, and reduced level of consciousness. All patients were hospitalized. The disease severity at the lumbar puncture was moderate to severe. Disease severity, comorbidities, and other demographics are presented in Supplementary Table 1. Written informed consent was obtained from each patient, or next-of-kin if a patient was unable to give consent.

The collection and analysis of neuro-COVID and healthy control samples was approved by the Swedish Ethical Review Authority (No. 2017-043 with amendments 2019-00169, 2020-01623, 2020-02719, 2020-05730, 2021-01469, and 2020-01883; and No. 2022-00526-01). The Declaration of Helsinki and its subsequent revisions were followed.

SARS-CoV-2 serology

Serological classifications for sample selection were obtained from previous studies of the cohorts^29,30,31. For increased resolution in the upper ranges of the data, samples were re-analyzed at a higher dilution (1:5000 vs 1:50) while otherwise following the same procedure²⁸.

Planar arrays

Initial exploration of autoantibody repertoires was performed using two sets of in-house developed protein arrays. The Proteome-wide planar array contains 42,000 protein fragments from the Human Protein Atlas (proteinatlas.org) that represent 18,000 proteins and cover ~40% of the amino acid residues of the human proteome^33,58. The Secretome array contains 1522 full-length secreted or extracellular proteins from 1482 genes representing 58% of the human secretome, i.e., the proteins secreted by human cells³². The experimental procedure has been described previously¹⁴. In brief, plasma samples were pooled within the described groups and diluted 1:100 in before applying them to planar arrays. After incubation and washes, fluorescently labeled secondary detection antibodies were applied to the array for detection of autoantibody binding and detection of array microspots. Readout was performed in a microarray scanner, where the red channel readout corresponds to autoantibodies and the green channel readout to array microspots. Grid alignment and data acquisition from the scanned images was performed using GenePix Pro 5.1. Selected protein fragments are available on reasonable request.

Bead arrays

Bead arrays were used for the investigation of new-onset autoantibodies in the COMMUNITY cohort and for epitope mapping in the COMMUNITY and validation cohorts. Bead array construction and assays were performed as previously described^14,59. For an exploration of new-onset autoantibodies, the coupled antigens had been selected from planar arrays as described below. For epitope mapping, the coupled antigens were 93 custom-synthesized biotinylated peptides of 14 to 15 amino acid residues (GenScript Biotech). Assay readout was performed using Luminex FLEXMAP 3D® instruments with xPONENT® software (Luminex Corp.) and responses recorded as median fluorescent intensity (MFI), in arbitrary units.

Analysis of planar array data

Planar array data were processed as previously described³³.

Proteome-wide array data were analyzed per subarray. Spots that were flagged in image acquisition, that were smaller than 30 pixels, or that had a green channel signal lower than 4 SD above the mean local background were removed. Duplicate spots were deduplicated by selecting the spot with highest signal in the green channel. Finally, red channel data were Z-scored, and antigens with Z ≥ 12 were classified as reactive.

Secretome array data was analyzed per subarray. Microarray spots that were flagged in image acquisition, that were smaller than 40 pixels, or that had a green channel signal at or below the local background level were removed. Duplicate spots were filtered if their CV exceeded 50, and remaining spot pairs were deduplicated by taking the mean. Finally, local background was subtracted from red channel data and resulting values were Z-scored. Antigens with Z ≥ 8 were classified as reactive.

Z-scored planar array data were used for the selection of antigens for further investigation in the full HCW and hospitalized patient cohorts. First, antigens reactive in single samples on single arrays were selected. Second, lowered selection thresholds and detection in multiple samples were used to diversify the selection. The following selection criteria were used: protein fragments meeting Z ≥ 8 in multiple samples; full-length proteins meeting Z ≥ 4 in multiple samples; protein fragments and full-length proteins meeting Z ≥ 4 in single or multiple samples on both array types. Third, antigens from the literature and in-house studies were selected: protein fragments noted in both the literature and in-house studies; protein fragments meeting Z ≥ 8 and noted in either the literature or in-house studies; full-length proteins noted in multiple publications; full-length proteins meeting Z ≥ 2 and noted in the literature; handpicked protein fragments and full-length proteins noted in either the literature or in-house studies. Selected antigens and their matching selection criteria are listed in Supplementary Data 1.

Classification of new-onset autoantibodies in individuals with a seronegative baseline sample

New-onset autoantibodies were classified by applying Partitioning Around Medoids (PAM) clustering⁶⁰ to the bead array data of baseline seronegative individuals. The approach is documented in the companion R package abtract⁶¹.

The log₂ FC of each autoantibody trajectory was computed relative to the most recent seronegative sample. Individuals with autoantibody data at seroconversion and the samplings immediately before and after seroconversion were considered for PAM clustering (n = 374). Principal component analysis (PCA) was conducted on these data to identify and exclude outlying individuals with an increase in most autoantibody trajectories at seroconversion (n = 5, Supplementary Fig. 7a). Outliers were defined as PC1 ≥ mean + 3 × SD of PC1 (PC1 ≥8.85). To reduce the running time of the model, trajectories that never exceed background MFI levels were pre-classified as not new-onset. In each antigen class (protein fragments and full-length proteins), background MFI level was defined as robust Z-score = 3 in the seronegative samples from September 2020 and January 2021 from the 40 individuals that seroconverted in May 2021 (MFI_fragment ≥73, MFI_full-length ≥322; n filtered trajectories = 619894/910041 (68%)).

The PAM clustering was performed with the custom cosine × euclidean distance metric. To evaluate the parameter k (number of clusters), the PAM clustering was run with k ranging from 1 to 20, each with 10 random starts of the clustering algorithm. Using the silhouette method, the optimal cluster number was determined to be 7. Based on a median log₂ FC ≥2, clusters 5, 6, and 7 were classified as new-onset.

Classification of new-onset autoantibodies in individuals without any seronegative baseline sample

New-onset autoantibodies were assessed in individuals without seronegative baseline samples using multinomial linear regression (MNL)^62,63.

The MNL model was built using classifications and autoantibody trajectories of the 22 prevalent new-onset autoantibodies in the baseline seronegative individuals that had follow-up samples at 4 and 8 months after seroconversion. The MNL model specification was \({\bf{New}}{\boldsymbol{\hbox{-}}}{\bf{onset\; type}}\, \sim \,{\log }_{10}{{\bf{MFI}}}_{t\left(0\right)}+{{\bf{FC}}}_{t\left(4\right)}+{{\bf{FC}}}_{t\left(8\right)}\), where t is months from seroconversion, and FC is fold change relative to time of seroconversion (t(0)). The trajectories were split into training (0.8) and testing (0.2) sets using randomized stratified sampling on the trajectory outcome (target antigen and simplified new-onset classification; acute, delayed, or not new-onset autoantibody). The MNL model was trained using repeated stratified K-fold cross-validation (repeats = 10, K = 5, stratification on trajectory outcome). Grid search was used to optimize the penalizing decay parameter, which was set to 0 based on model accuracy.

Individuals that did not have seronegative baseline samples, but that did have autoantibody data at seroconversion and the two samplings immediately after seroconversion, were considered for MNL classification (n = 135). PCA was used to exclude outlying individuals with an increase in most autoantibody trajectories in the considered time points (n = 2, Supplementary Fig. 7b). Outliers were defined as PC1 ≥ mean + 3 × SD of PC1 (PC1 ≥2.1). Trajectories of the 22 prevalent new-onset autoantibodies in the remaining 133 individuals were classified using the MNL model (n trajectories = 2926).

Annotation of antigen location

Antigen location was determined based on Cellular Component Gene Ontology (GO) terms for each corresponding gene^64,65. The GO terms were simplified using the Generic GO subset⁶⁶ and further subdivided into broad location categories as detailed in Supplementary Table 3. For antigens with multiple locations, the annotation was selected in the following order: Extracellular > Plasma membrane > Nuclear > Intracellular. Antigens may lack annotation of location if no parent terms are included in the GO subset.

Association of new-onset autoantibodies with symptoms post-COVID-19

Association of new-onset autoantibodies with symptoms post-COVID-19 was performed using data from HCW displaying one or more of the 22 prevalent new-onset autoantibodies, or none (n = 403). Among these, neuropsychiatric symptoms were defined as reporting anxiety, brain fatigue, difficulties concentrating, depressive disorders, or impaired memory that lasted for at least 2 months after COVID-19 (n = 384). The highest reported severity was used (severe = 20 (5%); moderate = 58 (15%); mild = 23 (6%); no neuropsychiatric sx = 283 (74%)). For other symptoms post-COVID-19, reported mild symptoms were excluded (Supplementary Table 4).

Identification of main epitopes

The peptide bead array was used for epitope mapping. Data were acquired from 142 baseline samples and 150 samples after autoantibody onset in the 142 individuals who had one or more new-onset autoantibodies whose antigens were represented on the peptide bead array. (Eight individuals had both acute and delayed new-onset autoantibodies, increasing the number of samples after autoantibody onset.) The FC was computed relative to baseline, both for beads with coupled peptides as well as the negative control bead with coupled biotin. To adjust for individual background levels, the FC of biotin was subtracted from the FC of the peptides, and the median was set to 1. Any resulting negative values (n = 2) were imputed with the smallest positive value present (0.07).

Main epitopes were defined in individuals having the corresponding new-onset autoantibody by taking the mean FC of each peptide. Peptides, where the mean exceeded the background FC cutoff, were defined as main epitopes that were common across individuals. The background FC cutoff was defined as mean + 5 × SD of the FC of the negative biotin control (FC ≥ 2.52).

Sequence alignment

Sequence alignment was performed using NCBI BLASTP adjusted for short input sequences: E 

Statistical analysis

Correlations were performed using Spearman correlation. Group differences of continuous variables were investigated using the Mann–Whitney U-test. Binary variables were investigated using logistic regression, categorical variables using multinomial logistic regression, and ordinal variables using proportional odds logistic regression. Regression models used correction for age and sex in group comparisons. False discovery rate correction was performed using the Benjamini–Hochberg procedure, and resulting values are reported as q values. All statistical tests were two-sided. All measurements were taken from distinct samples.

All data analysis and statistical analysis was performed in RStudio with R 4.1.2 and 4.2.1 and R packages tidyverse, lubridate, rlang, scales, knitr, pander, httr2, readxl, rstatix, proxy, cluster, caret, MLmetrics, MASS, broom, Omixer, ggpubr, ragg, patchwork, cowplot, GGally, ggsignif, ggfortify, ggdist, ggbeeswarm, ggrepel, ComplexHeatmap, and circlize.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.