Type I interferon autoantibodies are associated with systemic immune alterations in patients with COVID-19.

Neutralizing autoantibodies against type I interferons (IFNs) have been found in some patients with critical coronavirus disease 2019 (COVID-19), the disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). However, the prevalence of these antibodies, their longitudinal dynamics across the disease severity scale, and their functional effects on circulating leukocytes remain unknown. Here, in 284 patients with COVID-19, we found type I IFN–specific autoantibodies in peripheral blood samples from 19% of patients with critical disease and 6% of patients with severe disease. We found no type I IFN autoantibodies in individuals with moderate disease. Longitudinal profiling of over 600,000 peripheral blood mononuclear cells using multiplexed single-cell epitope and transcriptome sequencing from 54 patients with COVID-19 and 26 non–COVID-19 controls revealed a lack of type I IFN–stimulated gene (ISG-I) responses in myeloid cells from patients with critical disease. This was especially evident in dendritic cell populations isolated from patients with critical disease producing type I IFN–specific autoantibodies. Moreover, we found elevated expression of the inhibitory receptor leukocyte-associated immunoglobulin-like receptor 1 (LAIR1) on the surface of monocytes isolated from patients with critical disease early in the disease course. LAIR1 expression is inversely correlated with ISG-I expression response in patients with COVID-19 but is not expressed in healthy controls. The deficient ISG-I response observed in patients with critical COVID-19 with and without type I IFN–specific autoantibodies supports a unifying model for disease pathogenesis involving ISG-I suppression through convergent mechanisms.

INTRODUCTION

The coronavirus disease 2019 (COVID-19) pandemic has led to the infection of at least 192 million individuals worldwide and more than 4.1 million deaths. A perplexing aspect of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pathogenesis is the extreme clinical heterogeneity of infected individuals, with about 15% of symptomatic patients and less than 10% of infected individuals presenting with severe forms of the disease, as defined by dyspnea, pulmonary infiltrates on lung imaging, and low blood oxygen saturation (–). Overall, 26.8% of patients who are hospitalized develop critical disease defined as category 7 on the National Institutes of Health (NIH) ordinal scale, requiring mechanical ventilation (). These patients are at the greatest risk for poor outcome and place the largest burden on the health care system. Despite increasing vaccine availability, some vulnerable individuals may develop critical disease before and even perhaps despite vaccination, especially in the context of emerging highly transmissible, more virulent, and antigenically distinct variants of SARS-CoV-2 isolates (–). Thus, there is a need to disentangle the immunological consequences of SARS-CoV-2 infection and the underlying immunological causes of critical COVID-19 to stratify patients early in their disease course and to target treatments accordingly.

Evidence is emerging that genetic and immunological features that predate SARS-CoV-2 infection could play an unexpected pathogenic role in severe disease (). Among patients with critical COVID-19, these features include inborn errors of type I interferon (IFN)–mediated immunity (, ) and the production of autoantibodies against type I IFNs (, ). These autoantibodies, which seldom occur in healthy controls (with frequencies less than 0.3%) and have not been found in asymptomatically infected individuals, are observed in at least 10% of patients with critical COVID-19 (, –). The causal effects of autoantibodies against type I IFNs on COVID-19 severity has been supported by their documentation before infection and their frequent occurrence in patients with genetic disorders, such as autoimmune polyglandular syndrome type 1 (APS-1) (–).

However, it remains to be determined whether autoantibodies to type I IFNs occur in patients with COVID-19 who do not require mechanical ventilation, whether they fluctuate longitudinally during the disease course, and what their consequences are on the composition and phenotypes of circulating leukocyte subsets. Furthermore, few studies have examined circulating leukocytes over the course of SARS-CoV-2 infection (, ) or have compared with patients presenting with similar respiratory manifestations requiring hospitalization due to other causes (). Insights into how the natural innate and adaptive immunity longitudinally evolve in response to SARS-CoV-2 infection, in both anti–type I IFN autoantibody-positive and autoantibody-negative cases, may enable the early identification of patients who are likely to develop life-threatening COVID-19 and the discovery of mechanisms that can be targeted by therapy.

RESULTS

Anti–IFN-α2 antibodies are reproducibly detected in patients with COVID-19 but seldomly detected in the general population

Given the recent description of neutralizing type I IFN autoantibodies in greater than 10% of critical COVID-19 cases, we sought to determine the frequency of these antibodies in San Francisco in a total of more than 4,500 individuals divided over the following: (i) SARS-CoV-2–positive individuals that span the NIH COVID-19 severity scale (), (ii) a largely asymptomatic community population, and (iii) convalescent serum samples from patients previously infected with SARS-CoV-2. We first determined the frequency of autoantibodies to the type I IFN, IFN-α2, in 284 patients with confirmed SARS-CoV-2 infection using a radioligand binding assay (RLBA). These patients were categorized using the NIH ordinal scale, with those scoring between 1 and 4 classified as moderate, those scoring 5 or 6 as severe, and those scoring 7 as critical (Table 1). Of the 284 patients, 53 were enrolled in the COVID-19 Multi-Phenotyping for Effective Therapies (COMET) cohort (Table 2 and data file S1) and 231 were enrolled from Zuckerberg San Francisco General Hospital and Trauma Center (ZSFG) (Table 3 and data file S2). We also tested 4 samples from individuals negative for COVID-19 with an APS-1 diagnosis as positive controls and 14 samples from individuals negative for COVID-19 without an APS-1 diagnosis as negative controls (Fig. 1A). The 284 patients ranged in age from 0 to 90+ years, were 69% male, had at least one positive SARS-CoV-2 polymerase chain reaction (PCR) test, and varied in disease severity (Tables 2 and 3 and data files S1 and S2). We found the prevalence of anti–IFN-α2 autoantibodies to be 5 of 26 (19%) in critical, 6 of 102 (6%) in severe, and 0 of 156 in moderate disease (Fig. 1A; see fig. S1 for the neutralizing capacity of these antibodies in the anti–IFN-α2–positive patients; P < 0.01, Fisher’s exact test comparing severe-to-critical to moderate cases). The positive patients were aged 28 to 72 years (mean = 55.7, std = 12.2), and 9 of 11 (82%) were male. The prevalence of anti–IFN-α2 in critical COVID-19 and the tendency of positive patients to be of the male sex and of advanced age are consistent with previously published descriptions ().

Table 1.

Categorization of moderate, severe, and critical patients in the COMET and ZSFG cohorts.

ECMO, extracorporeal membrane oxygenation; n/a, not applicable.

NIH ordinal scale	WHO criteria	Categorization: COMET cohort	Categorization: ZSFG cohort
1	Not hospitalized, no limitations onactivities	n/a	Moderate
2	Not hospitalized, limitation onactivities or requiring home oxygen	n/a	Moderate
3	Hospitalized, not requiringsupplemental oxygen—no longerrequires ongoing medical care	Moderate	Moderate
4	Hospitalized, not requiringsupplemental oxygen—requiringongoing medical care	Moderate	Moderate
5	Hospitalized, requiring supplementaloxygen	Severe	Severe (oxygenation information notavailable)
6	Hospitalized, on noninvasiveventilation or high-flow oxygendevices	Severe
7	Hospitalized, on mechanicalventilation or ECMO	Critical	Critical

Table 2.

Demographics of the COMET cohort.

Demographics and clinical characteristics of patients from the COMET cohort, broken down by anti–IFN-α2+ COVID-19+, anti–IFN-α2− COVID-19+, and COVID-19− patients. RLBA, radioligand binding assay.

Characteristic	COVID-19+ anti–IFN-α2− (n = 50)*	COVID-19+ anti–IFN-α2+ (RLBA value > 6 SDs above healthy controls) (n = 4)	COVID-19− (n = 15)†
Mean age (SD)	58.1 (16.8)	56.3 (11.8)	61.5 (14.2)
Sex
Male	37	4	7
Female	13	0	8
Race
White or European	13	0	7
Asian	9	1	3
Native Hawaiian/Other PacificIslander	3	0	0
Black/African-American	1	0	1
Other/Multiple races	22	3	2
Unknown	2	0	2
Ethnicity
Hispanic/Latino	20	3	1
Non-Hispanic/Latino	30	1	11
Unknown	0	0	3
Duration hospitalization (SD)	18.3 (15.6)	22.3 (8.4)	7.5 (5.1)
Died
Yes	6	1	4
No	44	3	11

