Human SARS-CoV-2 challenge resolves local and systemic response dynamics

Study participants and design

16 healthy adults aged 18-30 years, with no evidence of a previous SARS-CoV-2 infections or vaccinations (seronegative), were included for this study from the wider cohort (34 participant) enrolled as part of the Human COVID-19 Challenge study, pioneered by the government task force, Imperial college london, Royal Free London NHS Foundation Trust and hVIVO⁶. These participants were enrolled as part of the cohort 5 and 6, from June - August 2021. Reported patient identifiers have been de-identified and are not known by the participants in this study. Volunteers were tested for the presence of anti-SARS-CoV-2 protein antibodies via the MosaiQ COVID-19 antibody microarray (Quotient) prior to enrollment and excluded based upon a positive test, as well as, upon risk factors assessed by clinical history, physical examinations and screening assessments. See Killingley, Mann et al, (2022)⁶ for the full list of inclusion and exclusion criteria and for further details regarding the challenge set up and ethics. In short, written informed consent was obtained from all volunteers before screening and study enrollment. The clinical study was registered with ClinicalTrials.gov (identifier NCT04865237). This study was conducted in accordance with the protocol, the Consensus ethical principles derived from international guidelines including the Declaration of Helsinki and Council for International Organizations of Medical Sciences (CIOMS) International Ethical Guidelines, applicable ICH Good Clinical Practice guidelines, applicable laws and regulations. The screening protocol and main study were approved by the UK Health Research Authority – Ad Hoc Specialist Ethics Committee (reference: 20/UK/2001 and 20/UK/0002).

Participant 674007, who fulfilled enrollment criteria, was later found to have low levels of neutralizing and spike-binding antibodies on admission to the quarantine upon more sensitive post-study experiments and analysis. This patient was classified as an abortive infection based on the virus kinetics (see virology method below), with the exclusion of this individual found not to alter any of our conclusions.

The participants were followed for 1 year after inoculation, with continued samples and metadata collected for the use in future/further studies and to benefit the research community. This study however focused primarily on the first 28-days post-inoculation (with the exception of 46 days for one participant as noted below, see Sample collection).

Challenge virus

Participants were inoculated intranasally with an wild-type pre-alpha SARS-CoV-2 challenge virus (full formal name: SARS-CoV-2/human/GBR/484861/2020) at dose 10 TCID50 at day 0. 100 µl per naris was pipetted between both nostrils and the participant was asked to remain supine (face and torso facing up) for 10 minutes, followed by 20 minutes in a sitting position wearing a nose clip after inoculation to ensure maximum contact time with the nasal and pharyngeal mucosa. Mid-turbinate nose and throat samples were collected twice daily using flocked swabs and placed in 3 ml of viral transport medium (BSV-VTM-001, Bio-Serv) that was aliquoted and stored at −80 °C in order to evaluate viral kinetics (infection status) as described in Virology method section below. Participants remained in quarantine for a minimum of 14 days post-inoculation until the following discharge criteria were met: two consecutive daily nose and or throat swabs with no viral detection or a qPCR Ct value >33.5 and no viable virus by overnight incubation viral culture with detection by immunofluorescence. For protocol/full details and ethics used within the Human SARS-CoV-2 challenge study see the Challenge virus section of methods in Killingley and Mann., et al (2022)⁶.

Sample collection

Clinical assessments

Participants were carefully monitored and assessed daily using an array of blood tests, spirometry, electrocardiograms and clinical assessments (vital signs, symptom diaries and clinical examination). Full details of all the safety and clinical data collected with the human SARS-CoV-2 challenge study can be seen in the methods in Killingley and Mann., et al (2022)⁶, with an overview of select metadata for the 16 participants enrolled in this study in Extended Data Table 1g.

Virology

Longitudinal measures from the nose and throat (pharyngeal) were carried out daily in order to assess and quantify the viral kinetics of each participant pre- and post-inoculation. These were measured using two independent assays: (1) qRT–PCR with N gene primers/probes adapted from the Centers for Disease Control and Prevention (CDC) protocol²⁷ (updated 29 May 2020) and (2) quantitative culture by Focus Forming Assay (FFA). For full details of each assay and statistical analysis refer to the methods in Killingley and Mann., et al (2022)⁶.