*One of 50 not tested for anti–IFN-α2.

†One of 15 not tested for anti–IFN-α2.

Table 3.

Demographics of the ZSFG cohort.

Demographics and clinical characteristics of patients from the ZSFG cohort, including comparison across anti–IFN-α2–positive and anti–IFN-α2–negative patients. Significance values were determined using Fisher’s exact test (FE), except in the case of continuous distributions (length of stay), where a Kolmogorov-Smirnov (KS) test was used.

Characteristic	COVID-19+ anti–IFN-α2− (n = 224)	COVID-19+ anti– IFN-α2+ (RLBA value > 6 SDs above healthy controls) (n = 7)	P value
Mean age (SD)	50.3 (15.9)	55.3 (13.3)
Sex
Male	144	6	0.46 (FE)
Female	78	1
Unknown	2	0
Race
Black or African-American (8%)	19	0	1 (FE)
White (6%)	12	0	1 (FE)
Asian (6%)	12	1	0.35 (FE)
Other, includingHispanic or Latinoethnicity (78%)	178	5	1 (FE)
Not reported (2%)	3	1	0.13 (FE)
Ethnicity
Hispanic or Latino(75%)	170	5	0.68 (FE)
Non-Hispanic orLatino (25%)	54	2
Hospital admission	103	7	0.15 (FE)
Noninvasiveoxygenation ormore	79	6	0.12 (FE)
ICU	48	6	0.02 (FE)
Mechanicalventilation	31	4	0.04 (FE)
Pressors	25	5	0.006 (FE)
Death	7	2	0.04 (FE)
Mean length of stay (SD)	5.6 (10.7)	28 (30.0)	0.0005 (KS)
Men	5.5 (9.6)	30.2 (32.3)	0.003 (KS)
Women	5.4 (12.3)	15 (n = 1)	0.3 (KS)

Fig. 1.

Detection of anti–IFN-α2 antibodies is positively associated with COVID-19 disease severity.

(A) Anti–IFN-α2 antibody (Ab) index values (y axis) for 4 APS-1 patients, 26 critical COVID-19+ (C19+) cases, 102 severe COVID-19+ cases, 156 moderate COVID-19+ cases, and 14 moderate-severe COVID-19− (C19−) cases, separated by disease severity and colored by hospitalization status. Positive samples were tested for neutralization against IFN-α2 and IFN-ω (see fig. S1), with arrows indicating those samples with partial or full neutralization ability. The dotted line indicates 6 SDs above healthy control mean. (B) Distribution of the anti–IFN-α2 index values across 4041 individuals in a community cohort from the San Francisco Mission District. (C) Anti–IFN-α2 index values (y axis) for five additional APS-1 patients and 175 CCP donors from the Vitalant Blood Center. (D) Timeline of blood draws for participants in the COVID-19 Multi-Phenotyping for Effective Therapies (COMET) cohort. Disease status, severity, and gender breakdowns are also shown. (E) Anti–IFN-α2 Ab index values over days since first hospitalization for 53 hospitalized COVID-19+ and 14 COVID-19− COMET samples. For 2 of 69 COMET samples, anti–IFN-α2 titers were not assessed.

We next examined a community cohort collected during a study of SARS-CoV-2 transmission in San Francisco (Fig. 1B, Table 4, and data file S3) (). The cohort consisted of 4,041 individuals aged 4 to 90 years of European (36%), Hispanic/LatinX (33%), Asian/Pacific Islander (9%), Black/African American (2%), and other or unknown (20%) ancestry. In this cohort, a total of 13 anti–IFN-α2–positive individuals (0.32%) were identified. Of these, five were male, six were female, and two were of unknown gender. Positive samples were identified across all represented ethnic groups. None of the participants who were confirmed positive for past or present SARS-CoV-2 infection (117 of 3,851 by serology and 64 of 3,758 by PCR) were positive for anti–IFN-α2 antibodies, and all were ambulatory or asymptomatic at the time of testing. These data are consistent with the previously reported absence of autoantibodies in patients with ambulatory COVID-19 (). Our results also confirm the low frequency of anti–IFN-α2 antibodies in individuals independent of, and likely before, infection with SARS-CoV-2. In addition to assessing the presence of anti–IFN-α2 antibodies in San Francisco community cohorts, we analyzed aliquots of convalescent plasma (CCP) from a central blood bank supplier, encompassing 175 unique plasma donors who had recovered from SARS-CoV-2. Compared with five additional individuals diagnosed with APS-1, we found that none of the donors tested positive for anti–IFN-α2 autoantibodies (Fig. 1C). This latter cohort suggests that these potentially harmful autoantibodies are rare or absent in the supply from convalescent donors.

Table 4.

Demographics of the community cohort.

Demographics of patients from the community cohort, broken down by anti–IFN-α2–positive and anti–IFN-α2–negative individuals. Significance values were determined using Fisher’s exact test.

Characteristics	Anti–IFN-α2− (n = 4028)	Anti–IFN-α2+ (z score > 3.3, P = 0.005) (n = 13)	P value
Mean age (SD)	41.0 (15.9)	51.5 (20.3)
Sex
Male	1811	5	0.77
Female	1731	6
Unknown	486	2
Race/ethnicity
Hispanic or Latino (33.6%)	1354	4	1
White or European (36.2%)	1458	2	0.15
Asian or Pacific Islander (8.6%)	347	1	1
Two or more races (2.8%)	113	2	0.05
Black or African-American (2.0%)	79	1	0.23
Other (3.6%)	143	1	0.38
Unknown (13.2%)	534	2	0.69
Antibody test result
Positive	117	0
Not detected	3723	11
Unknown	188	2
PCR test result
Positive	64	0
Not detected	3684	10
Unknown	280	3
Antibody + PCR testing
Positive for one or both tests	154	0
Negative or unknown for both	3874	13

Single-cell epitope and RNA sequencing of peripheral blood mononuclear cells enables assessment of systemic immune associations to COVID-19 and the presence of type I IFN autoantibodies

We next sought to identify the associations of COVID-19 status and the presence of anti–IFN-α2 autoantibodies with the composition, transcript abundance, and surface protein abundance of circulating leukocytes. For this, we leveraged the COMET cohort where peripheral blood mononuclear cells (PBMCs) and serum were longitudinally collected from 69 hospitalized patients presenting with COVID-19 symptoms (67 of 69 patients overlapped with those for whom anti–IFN-α2 autoantibodies were assessed; Fig. 1A), of whom 54 were positive (COVID-19+) and 15 were negative (COVID-19−) for SARS-CoV-2, in addition to 11 healthy controls (Fig. 1D). Of the COVID-19+ cases, 18 presented with moderate disease, 17 with severe disease, and 19 with critical disease according to the NIH severity scale () at the time of hospitalization (Tables 1 and 3 and data files S1 and S4). For 8 of 54 COVID-19+ patients, the severity changed over the course of hospitalization, of whom 6 improved and 2 worsened. The studied hospitalized patients were ethnically diverse, skewed older than the general population (mean = 59, range = 25 to 90), and were predominantly male (47 men and 22 women) (Fig. 1D). Although all COVID-19− cases presenting with symptoms concerning for COVID-19 tested negative for SARS-CoV-2, many were infected with common respiratory pathogens confirmed by metagenomic sequencing (data file S1). Among the 67 of 69 COMET patients assessed for anti–IFN-α2 autoantibodies, 4 of 19 (21%) of the critical COVID-19 cases, none of the moderate to severe cases, and none of the COVID-19− cases were positive for anti–IFN-α2 antibodies (Fig. 1A). All four cases had anti–IFN-α2 antibodies at the earliest time of sampling, and the concentration of anti–IFN-α2 antibodies remained stable for three of four cases across their disease course (Fig. 1E).