The lower limit of quantification (LLOQ) for RT-qPCR was 3 log10 copies per milliliter, with positive detections less than the LLOQ assigned a value of 1.5 log10 copies per milliliter and undetectable samples assigned a value of 0 log10 copies copies per milliliter. Only samples where participants presented with a positive RT-qPCR were further tested using the FFA assay. In the FFA the LLOQ was 1.27 FFU ml−1; viral detection less than the LLOQ was assigned 1 log10 FFU ml−1; and undetectable samples were assigned 0 log10 FFU ml−.

A sustained laboratory-confirmed infection was defined as quantifiable RT–qPCR detection greater than the LLOQ from mid-turbinate and/or throat (pharyngeal) swabs on two or more consecutive 12-hourly time points, starting from 24 hours after inoculation and up to discharge from quarantine. Participants where only stand alone RT–qPCR tests returned quantifiable results (> LLOQ) were classified as transient infections. Participants where no RT–qPCR tests returned quantifiable results (> LLOQ) were classified as abortive infections (See Extended Data Fig. 1b and Extended Data Table 1a,b,h,i).

Infection intervals for each participant were calculated based on the time of the first and last RT-qPCR test with detectable virus (across the nose and/or throat), where timepoints in which tests below the LLOQ (1.5) were also counted if they occurred < 2 days of a quantifiable (> LLOQ) test result.

Nasopharyngeal swab dissociation and processing

Following freezing, nasopharyngeal swabs were transferred to a category level 3 facility at University College London were stored and processed in batches of 7-to -8 samples at a time to a single cell suspension. All work was carried out in a MSC class I hood in compliance with standard category level 3 safety practices. The dissociation and collection of cells from nasopharyngeal swab was carried out in accordance with the previously described protocol^{28, 29}, with minor modifications. This approach involves multiple parallel washes and digestion steps using both the nasopharyngeal swab and collected freezing/ wash media to help ensure the maximum cells and cellular material is collected. First samples are exposed to DTT for 15 mins, followed by an accutase digestion step for 30 mins before cells from the same sample (collected directly from the swab or the freezing media/washes from that swab) are quenched, pooled and filtered prior to checking the cell number and viability.

Briefly, samples were rapidly thawed and the liquid collected in an empty 15mL falcon tube (Tube B). The cryovial, lid and swab was then carefully rinsed with 3x 1mL warm RPMI 1640 medium which was added dropwise to the 15 mL tube whilst gently swirling the tube, in order to slowly dilute the DMSO from the freezing media to help present the cells bursting. After waiting 1 min, the tube (Tube B) was then topped up with an extra 2 mL of warm RPMI 1640 media and centrifuged at 400 g for 5 min at 4°C. The cell pellet was then resuspended in RPMI 1640/10 mM DTT (Thermofisher, R0861), and incubated for 15 min on a thermomixer (37°C, 700 rpm), centrifuged as above and the supernatant was aspirated and the cell pellet was resuspended in 1 mL Accutase (Merck, A6964-500ML). This was then incubated for a further 30 min on the thermomixer (37°C, 700 rpm).

In parallel to the processing of the cell freezing media/washes above, the swab was moved to an new 1.5 mL eppendorf tube (Tube C) containing 1 mL RPMI 1640/10 mM DTT and placed on the thermomixer (37°C, 700 rpm) for 15 minutes. In accordance with the steps above, the swab was next transferred to a new 1.5 mL eppendorf (Tube D) containing 1 mL Accutase and incubated with agitation (700 rpm) at 37°C. The 1 mL RPMI 1640/10 mM DTT from the nasopharyngeal swab incubation (in Tube C) was centrifuged at 400 g for 5 min at 4°C to pellet cells, the supernatant was discarded, and the cell pellet was resuspended in 1 mL Accutase and incubated for 30 min at 37°C with agitation (700 rpm).