To profile the circulating immune response, we collected about 200 PBMC samples from up to four longitudinal time points: 0, 4, 7, and 14 days since hospitalization. Multiplexed single-cell epitope and transcriptome sequencing (muxCITE-seq) was performed across nine pools of genetically distinct samples to simultaneously measure mRNA abundances transcriptome-wide and surface protein abundances of 189 markers from the same cell (fig. S2 and data file S5). A total of 971,550 cell-containing droplets were sequenced, and 600,929 cells remained in the final dataset after quality control and removal of doublets, platelets, and red blood cells. Genetic demultiplexing using Freemuxlet resulted in an average of 3020 cells per sample (fig. S3A).

Critical COVID-19 is characterized by increased frequencies of plasmablasts and classical monocytes

We compared the frequencies of 11 cell types between COVID-19+ cases, COVID-19− cases, and healthy controls, as well as within COVID-19+ cases separated by severity. The assessed cell types were defined by their transcriptomic profiles and included plasmablasts (PBs), B cells (B), CD4+, CD8+, and γδ T cells (CD4T, CD8T, and γδT), natural killer cells (NKs), conventional and plasmacytoid dendritic cells (cDC and pDC), classical and nonclassical monocytes (cM and ncM), and hematopoietic progenitor cells (Progens) (Fig. 2A). We first confirmed that muxCITE-seq–derived estimates of lymphocyte and monocyte frequencies were correlated with complete blood count measurements reported in the electronic health record from the same donor within 2 days (Pearson Rmonocyte = 0.59, P = 7.2 × 10−105; Pearson Rlymphocyte = 0.57, P = 8.2 × 10−55; fig. S3B). Qualitatively, COVID-19+ cases exhibited shifts in the uniform manifold approximation and projection (UMAP) space of circulating leukocytes, particularly of myeloid cells, that were not confounded by processing batch and pool (Fig. 2B and fig. S3C). Comparing critical COVID-19+ cases to healthy controls, we observed statistically significant changes in frequencies for every cell type, including prominent increases in the frequencies of cM, B, and PB (cM: median change +10.0%, differential proportion analysis permutation P < 10−5; B: +2.7%, P = 2.1 × 10−3; PB: +2.1%, P < 10−5) and decreases in the frequencies of CD8T, γδT, NK, cDC, and pDC (CD8T: −15.4%, P < 10−5; γδT: −3.9%, P < 10−5; NK: −4.7%, P < 10−5; cDC: −2.2%, P < 10−5; pDC: −0.6%, P = 3.0 × 10−4). The frequencies of CD8T, γδT, pDC, and PB in moderate and severe cases were between those observed in critical cases and healthy controls (Fig. 2C, fig. S3D, and data file S6). The frequency of CD8Ts was even lower, and the frequency of cMs was even higher in critical COVID-19+ cases with detectable anti–IFN-α2 antibodies than those without (CD8T: −5.8%, P < 10−5; cMs: +8.3%, P = 0.034) (Fig. 2C and data file S6). Changes in frequencies of B, PB, and ncMs were also significantly different between critical COVID-19+ cases and COVID-19− hospitalized patients (B: +10.4%, P < 10−5; PB: +2.1%, P < 10−5; ncM: −2.4%, P < 10−5), suggesting that these effects are not likely explained by hospitalization in general (Fig. 2C and data file S6). For the 14 COVID-19+ donors for whom all four time points were available, we observed decreases in the frequencies of B and PB cells over time (median change D0 versus D14; B: −3.7%, P = 6.0 × 10−5; PB −0.8%, P < 10−5) and increases in the frequencies of cM and ncMs (D0 versus D14: cM +5.7%, P < 10−5; ncM +2.5%, P < 10−5) for both days since hospitalization and days since onset of first symptoms (Fig. 2D; fig. S3, E and F; and data file S6). These longitudinal changes normalized toward frequencies observed in healthy controls. The exception was the frequency of cMs, which appeared to further increase from those observed in healthy controls. Previously, the frequency of PBs had been observed to correlate with COVID-19 disease severity (, ) and to diminish upon recovery (). We observed a positive correlation between viral RNA abundances from tracheal aspirate samples, collected from mostly critical cases, and PB frequency (in critical cases: Pearson R = 0.47, Padjusted = 0.0083) (Fig. 2E); both the PB frequencies (Fig. 2D) and viral RNA abundances decreased over time (fig. S4A). As expected from the high correlation between viral RNA abundances measured in tracheal aspirates and nasal swabs from the same donors (fig. S4B; Pearson R = 0.94), we replicated the high correlation between viral RNA and PB frequency in nasal swab samples of critical cases (Pearson R = 0.58, Padjusted = 0.023) (Fig. 2E). However, in moderate and severe cases, no association was found (Pearson R = 0.08, Padjusted = 0.76) (Fig. 2E). On the contrary, for none of the other cell types, a correlation between cell type frequency and viral titers was identified, independent of disease severity (fig. S4, C and D). Together, these results suggest coordinated dynamic changes of host humoral immunity and viral load over the course of hospitalization in critical cases but do not implicate a specific cause. Nevertheless, a recent study by Stephenson et al. () confirmed that viral RNA abundance is an independent contributing factor to cell type proportion changes observed in patients with COVID-19. Overall, these analyses revealed shifts in cell type composition specific to COVID-19, between patients of varying disease severity, and over time. The general comparable composition of circulating leukocytes among patients with COVID-19+ of varying severity with and without anti-IFN autoantibodies (fig. S3D) suggests the presence of a broader, conserved mechanism underlying critical disease, such as additional IFN-related pathology, particularly in the autoantibody-negative patients.

Fig. 2.

Shifts in circulating leukocyte composition are observed in samples isolated from individuals with critical COVID-19.

(A) Marker genes are shown for each of the 11 cell types identified including CD4+ T cells (CD4T), CD8+ T cells (CD8T), γδ T cells (γδT), natural killer cells (NK), B cells (B), plasmablasts (PBs), classical monocytes (cM), nonclassical monocytes (ncM), conventional DCs (cDC), plasmacytoid DCs (pDC), and CD34+ hematopoietic progenitors (Progen). (B) UMAP projections of PBMCs are shown for donors separated by COVID-19 status and severity. Cells are colored by type. (C) Box plots of cell type percentages (y axis) by COVID-19 status and severity on day of hospital admission (D0). Each dot represents the percentage of a specific cell type per donor. Statistical comparisons between cells from all COVID-19+ critical donors (including the anti–IFN-α2 autoantibody donors) and healthy controls, COVID-19− donors, or combined COVID-19+ moderate-severe donors are shown. (D) Box plots show the percentages of PB and B (y axis) in patients with COVID-19 over days 0, 4, 7, and 14 since hospitalization (D0, D4, D7, and D14). Data in (C) and (D) are presented as median and 25th and 75th percentiles. Statistically significant differences are presented using Holm’s multiple-testing corrected, permutation-based P values: *P < 0.05, **P < 0.01, and ***P < 0.001; ns, not significant. (E) Scatter plot of normalized SARS-CoV-2 viral RNA abundances (x axis) versus percentage of PB (y axis) was quantified in donor-matched single-cell PBMC data. RNA was measured by RT-qPCR in tracheal aspirates (left) or nasal swabs (right) as inverse ΔCT (inverse Ct difference between viral E gene and human RNase P gene). R, Pearson correlation in critical (RCr) or moderate-severe (RMS) cases.

Critical COVID-19 is marked by deficient type I IFN–stimulated gene expression in myeloid cells early in disease

To further characterize cell type intrinsic differences associated with COVID-19 and the presence of anti–IFN-α2 antibodies, we compared mRNA and surface protein abundances between COVID-19+ cases, COVID-19− cases, and healthy controls for each cell type. We identified 161 genes [false discovery rate (FDR) < 0.05, log2FC > 1] whose transcripts were differentially up-regulated between COVID-19+ cases at day 0 and healthy controls in at least one cell type (Fig. 3A, fig. S5A, and data file S7). K-means clustering of the 161 differentially expressed genes aggregated for each of the 11 cell types at day 0 identified five clusters, including a cluster (cluster 1) of genes enriched for type I IFN signaling and viral response primarily expressed in myeloid cells [gene set enrichment analysis (GSEA): type I IFN signaling pathway, permutation P < 10−5], a cluster (cluster 2) enriched for neutrophil degranulation (GSEA: neutrophil degranulation, P < 10−5), a cluster (cluster 3) of immunoglobulins (Igs) and PB activation markers, and a cluster (cluster 4) enriched for complement activation in ncMs (GSEA: complement activation, P = 0.026) (data file S8). Given the heterogeneous expression of the IFN signaling cluster (cluster 1) within COVID-19 samples, we further compared the expression of type I–specific (ISG-I) and type II–specific IFN-stimulated genes (ISG-II) between COVID-19+ cases and healthy controls and within COVID-19+ cases of varying severity. To differentiate ISG-I and ISG-II, we compared healthy donor PBMCs stimulated with recombinant IFN-β or IFN-γ from an independent single-cell RNA sequencing (scRNA-seq) dataset to identify genes specifically up-regulated by either IFN (Fig. 3B). In myeloid cells (cM, ncM, cDC, and pDC), the average expression of ISG-I, and to a lesser extent ISG-II, in critical cases on day 0 of hospitalization was significantly lower compared to moderate and severe cases [ISG-I: log2FC = −0.51 to −0.82, P = 8.6 × 10−4 to 6.9 × 10−3; ISG-II (cM only): log2FC = −0.46, P = 3.1 × 10−3] (Fig. 3C and fig. S5, B and C). Furthermore, in DCs (cDC and pDC), ISG-I expression for the four critical COVID-19+ cases with anti–IFN-α2 autoantibodies was significantly lower than that for critical cases without autoantibodies [ISG-I: log2FC = −0.75 to −1.05, P = 2.0 × 10−2 to 2.8 × 10−2; ISG-II (pDC only): log2FC = −0.92, P = 2.1 × 10−3; Fig. 3C]. Through the disease course, average expression of ISG-I in moderate cases was high at the time of hospitalization but quickly diminished, whereas in critical cases, especially those with anti–IFN-α2 autoantibodies, average expression of ISG-Is remained low [linear mixed model (LMM); moderate: slope = −4.8 × 10−3, P < 0.05; other groups, P > 0.05; Fig. 3D]. This deficient ISG-I response in critical versus moderate-severe cases independent of days since symptom onset is consistent with previous literature (). These findings suggest that there may be a shared causal mechanism of critical disease in patients with and without autoantibodies to type I IFNs.

Fig. 3.

The degree of activation of the ISG-I response differs among the COVID-19 severity scale.

(A) The heatmap shows 161 differentially expressed genes at day 0 [FDR < 0.01, |log(fold change)| > 1] in at least 1 of 11 cell types. CD4+ T cells (CD4T), CD8+ T cells (CD8T), NK cells, B cells (B), PB, cM, ncM, and cDC are shown. Each row represents a gene, and each column is the average expression of the genes in a particular sample across all cells of a specific type. Samples are grouped by both case control status and COVID-19+ severity. Expression is row standardized. Genes are grouped by cluster, with the enriched clusters annotated. (B) Matrix plots showing the expression of shared, type I–specific, and type II–specific ISGs in cM. The shared and specific ISGs were defined using an orthogonal scRNA-seq dataset containing PBMCs stimulated with IFN-β and IFN-γ or left unstimulated (Ctrl) (left). The expression of these same gene set was then plotted in the COMET cohort separated by case control status and disease severity (right). (C) ISG-I and ISG-II scores (y axis) at day 0 across four myeloid cell types and pseudobulk of all other cell types, separated by case control status and disease severity (x axis). Box plots show median and 25th and 75th percentile. (D) ISG-I scores (y axis) over the course of disease for healthy controls and COVID-19− and COVID-19+ cases in cMs. COVID-19+ cases are separated by severity and the presence of anti–IFN-α2 antibodies. Inset shows significant results of linear mixed model (LMM) testing for changes over time. *P < 0.05, **P < 0.01, and ***P < 0.001.

ISG-I deficiency is associated with elevated surface expression of leukocyte-associated Ig-like receptor 1 in myeloid cells

We next sought to identify changes in the expression of surface proteins in patients with COVID-19 that may be correlated with ISG-I expression. The correlation of surface protein and transcript abundance varied across the 189 targeted genes, with lineage-specific surface markers exhibiting the highest correlation (fig. S6A). Comparing COVID-19+ cases at day 0 of hospitalization to healthy controls for each cell type separately, we identified 5 of 189 differentially expressed surface proteins in cMs and an additional 14 of 189 in other cell types (|log2FC| > 0.5, FDR < 0.05; Fig. 4A and fig. S6B). Of the five proteins differentially expressed in cMs, four [transferrin receptor (TFRC), sialic acid binding Ig-like lectin 1 (SIGLEC1), Fc fragment of IgG receptor 1a (FCGR1A), and leukocyte-associated Ig-like receptor 1 (LAIR1)] were higher expressed in COVID-19+ cases. SIGLEC1 is a known up-regulated ISG () whose pattern of surface expression is consistent with the expression of ISG-I present in moderate and severe cases (Figs. 3C and 4B). In contrast, LAIR1, a collagen receptor (, ) and a type 1 Ig containing two immunoreceptor tyrosine-based inhibitory motifs (), was most differentially up-regulated in myeloid cell types of critical COVID-19+ cases compared to healthy controls (cM, ncM, and cDC: log2FC = 0.51 to 0.97, Padjusted = 4.8 × 10−3 to 3.9 × 10−5) and compared to moderate-severe COVID-19+ cases (cM and ncM: log2FC = 0.47 to 0.50, Padjusted = 1.2 × 10−2 to 7.9 × 10−3; Fig. 4B and fig. S6B). LAIR1 expression in cMs from critical COVID-19+ cases was high early in the disease course and significantly decreased over time in those patients (LMM; COVID-19+ critical: slope = −1.3 × 10−4, P < 0.05; COVID-19+ critical anti-IFN: slope = −1.5 × 10−4, P < 0.05; Fig. 4C and fig. S6C). LAIR1 transcripts were also differentially expressed by cM, ncM, cDC, and NK cells in COVID-19+ cases when compared to healthy controls (log2FC = 0.65 to 1.06, FDR = 5.0 × 10−6 to 3.0 × 10−3; other cell types: log2FC = −1.2 × 10−2 to 1.6 × 10−2, FDR = 0.64 to 0.93), and the transcript and protein abundances were correlated (mean across all cell types: Pearson R = 0.89, P = 2.3 × 10−4). In COVID-19+ cases at day 0, the expression of surface LAIR1 was inversely correlated with expression of ISG-Is in cMs (Pearson R = −0.47, P < 0.01) and ncMs (Pearson R = −0.41, P < 0.01) (Fig. 4, D and E, and fig. S6D). Furthermore, LAIR1 surface protein expression for 40 samples was significantly correlated with those obtained from flow cytometry (Pearson R = 0.43, P = 5.5 × 10−3) (fig. S6E). Unlike known IFN-repressed surface proteins [such as CD244 ()] that were also inversely correlated with the ISG-I score (Fig. 4D), LAIR1 is not expressed in healthy samples, suggesting that it is not an IFN-repressed gene. These results demonstrate high surface expression of LAIR1 in myeloid cells to be a biomarker predictive of deficient type I–specific IFN response.

Fig. 4.

Surface protein abundance changes are observed in leukocyte subsets of patients with critical COVID-19.

(A) Volcano plot of log fold change between COVID-19+ and healthy controls (x axis) versus −log10(P value) (y axis) in cM cells. Proteins that are significantly (FDR < 0.05) differentially abundant and have a log2(fold change) > 0.5 are highlighted. (B) Normalized (Norm.) LAIR1 and SIGLEC1 surface expression (exp.; y axis) at day 0 across eight cell types separated by case control status, severity, and presence of anti–IFN-α2 antibodies (x axis). Box plots show median and 25th and 75th percentile. (C) Normalized LAIR1 surface expression (y axis) in cMs over the course of disease for healthy controls, COVID-19− controls, and COVID-19+ cases. COVID-19− controls and COVID-19+ cases are separate by severity and the presence of anti–IFN-α2 antibodies. Insets show significant results of LMM testing for changes over time. (D) The bar plot shows correlation between the surface expression of 52 statistically significantly proteins to ISG-I score in cMs at day 0. Proteins are colored by their log2(fold change) expression between COVID-19+ cases and healthy controls. Red indicates higher expression in COVID-19+ cases. Blue indicates lower expression in COVID-19+ cases. (E) Scatter plot of normalized LAIR1 expression (y axis) versus the ISG-I score (x axis) for COVID-19+ cases colored by severity and anti–IFN-α2 status. *P < 0.05, **P < 0.01, and ***P < 0.001.