Following the Accutase digestion step, all cells were combined (Tube B, C and D) and filtered using a 70 µm nylon strainer (pre-wetted with 3 mL quenching media: RPMI 1640/10% FBS/1 mM EDTA (Invitrogen, 1555785-038) in a 50 mL conical tube (Tube E).

The filter, tubes and swab were then further thoroughly rinsed with quenching media in order to collect all cells and the washes combined. The dissociated, filtered cells (Tube E) were then centrifuged at 400 g for 5 min at 4°C, and supernatant discarded. The cell pellet was resuspended in residual volume (∼500 µL) and transferred to a new 1.5 mL eppendorf tube (Tube F). Tube E was then washed with a further 500 µL of RPMI 1640/10% FBS and combined with Tube F, centrifuged as above, supernatant removed and cells resuspended in 20 µL of RPMI 1640/10% FBS. Using Trypan Blue, total cell counts and viability were assessed. The cell concentration was adjusted for 7,000 targeted cell recovery according to the 10x Chromium manual before loading onto the 10x chip (between 700–1,000 cells per µl) and processing immediately for 10x 5′ single-cell capture using the Chromium Next GEM Single Cell V(D)J Reagent Kit v1.1 (Rev E Guide). For samples where fewer than 13,200 total cells were recovered, all cells were loaded.

Note: Due to the sample type, necessary freezing process and no access to a class 3 flow facility to sort out viable cells, the majority of the samples processed were seen to have low viability (ranging from 5.4% viability to 57.85, with the average viability of samples processed of 26.89%).

PBMC CITE-seq staining for single-cell proteogenomics

Frozen PBMC samples were thawed and processed in batches of 16 to allow for a carefully designed pooling strategy. Here each sample was pooled twice into two unique pools containing up to four PBMC samples per pool from mixed timepoints. Note: Only one sample from each donor was even pooled together at a time to assist with the demultiplexing later. This pooling strategy was used to help remove and correct for any protocol-based batch effects.

In brief, PBMC samples were thawed quickly at 37 °C in a water bath. Warm RPMI 1640 medium (20– 30 mL) containing 10% FBS (RPMI 1640/FBS) was added slowly to the cells before centrifuging at 300g for 5 min. This was followed by a wash in 5 mL RPMI 1640/FBS. The PBMC pellet was collected, and the cell number and viability were determined using Trypan Blue.

PBMCs from 4 different donors were then pooled together (1.25 × 105 PBMCs from each donor) to make up 5.0 × 105 cells in total. The remaining cells were used for DNA extraction (Qiagen, 69504). The pooled PBMCs were resuspended in 22.5 µl of cell staining buffer (BioLegend, 420201) and blocked by incubation for 10 min on ice with 2.5 µl Human TruStain FcX block (BioLegend, 422301). The PBMC pool was then stained with TotalSeq-C Human Cocktail, V1.0 antibodies (BioLegend, 399905) according to the manufacturer’s instructions (1 vial per pool). For a full list of TotalSeq-C antibodies (130 abs + 7 isotype controls) refer to Extended Data Table 1j. Following a 30 min incubation period with the TotalSeq-C Human Cocktail V1.0 antibodies (at 4 °C in the dark) the PBMCs were topped up using cell staining buffer and centrifuged down to a pellet (500g for 5 min at 4 °C) and discarding the supernatant. The pellet was then resuspended and washed in the same manner two more times using the resuspension buffer (0.05% BSA in HBSS), before finally being resuspended in a 20-30 µL resuspension buffer and counted again. The PBMC pools were then processed immediately for 10x 5′ single cell capture (Chromium Next GEM Single Cell V(D)J Reagent Kit v1.1 with Feature Barcoding technology for cell Surface Protein-Rev D protocol). 25,000 cells were loaded from each pool onto a 10x chip.

PBMC Dextramer staining for SARS-CoV-2 antigen specific T cell enrichment and single-cell sequencing

In order to further validate and investigate the SARS-CoV-2 antigen-specific T-cell populations in our single cell dataset Day 10, 14 and 28 post-inoculation PBMCs samples, from all 16 participants, were further enriched and processed for single cell sequencing using a multi-allele panel of 44 SARS-CoV-2 antigen specific dCODE™ Dextramer® (10x compatible) (Immudex, see Extended Data Table 1k for full panel). This panel includes; five antigen-specific T-cell populations, spanning four MHCI and one MHCII alleles (covering a total of 15 participants, see Extended Data Table 1l) and several negative controls. Samples were then stained with several FACS antibodies (for monocyte and T cells) and sorted using MACSQuant® Tyto® Cell Sorter (Miltenyi Biotec), where PE-dCODE Dextramer®-positive cells were collected and processed for 10x 5′ single cell capture. This allowed for the quantification of paired clonal TCR sequence and TCR specificity by overlaying single-cell V(D)J expression onto dCODE Dextramer®-positive cell clusters.

The dextramer staining protocol was taken from Immundex and optimised/adapted to suit our samples and pooling/staining strategy. In brief, the PBMC samples were thawed in batches of 7 to 8 samples and the cell number and viability for each sample calculated using Trypan Blue as previously described above. All cells from each sample were then pooled together in a fresh 1.5 mL eppendorf tube. Note: The pooling strategy here was such that only one sample per participant/donor was used per pool in order to enable de-multiplexing by genotype later on and each pool containing a mixture of timepoints to help reduce batch effect. In order to ensure the collection of as many cells as possible, each of the original sample tubes was then washed with 200 µl of staining buffer (1X PBS pH 7.4 containing 5% Heat inactivated FBS (Thermo Fisher Scientific, 10500064) and 0.1g/l Herring sperm DNA (Thermo Fisher Scientific, 15634017) and added to the pool. The tube was then topped up to 1.4 mLs with staining buffer and centrifuged down to a pellet (400g for 5 min at 4 °C). The supernatant was carefully removed and the cell pellet gently resuspended in a total of 30-40 µl staining buffer depending on pellet’s size, ready for staining.

In parallel the dCODE™ Dextramer® master mix was prepared (in the dark) as per manufacturer protocol. To help avoid aggregates, each individual dextramer reagent was first microcentrifuge at full speed for 5 mins before adding 2 µL from each dCODE™ Dextramer® specificity to a low-bind nucleus free 1.5 eppendorf tube (Eppendorf, 30108051) containing 8.8 µl 100µM d-Biotin (Avidity science, BIO200) (0.2 ul d-Biotin per number of dCODE™ Dextramer® specificity i.e. 44).The dCODE™ Dextramer® master mix was mixed by gently pipetting before the total volume (96.8 µl) was added to the resuspended cells. The sample was then thoroughly mixed and incubated at room temperature for 30 mins in the dark. Following the addition of anti-human CD14-FITC (Biolegend, 325603) and CD3-APC (Biolegend, 300458) (at 1:50) the cells were incubated for a further 20 mins (at room temperature in the dark) before being topped up to 1.4 mL with wash buffer (1X PBS pH 7.4 containing 5 % heat inactivated FBS). The cells were centrifuged down to a pellet (400g for 5 min at 4 °C) and the supernatant discarded. The wash step was then repeated X2, with the latter using the addition of 1.4 mL wash buffer + 1:5000 DAPI (Sigma) as live/dead stain. The supernatant was removed and the cell pellet resuspended in 4 mL FACS buffer ; 1X PBS, 1% FBS, 25mM HEPES (Thermo Fisher Scientific, 15630-056) and 1 mM EDTA. The samples were then filtered (35 µm nylon mesh cell strainer) and PE dCODE Dextramer®-positive cells sorted using a MACSQuant® Tyto® Cell Sorter as per manufacturer guide (Settings; Mix speed= 800 rpm, Chamber temperature= 4 °C, Pressure= 150hPA, Noise Threshold =14.40, Trigger Threshold= off). Note: In order to collect as many cells as possible during sorting the entire sample was run on the MACSQuant Tyto, with the negative run through collected and re-run a second time to ensure no true positives were lost. See Extended Data Fig. 8a for gating strategy for sorting. The PE dCODE Dextramer®-positive cells were then collected, centrifuged (400g for 5 min at 4 °C) and resuspended in resuspension media before counting the cells. The entire sample was then processed for 10x 5′ single cell capture (Chromium Next GEM Single Cell V(D)J Reagent Kit v1.1 with Feature Barcoding technology for cell Surface Protein-Rev D protocol). Where over 25,000 cells were collected the sample was split equally and loaded over two lanes.

In order to provide additional controls, participants with non-compatible HLA types, including one volunteer (674700) matching none of the HLA types for the multi-allele dCODE Dextramer panel, were also processed and used to determine background noise.

Library generation and sequencing

The Chromium Next GEM Single Cell 5′ V(D)J Reagent Kit (V1.1 chemistry) was used for single-cell RNA-seq library construction for all nasopharyngeal swab samples, and the Chromium Next GEM Single Cell V(D)J Reagent Kit v1.1 with Feature Barcoding technology for cell surface proteins was used for PBMCs, both to process the PBMCs stained with CITE-sequencing antibodies panel and the dCODE™ Dextramer® (10x compatible) panel. GEX and V(D)J libraries were prepared according to the manufacturer’s protocol (10x Genomics) using individual Chromium i7 Sample Indices. Additional TCR γ/δ enriched libraries were also generated based upon an in-house protocol previously described in ³⁰. The cell surface protein libraries were created according to the manufacturer’s protocol with slight modifications used for the creation of libraries generated from the CITE-sequencing antibody panel. These included doubling the SI primer amount per reaction and reducing the number of amplification cycles to 7 during the index PCR to avoid the daisy chains effect. GEX, V(D)J and the CITE-sequencing derived cell surface protein indexed libraries were pooled at a ratio of 1:0.1:0.4 and sequenced on a NovaSeq 6000 S4 Flowcell (paired-end, 150 bp reads) aiming for a minimum of 50,000 paired-end reads per cell for GEX libraries and 5,000 paired-end reads per cell for V(D)J and cell surface protein libraries. The dextramer derived cell surface protein indexed libraries were submitted at a ratio of 0.1.

Single cell genomics data alignment

Single-cell RNA-seq and CITE-seq data from PBMCs was jointly aligned against the GRCh38 reference that 10X Genomics provided with CellRanger 3.0.0, and alignment was performed using CellRanger 4.0.0. CITE-seq antibody-derived tag (ADT) barcodes were aligned against a barcode reference provided by the supplier, which we annotated to add informative protein names and made available in our GitHub repository. Single-cell RNA-seq data from nasopharyngeal swab samples were aligned against the same reference using STARSolo 2.7.3a, and post-processed with an implementation of emptydrops extracted from CellRanger 3.0.2. To detect viral RNA in infected cells, we added 21 viral genomes including pre-Alpha SARS-CoV-2 (NC_045512.2) to the above mentioned reference genomes for RNA-seq alignment, as described in Yoshida et al ⁵. Single cell alpha/betaTCR and BCR data was aligned using CellRanger 4.0.0 with the accompanying GRCh38 VDJ reference that 10X Genomics provided. Single cell gamma/delta TCR data was aligned against the GRCh38 reference that 10X Genomics provided with CellRanger 5.0.0, using CellRanger 6.1.2.

Single cell genomics data processing

Both single cell RNA-seq and ADT-seq data were corrected using SoupX ³¹ to remove free-floating and background RNAs and ADTs. To correct ADT counts, SoupX 1.5.2 parameters soupQuantile and tfidfMin parameters were set to 0.25 and 0.2, respectively, and lowered by decrements of 0.05 until the contamination fraction was calculated using the autoEstCont function. SoupX on RNA data was performed using default settings. To confidently annotate SARS-CoV-2 infected cells, we used SoupX corrected viral RNA counts to remove false positives due to freely floating SARS-CoV-2 virions. However, when quantifying the amount of reads per cell in Fig. 2h and their distribution over the viral genome in Fig. 2f, we used the raw counts and sequencing data. To profile the distribution of viral reads, we removed PCR duplicates from the aligned BAM files that STARSolo produced with MarkDuplicates in picard (https://broadinstitute.github.io/picard/), and tallied the location within the SARS-CoV-2 genome using the start of each sequencing read. Aligned single cell RNA-seq data was imported from the filtered_feature_bc_matrix folder into Seurat V4.1.0 for processing, keeping only cells with at least 200 RNA features detected. Nasopharyngeal and PBMC cells with more than 50% and 10% of the counts coming from mitochondrial genes were excluded, respectively. SoupX corrected gene expression and ADT counts were normalized by dividing it by the total counts per cell and multiplying by 10 000, followed by adding one and a natural-log transformation (log1p).

Demultiplexing and patient id assignment

Each PBMC sample was pooled twice into two unique pools containing up to four PBMC samples per pool, followed by CITE-seq and single cell VDJ sequencing as described above. Souporcell V2.0 ³² was used to demultiplex each pools based on the genotype differences between the mixed samples. Souporcell analyses were performed with the skip_remap parameter enabled and using the common SNP database that was provided by the software. We used two complementary approaches to confidently assign participant identity to each souporcell cluster. First we compared the cluster genotypes with SNP array derived genotyping data, generated for all participants and performed using the Affymetrix UK Biobank AxiomTM Array kit by Cambridge Genomic Services (CGS). Second, the combinations of samples within each pool was unique, enabling assignment of participant identity based on the presence of unique participant-specific combinations of identical genotypes in two separate pools. This multiplexing and replication strategy furthermore enabled us to distinguish library specific batch effects from participant specific effects in downstream analyses.

Doublet detection

We used the output from souporcell to identify ground-truth doublets in PBMCs by selecting droplets that contained two genotypes from different participants. We then included these ground-truth doublets into the iterative rounds of subclustering and cell state annotation to look for doublet specific clusters that emerged, which we then subsequently removed. Doublets in the nasopharyngeal data were removed during iterative rounds of subclustering and cell state annotation by identifying cell clusters that expressed marker genes from multiple distinct cell types.

Clustering and cell type annotation

Principal component analysis was run on corrected gene expression counts from selected hypervariable genes and the first 30 principal components were selected to construct a nearest neighbour graph and UMAP embedding. We used harmony³³ to perform batch correction on the PBMC data on the sequencing library identity to remove technical batch effects. Leiden clustering³⁴ performed at resolutions of 0.5, 1, 4 and 32, on nearest neighbour graphs and embeddings created with 500, 1000, 2000, 4000, 6000 and 8000 selected hypervariable genes (excluding TCR and BCR genes), was used to perform iterative rounds of cell type annotation based on marker gene expression and subsetting of clusters to obtain a highly granular cell state annotation. We used cell type marker genes described in Yoshida et al⁵ and Stephenson et al⁴ to define cell types. Our cell type annotation was furthermore guided by predicted cell type labels using Celltypist’s³⁵ provided models and custom trained models based on annotations in Yoshida et al⁵ and Stephenson et al⁴.

Single cell TCR and BCR data processing

Aligned single cell BCR and alpha/beta TCR sequencing data was imported in scirpy³⁶ to obtain a cell by TCR or BCR formatted table, which was then added to Seurat objects containing gene expression data. Aligned single cell gamma/delta TCR data was reannotated using Dandelion V0.2.4³⁷.

Integration of five COVID-19 studies

All transcriptomic data was processed with the single cell analysis Python workflow Scanpy³⁸. Each data set was individually filtered following best practices outlined in³⁹ (Between 200 and 3500 genes per cell, less than 10% mitochondrial genes expressed per cell, genes expressed in at least 3 cells, other parameters at default). The gene sets were reduced to their intersection before combining data sets. Cells came from a total of 602 individuals, with 325 acute COVID-19 patients, 110 COVID-19 convalescent patients, 114 healthy and 53 hospitalized controls (Extended Data Table 1d). This resulted in an integrated embedding containing 946,584 T cells with resolved TCR from 494 samples, made up of 455 donors of which 240 were acute COVID-19 patients, 82 convalescent, 88 healthy and 45 hospitalized controls (Extended Data Table 1e). The total number of donors in the integrated object is smaller as only samples with matching VDJ sequencing data were kept. A probabilistic scVI model (2 hidden layers, 128 hidden nodes, 20-dimensional latent space, negative binomial gene likelihood, other parameters at default⁴⁰) was trained on the data to map cells to a shared latent space, and visualised using UMAP.

Identification of activated TCR clonotype groups using Cell2TCR

To identify TCR clonotype groups, we used tcrdist3⁴¹ with the provided human references to compute a sparse representation of the distance matrices for all identified TRA and TRB CDR3 sequences, with the radius parameter set to 150. We then summed the distances for TRA and TRB to obtain a combined distance matrix. Next, we iterated over possible TCR distance thresholds between 5 and 150 with increments of 5, to compute TCR clonotype groups at each threshold. We then generated a distance adjacency graph of TCRs from different T cells with a distance lower than the threshold, which was clustered to identify TCR clonotype groups using leiden³⁴ clustering through the igraph package⁴², at a resolution of 1 and using the RBConfigurationVertexPartition partition. To find the optimal distance threshold at which only TCRs that recognise the same antigen are grouped together, we quantified clonotype group contamination at each threshold using two approaches. First, we assumed that T cells that were annotated as naive should not participate in an expanded clonotype group, and quantified the proportion of naive T cells in each clonotype group to determine the largest threshold at which we observed minimal participation of naive T cells. Second, we assumed that CD4+ T cells and CD8+ T cells should never be part of the same TCR clonotype group, so we set out to quantify the proportion of CD4+/CD8+ mixing in each clonotype group to find the largest threshold where mixing is minimal. Both approaches revealed the same optimal threshold of 35 at which both naive T cell participation and CD4+/CD8+ mixing is minimal, which we then used for downstream analyses. To identify activated TCR clonotype groups we assumed that these groups should include activated T cells and that we should at least detect multiple independent TCR clonotypes that appear to be raised against the same antigen at the same time; we therefore selected clonotype groups that contained at least one participating activated T cell and that contained at least two unique CDR3 nucleotide sequences.

Identification of activated BCR clonotype groups

To identify BCR clonotypes groups that were activated during infection, we used a similar approach as described above for T cells. Instead of using tcrdist to compute distances, we used the levenshtein distance and iterated over possible thresholds between 1 and 20 to find an optimal threshold by quantifying naive B cell participation. This revealed that a levenshtein distance of 2 is optimal to identify BCR clonotype groups that only contain B cells that recognise the same antigen. To identify activated BCR clonotype groups, we assumed that these groups should include antibody secreting B cells (plasmablasts and plasma cells) and that we should at least detect multiple independent BCRs clonotypes that appear to be raised against the same antigen at the same time; we therefore selected clonotype groups that contained at least one participating antibody secreting B cell and that contained at least three unique CDR3 nucleotide sequences.

Generation of VDJ logos

TCR and BCR logos in Figure X, X and X were generated by providing the CDR3 amino acid sequences of each clonotype group to the ggseqlogo R package⁴³ or the logomaker Python package⁴⁴. When clonotype groups contained CDR3 amino acid sequences of variable lengths, we selected the sequences with the most frequently occurring length within each group for visualization purposes only.

Generalized linear mixed models of cell state compositional changes over time

The relative amount of cells per cell type in each sample was modeled using a generalized linear mixed model with a poisson outcome. When technical replicates were available (most of the PBMC samples), these were modeled as separate samples. We modeled participant identifiers, days since inoculation, and sequencing library identifiers (of multiplexed libraries), as random effects terms to overcome colinearity between these factors. The effect of each clinical/technical factor on cell type composition was estimated by the interaction term with the cell type. The glmer function in the lme4 package implemented on R was used to fit the model. The standard error of the variance parameter for each factor was estimated using the numDeriv package. The conditional distribution of the fold change estimate of a level of each factor was obtained using the ranef function in the lme4 package. The log- transformed fold change is relative to the pre-inoculation time point (day -1). The statistical significance of the fold change estimate was measured by the local true sign rate, which is the probability that the estimated direction of the effect is true, that is, the probability that the true log-transformed fold change is greater than 0 if the estimated mean is positive (or less than 0 if the estimated mean is negative). We calculated P values using a two-sample Z-test using the estimated mean and standard deviation of the distribution of the effect (log-transformed fold change). P values were converted into false discovery rates using the Benjamini & Hochberg method.

Gaussian processes regression and latent variable models to infer time since viral exposure

To infer time from cell state abundances, we first generated a linear regression model using CellTypist to predict PBMC or nasopharyngeal cell states based on the highly detailed manual cell state annotation presented in this work. Celltypist models were trained and used using default parameters, with check_expression set to false, balance_cell_type set to true, feature_selection set to true, and max_iter set to 150. We next set out to build a predictive model to infer time since viral exposure using the PBMC data presented in this work as a training dataset. We used the above mentioned publicly available PBMC data from five studies as a test dataset to predict time since viral exposure on. Because we were specifically interested in comparing time since viral exposure to reported time since onset of symptoms in varying disease severities, we excluded samples for which these features were unknown. To ensure that the cell state proportions in the training and test dataset were similar, we used our CellTypist model on both datasets to predict relative cell state frequencies, which was used as input for our time prediction model. To account for participant-to-participant heterogeneity and continuous variation in the timeline of immune responses, we first constructed a gaussian process latent variable model ⁴⁵ (GPLVM) to smooth the time since viral exposure in the training dataset. We applied the Pyro implementation of GPLVM ⁴⁶ across all predicted cell state abundances, and restricted the model to 2000 iterations and a single latent variable that was initialized on the square root transformed time since inoculation, which resulted in an accurate recapitulation of the mean time since inoculation while smoothing outliers. We next used each predicted cell state as input for a task to generate a multi-task gaussian process regression model⁴⁷ to predict the smoothened time since inoculation using GPyTorch ⁴⁸. We used the Adam optimiser and allowed for as many iterations for the loss in marginal log likelihood to reach zero. We next predicted the cell state compositions across the entire tested timeline (day -1 to day 28), and compared these cell state compositions to the cell state compositions in our query dataset as predicted by our CellTypist model. Last, we selected the time point whose predicted cell state composition had the lowest mean squared error compared to the observed cell state composition.

Matching clonotype groups to antigen-TCR database

Computation of fold change enrichment of SARS-CoV-2 specific TCRs in Activated T cell populations compared to other T cell populations: Median p for 10 random draws of n=5000 unique clones of both populations = 3.75, p = 1.68 × 10^-21

Bulk TCR sequencing and processing

Total RNA was extracted from whole blood samples collected in Tempus Blood RNA tubes (Thermofisher #4342792) using the manufacturer’s protocol for RNA extraction. TCR α and β genes were sequenced using a pipeline which introduces unique molecular identifiers attached to individual cDNA molecules using single-stranded DNA ligation. The unique molecular identifier allows correction for sequencing error PCR bias, and provides a quantitative and reproducible method of repertoire analysis. Full details for both the experimental TCRseq library preparation (^{49, 50}) and the subsequent TCR annotation (V, J and CDR3 annotation) using Decombinator V4⁵¹ are published. The Decombinator software is freely available at https://github.com/innate2adaptive/Decombinator.

Memory formation analysis

T cell phenotypes (naive, activated, effector, memory) were recorded for an antigen-specific TCR clone at different time points throughout infection. TCR clones were filtered by having an activated label at least once, being observed in at least two samples, one of which had to be at day 28. Unique TCR clones are distinguished by color and numbered with their clone_id identifier. A shaded area is drawn when the same clone appeared with several distinct cell type labels, and the size of the shaded area informs of their relative ratios.

Quantifying TCR diversity restriction in phenotypic clusters using coincidence analysis

To quantify the diversity of TCRs found within different phenotypic clusters we determined the probability with which two distinct clonotypes within a cluster share an identical CDR3 amino acid sequence ⁵². For visualization we normalized these probabilities by the same quantity calculated over the complete data regardless of phenotype. This ratio of probability of coincidences provides a stringent measure of convergent functional selection of distinct clonotypes that share the same TCR. The analysis is based on clonotypes defined by distinct nucleotide sequences of the hypervariable regions, and does not make direct use of clonal abundances as these can also reflect TCR-independent lineage differences. We focused our analysis on conventional T cells only, considered only cells with at least one valid functional alpha and beta chain, and kept only a single chain for each cell where there were multiple chains. We performed the analysis both on the alpha and beta chain separately, as well as on paired alpha and beta chain, in each instance requiring exact matching of the CDR3 amino acid sequences.