DISCUSSION

The marked clinical heterogeneity over the course of SARS-CoV-2 infection, ranging from asymptomatic to lethal, is a key observation and defining feature of the COVID-19 pandemic. It is important to understand what causes life-threatening COVID-19 pneumonia in a minority of infected individuals. Recent work has suggested that preexisting autoimmunity against type I IFN can underlie critical COVID-19 pneumonia in greater than 10% of cases (). Here, we have confirmed that neutralizing autoantibodies to type I IFNs are present in patients with severe to critical COVID-19 from two independent cohorts, showing a combined prevalence of about 9%. Both the COVID-19+ and asymptomatic cohorts studied here had substantial Hispanic representation, a population that has not been previously studied at scale. The presence of anti–type I IFN autoantibodies in this population indicates that this phenomenon may be widely conserved across a diversity of ancestries. In terms of age and gender, the majority of autoantibody-positive severe to critical COVID-19 cases were male and over 55 years of age, consistent with previous reports (). We observed roughly equal numbers of autoantibody-positive males and females in the general population of San Francisco, as assessed through our community survey. Further study will be required to determine whether differences exist between the autoantibodies in male versus female individuals before SARS-CoV-2 infection that could partially explain the downstream skewing of hospitalized patients toward male gender, such as differential neutralization ability or additional accompanying risk factors.

In four patients with critical COVID-19 and tested positive for anti–type I IFN autoantibodies, longitudinal analysis determined that these autoantibodies were present from the earliest time point in their clinical course (collected within 4 to 13 days from the start of their first disease symptoms). Given the time required for a detectable, stable humoral immune response to form (2 to 3 weeks) (), our data suggest that autoantibodies predate infection with SARS-CoV-2. Consistent with this, our survey of a community cohort in the San Francisco Mission District also revealed a subset of presumed uninfected individuals who were anti–type I IFN antibody positive (0.3%), suggesting that there are individuals who may be at higher risk for critical COVID-19 due to these preexisting autoantibodies, including both males and females across a broad range of ages. Moreover, in a community-based population study, we did not detect these autoantibodies in 154 patients with asymptomatic or ambulatory infection with SARS-CoV-2 (compared to 13 of 3821 uninfected individuals) or in CCP donor samples, suggesting that the penetrance of severe to critical COVID-19 in infected individuals with autoantibodies is potentially complete.

In addition to validating the presence of anti–type I IFN autoantibodies in patients with severe and critical COVID-19, we further showed, using single-cell transcriptomic analysis, that these antibodies are associated with impaired ISG-I response in several distinct myeloid populations. Although other similar high-dimensional immune profiling studies have found evidence of impaired ISG responses in monocytes (, ) and neutrophils (), we have now provided additional specificity and a potential mechanism of how this may unfold in a subset of individuals. We also found impaired myeloid ISG-I expression in additional individuals with critical disease without detectable anti–type I IFN autoantibodies. This observation suggests that impaired type I IFN immunity is a shared mechanism of more severe forms of the disease in patients with and without autoantibodies to type I IFNs (). It is possible that patients without detectable autoantibodies may have had lower titers of autoantibodies, autoantibodies that neutralized lower amounts of type I IFNs, or autoantibodies that were undetectable because they were bound to type I IFNs. Alternatively, these patients may carry inborn errors of the production and amplification of type I IFNs, as recently shown in other patients () or autoantibodies against other cytokines or proteins in type I IFN response (). It is also possible that other antibody-mediated mechanisms may exist that are independent of the direct binding to IFNs (). Genetic and immunological studies are underway in our cohort of patients. These findings, along with the observation of high ISG-I expression in moderate patients early during the disease course that quickly diminishes, further suggest that impaired type I IFN immunity during the first hours and days of infection may account for the protracted disease course, including pulmonary and systemic inflammation. A two-step model of life-threatening COVID-19 is emerging, with defective type I IFN intrinsic immunity in the first days of infection resulting in viral spread, in turn unleashing leukocyte-mediated excessive inflammation in the lungs and other infected organs during the second week of infection ().

Our analysis of 189 cell surface proteins identified the expression of LAIR1 in cMs to be elevated in patients with COVID-19 and correlated with the impaired ISG-I response. LAIR1 is an inhibitory surface protein first described to be expressed in lymphocytes and involved in inhibiting NK-mediated cell lysis and effector T cell cytotoxicity upon FcR-mediated cross-linking (, , ). More recently, it has also been shown in monocytes and pDCs that LAIR1 cross-linking or binding to cognate ligands collagen and C1Q can inhibit the production of IFNα in response to Toll-like receptor ligands in healthy controls and patients with systemic lupus erythematosus (, ). LAIR1 expression is highest in cMs at the time of initial hospitalization and decreases rapidly by day 4 among a subset of critical patients, including the four with anti–type I IFN autoantibodies. In a recent study, a large array of autoantibodies were characterized in patients with COVID-19, among which autoantibodies against LAIR1 that were found to be highly specific to severe-to-critical COVID-19 (). Whether LAIR1 plays a causal role in deficient type I IFN response would require further investigation. Nevertheless, the ability to use a highly cell type–specific surface protein to predict impaired type I IFN response in patients with critical COVID-19 early during disease provides an important biomarker.

Our study has several limitations. First, we do not formally have age-, sex-, or ancestry-matched individuals presenting with different disease severity or in our control groups for the single-cell analysis. However, we did not detect associations to age or sex on cell composition or ISG-I expression for critical patients with or without autoantibodies. Second, there is a slightly longer delay between symptom onset and hospitalization with increasing disease severity (moderate: 2 to 11 days, 6.5 days average; severe: 2 to 16 days, 10.4 days average; critical: 6 to 36 days, 13.2 days average) and a bias toward specific patients with comorbidities or more severe disease at later time points. Nevertheless, these differences in the dynamics of the disease course do not fully explain the dynamics of ISG-I expression because the responses in moderate and severe patients never reach the low expression observed in critical patients even after 14 days of hospitalization, something that is in line with previous literature (). Third, because we only had sequencing data from four anti–IFN-α2 autoantibody-positive patients, there was limited power to compare critical cases with and without anti–IFN-α2 autoantibodies.

In conclusion, our findings have several important implications for the ongoing COVID-19 pandemic and our understanding of patients with a critical clinical course. First, our results show that an impaired ISG-I response early in the disease course in multiple immune populations is associated with autoantibodies to type I IFNs, providing a glimpse into the immune dysregulation present in patients with a severe clinical course. In this regard, it is critical to be able to identify patients with an impaired ISG-I response early during disease course; a combination of the highly specific assays for autoantibodies against type I IFNs and biomarkers for deficient ISGs, such as LAIR1, could quickly allow triaging of patients during initial hospitalization. Second, treatment strategies with IFN-β might be particularly valuable for those with preexisting antibodies to type I IFNs because these appear to neutralize only IFN-α or IFN-ω (). The large immunological differences of severe to critical patients in the earliest time points additionally suggest that identification and treatment would likely need to happen early in the disease course. Third, we found that autoantibodies to type I IFNs in critical COVID-19 cases were present at the time of their presentation and precede the development of antibodies to SARS-CoV-2. This, along with the presence of healthy autoantibody-positive individuals in the community (), suggests that anti–type I IFN autoantibodies predate infection and that there exists an at-risk group for severe to critical disease in the general population. Going forward, strategic efforts to identify this high-risk population early in the disease course could have substantial impact on improving clinical outcomes including mortality rates, and identifying these individuals before infection could have a major impact on preventive measures.

MATERIALS AND METHODS

Study design

This study was designed to identify immunological features in the blood from hospitalized patients presenting with COVID-19 symptoms associated with disease status and severity. For this purpose, PBMCs and serum samples were collected from participants enrolled in the COMET cohort, which partially includes patients recruited to the Immunophenotyping Assessment in a COVID-19 Cohort (IMPACC) study (data file S4). The materials collected from each of the participants in this study were coded by donor ID numbers such that the experimentalist could not define the disease status of each participant. Hence, in each experiment, we randomized the processing of PBMCs from 22 to 25 participants, including 16 to 23 patients and 2 to 8 healthy individuals. All patients hospitalized for symptomatic COVID-19 infection at both the tertiary care center and the safety-net county hospital associated with the University of California San Francisco (UCSF) were eligible to participate in the COMET cohort study. Biospecimens may have been collected under an Institutional Review Board (IRB)–approved initial waiver of consent with subsequent attempts to consent surrogates and study participants for full study participation. We selected samples from 69 of 101 of participants enrolled in COMET between 8 April 2020 and 20 June 2020 (Table 2 and data file S1), which resulted in at least 15 to 20 donors per COVID-19 severity group. According to previous publications of similar cohort sizes, our study should have sufficient power to detect COVID-19–associated cell type compositional and gene expression changes of similar effect sizes as described (, ). The healthy anonymized donors were chosen within the age range of patients with COVID-19. Sample selection was prioritized in the patients that were hospitalized with longitudinal blood collections, and therefore, PBMCs were available across a 14-day period. This study is approved by the IRB of the UCSF Human Research Protection Program (IRB no. 20-30497).

A total of 231 patients with COVID-19 and 96 healthy controls were enrolled in the ZSFG cohort (Table 3 and data file S2). All ZSFG patients were collected between 1 March 2020 and 21 July 2020 and had positive results by SARS-CoV-2 reverse transcription PCR (RT-PCR) in nasopharyngeal swabs. Clinical data were extracted from electronic health records and included demographic information, major comorbidities, patient-reported symptom onset date, symptoms, and indicators of disease severity. Healthy, pre–COVID-19 control plasma samples were obtained as deidentified samples from the New York Blood Center. These samples were part of retention tubes collected at the time of blood donations from volunteer donors who provided informed consent for their samples to be used for research.

Details of the community-based cohort including 4041 participants were previously described by Chamie et al. () (Table 4 and data file S3). The four APS-1 patients in the study were collected at the NIH under protocol #11-I-0187 and were previously published by Ferre et al. (, ). CCPs were collected in the Vitalant system following U.S. Food and Drug Administration (FDA) guidance for donor eligibility. These criteria evolved throughout the study period because of testing availability and evolution of the pandemic in the United States. Evidence of COVID-19 was required in the form of a documented positive SARS-CoV-2 molecular or serologic test and complete resolution of symptoms initially at least 14 days before donation, but then, a minimum of 28 days was implemented. All CCP donors were also required to meet traditional allogeneic blood donor criteria. At the time of plasma collection, donors consented to use of deidentified donor information and test results for research purposes. All CCPs were tested for SARS-CoV-2 total Ig antibody using the Ortho VITROS CoV2T assay at our central testing laboratory (Creative Testing Solutions). CCP qualification requires the signal-to-cutoff (S:CO) ratio of this test to be at least 1.0. Retention samples of serum and plasma for all donations are archived at the Vitalant Research Institute Denver. Plasma samples are from 175 unique CCP donors, where some had repeated donations, for a total of 281 samples. These samples were selected solely on the Ortho VITROS CoV2T assay results to represent the entire range of high to low S:CO ratio. Collections were from across the U.S. Vitalant system from 8 April 2020 to 1 September 2020.

Isolation and preparation of PBMCs for scRNA-seq

Whole blood from 80 donors was drawn into plastic EDTA Vacutainer blood collection tubes (Becton, Dickinson and Company) at the time of hospital admission (D0) and 4 (D4), 7 (D7), and 14 (D14) days later. Of these donors, 69 were patients with high clinical suspicion of COVID-19 infection that were admitted at UCSF or Zuckerberg San Francisco General Hospital and Trauma Center (ZSFG) and 11 were healthy donors. COVID-19 status was assessed for all 69 patients by RT-PCR tests of nasal swab samples and confirmed that 15 patients were COVID-19 negative, whereas 54 patients were COVID-19 positive. During the hospitalization, the severity of each COVID-19–positive patient was assessed using the NIH COVID-19 severity scale (Table 1) (). For all analyses, we categorized patients on the basis of their severity at time of hospital admission (D0).

PBMCs were isolated at room temperature using SepMate PBMC Isolation Tubes (STEMCELL Technologies) by the UCSF Biospecimen Resource Program. Briefly, 7.5 ml of whole blood was centrifuged at 1000 relative centrifugal force (rcf) for 10 min with swinging-out rotor and brake off to separate 3.5 ml of plasma. The remaining blood was diluted with 8 ml of Dulbecco’s phosphate-buffered saline (DPBS) and slowly added to a SepMate-50 tube prefilled with 15 ml of Lymphoprep (STEMCELL Technologies). The tube was then centrifuged at 800 rcf for 20 min with brakes off. After centrifugation, the top layer including the buffy coat was gently and quickly poured into a 50-ml Falcon tube to centrifuge at 400 rcf for 10 min with brake on. The pellet was washed twice each time with 20 ml of EasySep buffer (STEMCELL Technologies) followed by centrifugation at 400 rcf for 10 min with brakes on. Washed PBMCs were counted and resuspended at 5 × 106 cells/ml in cold freezing media (fetal bovine serum with 10% dimethyl sulfoxide solution). Cells were aliquoted into cryovials at 1 × 106 to 5 × 106 cells per vial and transferred in a Mr. Frosty (Nalgene) to the −80°C freezer for 24 hours before cryopreservation in liquid nitrogen.

Single-cell multimodal immunophenotyping

Multiplexed single-cell sequencing was performed following a previously published protocol () and manufacturer’s user guide (Document CG000186 Rev D, 10X Genomics). The complete protocol is available on protocols.io (www.protocols.io/view/10x-citeseq-protocol-covid-19-patient-samples-tetr-bqnqmvdw). In each experiment, PBMCs from 22 to 25 participants were used, including 16 to 23 patients and 2 to 8 healthy individuals. Longitudinal samples of the same patient were used in different experiments to allow genetic demultiplexing. Each experiment used samples collected at different longitudinal time points to prevent experimental conditions from coinciding with potential batch effects.

Briefly, PBMCs were thawed in a 37°C water bath for 30 s and washed with 5 ml of complete RPMI medium followed by centrifugation at 350 rcf for 5 min at room temperature. Cell counts and viability were determined using a Vi-CELL XR automated cell counter (Beckman Coulter Life Sciences). An equal number of cells were aliquoted from each sample to create a pool of 1.5 × 106 cells with an average viability of 85% or higher. Pooled PBMCs were blocked with Human TruStain FcX (BioLegend) in cell staining buffer (BioLegend) for 10 min on ice, followed by staining with a customized TotalSeq-C human cocktail for 45 min on ice (data file S5). Cells were washed three times, resuspended in PBS with 0.04% bovine serum albumin, filtered through a 40-μm filter, and counted with Countess II Automated Cell Counter (Thermo Fisher Scientific). Single-cell suspensions were loaded into a Chromium Single Cell Chip A for single cell encapsulation using the 10X Chromium Controller according to the manufacturer’s user guide (Document CG000186 Rev D, 10X Genomics) and as previously described (). In each experiment, the pooled cells of 22 to 25 participants were loaded into four to six individual lanes aiming for 7 × 104 loaded cells per lane.

Single-cell library preparation and sequencing

Single-cell libraries were constructed following the manufacturer’s user guide (Document CG000186 Rev D). Complementary DNA (cDNA) libraries were generated using the Chromium Single Cell 5′ Library & Gel Bead Kit and i7 Multiplex Kit. Surface protein Feature Barcode libraries were generated with Chromium Single Cell 5′ Feature Barcode Library Kit and i7 Multiplex Kit N, Set A. In total, libraries generated from 971,550 cells were PE150-sequenced at the Chan Zuckerberg Biohub on 18 lanes of an Illumina NovaSeq 6000 sequencer using a NovaSeq 6000 S4 Reagent Kit v1.

Genotyping, sample demultiplexing, and doublet detection

To assign cells to donors of origin in our multiplexed design, we used the genetic demultiplexing tools Freemuxlet and vcf-match-sample-ids, each a part of the Popscle suite of statistical genetics tools (https://github.com/statgen/popscle). Freemuxlet leverages the single-nucleotide polymorphisms (SNPs) present in transcripts and performs unsupervised clustering on the droplet barcodes to assign each to a nameless donor or assign them as doublets between genetically distinct nameless donors. The algorithm takes in a list of candidate loci throughout the genome at which to scan for SNPs and returns droplet barcodes with donor assignments and a set of observed variants per donor. These sets of variants are then matched using genotypic similarity to those from an orthogonal bulk RNA-seq assay, done on an individual basis, to determine which donor is which patient. Once nameless donors are matched to uniquely identifiable patients, droplet data are then joined with the other clinical covariates available for the patients, including age, sex, race, and disease status. Freemuxlet was run on each of the nine droplet reaction runs separately, using a list of exonic SNPs that were expected to be found in the 5′-end scRNA-seq data and that have a minor allele frequency >0.05, based on data from the 1000 genomes project ().

Bulk RNA sequencing

Bulk RNA-seq data were generated to extract genotype information so that single cells could be demultiplexed and matched to single donors. For each donor, RNA was extracted from PBMCs using the Quick-RNA MagBead Kit (Zymo Research) on a KingFisher Flex system (Thermo Fisher Scientific) according to the manufacturer’s protocol. RNA integrity was measured with the Fragment Analyzer (Agilent), and library generation was continued when integrity was at least 6. Total RNA-seq libraries were depleted from ribosomal and hemoglobin RNAs and generated using FastSelect (Qiagen) and Universal Plus mRNA-Seq with NuQuant (Tecan) reagents. Pooled libraries were PE100-sequenced on a HiSeq 4000 or PE150-sequenced on an Illumina NovaSeq 6000 S4 flow cell at the Chan Zuckerberg Biohub.

Single-cell epitope and RNA-seq preprocessing and alignment

CellRanger v3.1.0 (run 1 to 7, cDNA library generation of run 6 failed) or v4.0.0 (run 8 to 10) software with the default settings was used to demultiplex the sequencing data and generate FASTQ files (Cellranger mkfastq), align the sequencing reads to the hg38 reference genome, and generate a unique molecular identifier (UMI)–filtered gene and protein expression count matrix for each lane (Cellranger count for scRNA-seq and CITE-seq data). Count matrices were then concatenated across all 50 lanes to generate two matrices: one mRNA matrix with 971,550 cells and 36,601 genes, and one surface protein matrix with 971,550 cells and 189 proteins.

Single-cell epitope and RNA-seq processing and quality control

Resulting gene and protein expression count matrices were further processed in the Python package Scanpy v1.5.1 (). Processing of the concatenated mRNA count matrix was done using a two-step process. We have found empirically that traditional workflows for the processing of droplet-based RNA-seq data for PBMCs can sometimes create unwanted effects in the downstream end points. In particular, filtering of cells with a high percentage of mitochondrial cells may create visual artifacts in UMAP projections, and filtering of the matrix to only a few hundred highly variable genes and reducing the memory footprint of the data can sometimes lead to spurious clustering of cells based on only a few genes. Our iterative process yields a UMAP projection that captures all available heterogeneity with minimal filtering in the first iteration, then removes nontarget cells, and corrects for nonbiological signal in the second iteration. By doing this, we use a relatively large number of components to inform projections and clustering but observe that the outputs in our dataset match the known biology better (such as proximity of similar cell types and states in UMAP space) and yield higher-confidence annotations.

In the first step, the mRNA matrix was filtered to remove doublet droplets, as annotated by Freemuxlet, and very lowly expressed genes with less than 100 UMIs across the nine runs. The matrix was then normalized to yield a constant UMI sum per cell and log-transformed. Matrix values were scaled to yield a mean of zero and SD of 1 per gene. Principal components analysis (PCA) was performed, followed by nearest neighbors, UMAP projection, and Leiden clustering, using an input of the 150 principal components with the highest variance explained and otherwise default Scanpy parameters. At this stage, Leiden clustering resolution was adjusted and restricted to certain clusters to separate out clusters of cells that projected into similar UMAP space. These clusters were subsequently annotated to mark those with a high percentage of mitochondrial content, which typically represent dead or dying cells, as well as mark clusters with high expression of hemoglobin and platelet factor 4, representing the nontarget cell types of red blood cells and platelets. At this stage, we also observed that there were prominent batch effects in UMAP space that required correction.

In the second iteration, nontarget cell types marked in the first step were removed before processing. Then, the same processing was followed starting from the raw data, with the exception that ComBat batch correction () was performed (to correct for the “run” covariate) after scaling and before performing PCA. Last, further filtration of a relatively small number of cells with high expression of platelet, red blood cell, and mitochondrial genes was performed, as well as removal of donors that declined study participation. After processing, 600,929 cells and 18,262 genes remained in the mRNA matrix.

The surface protein matrix was filtered to the cells found in the mRNA matrix. One protein was removed from the matrix, as it appeared to have very low counts relative to the other surface proteins. The remaining proteins were then normalized using the centered-log ratio (CLR), computed for each gene independently. The CLR has typically been used for CITE-seq data with the recognition that antibody counts are typically not zero-inflated and flow cytometry–like Gaussian distributions are achievable when treating the data as compositional. However, we recognized that the CLR-normalized distributions were affected by a relatively small number of cells that had extremely high or low expression, skewing the visualization of the Gaussian mixture distributions. To remedy this effect, we identified boundary values of the distributions for each gene using a bin height threshold when values were plotted on a histogram, clipped the values to these boundaries, and scaled the remaining values between 0 and 1.

Cell type classification

After processing, Leiden clustering was adjusted to match the clustering of cells projected into UMAP space. Cell type annotation was performed at three levels of granularity based on known marker gene and protein expression, as well as differentially expressed genes between clusters using a Wilcoxon rank-sum test. At the lowest level of granularity, we identified 11 cell types corresponding to the known major cell types present in PBMCs: CD4T, CD8T, γδT, cMs, ncMs, NK, B, PBs, cDC, pDC, and Progen cells. At the next level of granularity, we separate out memory, naïve, and proliferating subtypes in the lymphocytes; two different subtypes known in each of the NK cells and conventional DCs; regulatory T cells from the CD4T group; mucosal-associated invariant T cells from the CD8T group; subpopulations of lineage-committed progenitor cells; and a subpopulation of B cells that seemed to be committed to the PB lineage. At the highest level of granularity, we further separate out a CD8+ effector memory T cell population, an NK population with CD3 transcript expression, early and late proliferating subpopulations in the lymphocytes, and a few subpopulations that were either donor specific (patient 1002 in the B cells and patient 1006 in the monocytes) or run specific [such as cells from run 3, which exhibited some processing issues and for which ComBat () was unable to correct].

Differential proportion and expression analysis

Differences in cell type proportions were assessed in COVID-19+ versus COVID-19− cases or healthy controls by aggregating cell type observations per COVID-19 status at time point D0. In addition, in the COVID-19+ cases for which all four time points were available, cell type proportion changes were assessed over time. Differences in gene expression were determined for each of the myeloid cell types between COVID-19+ cases at D0, D4, D7, or D14 versus healthy controls. To assess whether these changes are specific for COVID-19+ cases or are a more general phenomenon due to acute respiratory distress syndrome, we also compared the up-regulated genes with the COVID-19− cases. Differential expression analysis was performed per run on the raw gene expression matrix using Memento v0.0.4 (doi: 10.5281/zenodo.5172943), after which results were meta-analyzed over all runs. Genes were prefiltered on the basis of a minimum raw mean expression of 0.07 within at least 90% of both comparison group.

Differential expression heatmaps

Heatmaps show the pseudobulked, z-scored expression values of the donors present at each time point for the top significantly up-regulated genes. To generate the heatmaps, cells were first subsetted from the larger mRNA matrix to only those of a given cell type. Counts were pseudobulked across all genes by patient present at each time point, yielding a single gene-by-sample matrix, with 179 unique donor–time point samples. The genes in this matrix were subsetted to the union of the top 150 genes with the highest differential expression coefficient at each time point, using the one-dimensional memento results that tested gene counts in COVID-19+ cases versus healthy controls. Genes were further filtered to remove those that had high variance in healthy controls (SD > 0.5) because these were enriched for what seemed to be a nonbiological signal (such as ribosome-associated genes). The matrix, now with 204 genes, was then z-scored and separated by time point to four matrices, with healthy samples being distributed to each matrix.

Ordering of the rows and columns were computed such that they would be consistent among the heatmaps. Genes were clustered by k means using only the values of the day 0 time point, with k = 6 chosen by the “elbow” point of the graph plotting distortion (using a sum of square errors cost function) with increasing numbers of clusters. Columns of each heatmap were determined by taking the 80 columns across all heatmaps that had the earliest time point for each donor, subsetting according to their disease status (healthy, COVID-19−, and COVID-19+), and then hierarchically clustering within each of those groups. With this ordering, each donor then had a unique position along the horizontal axis, which was then applied to all the heatmaps, omitting those samples that were absent from a given time point. GSEA was done using the GOATOOLS Python package (), filtering to terms with at least two associated genes.

ISG score method

An orthogonal scRNA-seq dataset containing PBMCs stimulated with IFN-β and IFN-γ was used to identify the specific and shared type I and type II ISGs in the cMs (Gene Expression Omnibus accession number GSE181897). The gene list was used to calculate a type I, type II, and shared ISG score based on the average gene expression count of the unique or shared type I and type II ISGs. These ISG scores were calculated for each unique combination of cell type, donor, and time point. Subsequently, ISG scores were averaged over each of the disease categories (COVID-19+ moderate-severe, COVID-19− moderate-severe, healthy control) and then log2-transformed.

Flow cytometry validation of elevated surface LAIR1 in cMs for patients with COVID-19

Frozen PBMCs from patients with available samples (8 critical, 9 moderate, 8 severe, 10 negative controls, and 5 healthy controls) were thawed as previously mentioned. Each sample (5 × 105 cells per donor) was then resuspended to 75 μl of cell staining buffer and aliquoted to a 96-well U-bottom plate (Genesee Scientific catalog no. 25-221) before blocking with 5 μl of TruStain FcX for 10 min on ice. Cells were then stained with a cocktail of CD14–Brilliant Violet 421 (BioLegend, clone M5E2, catalog no. 301830), CD16–Alexa Fluor (AF) 488 (BioLegend, clone 3G8, catalog no. 302019), SIGLEC1-phycoerythrin (BioLegend, clone 7-239, catalog no. 346004), and LAIR1-AF647 (BioLegend, clone NKTA255, catalog no. 342802) antibodies in a total volume of 100 μl (5 μl per antibody) for 30 min on ice in the dark. Cells were washed twice in 200 μl of cell staining buffer before resuspension in 300 μl of PBS. Cells were stained with propidium iodide (PI; final 1 μg/ml) before flow analysis on an LSR II (UCSF Parnassus Flow Cytometry Core). Analyzed cells were gated for singlets, live (PI−), and CD14+CD16− before mean fluorescence intensity (MFI) calculations.

SARS-CoV-2 detection by clinical RT quantitative PCR

Viral titers were quantified in a subset of the COVID-19+ patients in the UCSF CLIAHUB Clinical Microbiology Laboratory. For this, RNA was extracted from tracheal aspirate and nasopharyngeal swabs and used for RT-qPCR (quantitative PCR) as previously described (). In short, viral titers were assessed using primers targeting the SARS-CoV-2 N gene (Ct1) and E gene (Ct2) and human ribonuclease (RNase) P gene (Ct_host, positive control). The Ct value of the viral N or E gene was subtracted from the human RNase P gene (ΔCt: ΔCt1 and ΔCt2) (data file S9), and number signs were reversed to obtain a measurement for normalized viral RNA abundance. As there was an almost perfect correlation between ΔCt1 and ΔCt2 values (Pearson R = 0.99, P = 1.9 x 10−283) and ΔCt2 had the least missing values, viral RNA abundance is represented by the ΔCt2 values. ΔCt2 values as measured in the tracheal aspirate and nasopharyngeal swab samples were linked to the scRNA-seq PBMC data of the same donor and the closest possible time point (up to 2 days apart).

RLBA for anti–IFN-α2 autoantibody detection

A DNA plasmid containing full-length cDNA sequence with a Flag-Myc tag (OriGene, #RC221091) was verified by Sanger sequencing and used as template in T7-promoter–based in vitro transcription/translation reactions (Promega, #L1170) using [S35]-methionine (PerkinElmer, #NEG709A). IFN-α2 protein was column-purified using NAP-5 columns (GE Healthcare, #17-0853-01); incubated with 2.5 μl of serum, 2.5 μl of plasma, or 1 μl of anti-myc–positive control antibody (Cell Signaling Technology, #2272); and immunoprecipitated with Sephadex protein A/G beads (4:1 ratio; Sigma-Aldrich, #GE17-5280-02 and #GE17-0618-05) in 96-well polyvinylidene difluoride filtration plates (Corning, #EK-680860). The radioactive counts [counts per minute (cpm)] of immunoprecipitated protein were quantified using a 96-well MicroBeta TriLux liquid scintillation plate reader (PerkinElmer). Antibody index for each sample was calculated as follows: (sample cpm value – mean blank cpm value)/(positive control antibody cpm value – mean blank cpm value). For the COVID-19 patient and CCP cohorts, a positive signal was defined as greater than 6 SDs above the mean of pre–COVID-19 blood bank noninflammatory controls. For the large asymptomatic San Francisco community population cohort, a positive signal was defined as having a z score greater than 3.3 (P = 0.0005) relative to the whole cohort.

Luciferase reporter assays

The blocking activity of anti–IFN-α and anti–IFN-ω autoantibodies was determined by assessing a reporter luciferase activity. Briefly, human embryonic kidney (HEK) 293T cells were transfected with the firefly luciferase plasmids under the control human ISRE promoters in the pGL4.45 backbone and a constitutively expressing Renilla luciferase plasmid for normalization (pRL-SV40). Cells were transfected in the presence of the X-tremeGENE 9 transfection reagent (Sigma-Aldrich) for 36 hours. Dulbecco’s modified Eagle’s medium (Thermo Fisher Scientific) was supplemented with 10% healthy control or patient serum or plasma and was either stimulated with IFN-α or IFN-ω (10 ng/ml) or left unstimulated for 16 hours at 37°C. Each sample was tested once. Last, luciferase expression was measured with the Dual-Glo reagent according to the manufacturer’s protocol (Promega). Firefly luciferase values were normalized against Renilla luciferase values, and fold induction was calculated relative to controls transfected with empty plasmids.

Statistical analysis

Differential proportion analysis was performed using a permutation-based approach that compares observed cell type proportion differences with those calculated from a null distribution that is generated by randomly shuffling cell type labels (100,000 permutations) for a fraction (w = 0.1) of the total cells, as described previously (). Resulting P values were corrected for multiple testing using Holm’s correction, after which an adjusted P value of <0.05 was considered significant. For the differential expression analysis, an FDR of <0.05 was used to determine statistical significance. A Welch’s t test was performed to compare the ISG score between COVID-19+ patients and healthy controls. Significance was defined as Bonferroni-adjusted P value < 0.05. To detect changes in ISG-I score or LAIR1 protein expression over time, samples were first subsetted to COVID-19+ and COVID-19− patients for whom we had more than one time point (54 of 69). We then used an LMM for each group using the statsmodels Python package, with formula “ ~ time_point,” where was either type I score or LAIR1 expression. Slope deviations from 0 were considered significant when their P value was below 0.05. All correlations were calculated using the Pearson R. Significance was then defined as Holm’s corrected P value < 0.05.