Clone tracking through repeated malaria identifies high-fidelity CD4+ memory responses - PMC

Longitudinal scRNA/TCRseq enables the study of human CD4+ T cell memory in vivo

To study human CD4+ T cell memory in vivo in an unbiased way, we performed single-cell RNA and TCR sequencing on CD45RA- memory CD4+ T cells before, during, and at multiple timepoints after repeated episodes of malaria in Ugandan children (Fig. 1A). Prior to sequencing, cells were either left unstimulated or stimulated in vitro with αCD3/αCD28 coated beads for 6 hours or stimulated with Plasmodium-infected red blood cells (iRBCs) for 9 hours (Fig. 1A, Fig. S1). Our study included 34 samples collected longitudinally in the context of 13 infections across 8 children (ages 4–10) and 2 clinically immune adults (Fig. 1B, Fig. S2). By performing longitudinal sampling in this way, we aimed to track CD4+ T cell clones through multiple infections to characterize the longevity, functional stability, and memory potential of Pf-specific CD4+ T cells.

In total, we sequenced 504,041 high-quality CD4+ T cells and recovered TCR sequences for 450,078 of these cells. Using these data, we generated gene correlation matrices to identify gene modules (Fig. S3) and then annotated clusters based on gene module scores (Fig. S4). In doing so, we identifed 11 transcriptionally distinct memory CD4+ T cell populations from cells sequenced across all three in vitro stimulation conditions (Fig. 1C, Fig. S5). These annotated populations included T helper subsets (Th1, cytotoxic Th1, Th2, Th17, Tcm), regulatory subsets (Treg, Tr1), naïve cells (that either did not express CD45RA or represented a contaminating population), and 3 heterogeneous populations defined by their response to interferon stimulation, TCR activation, or proliferative signaling (Fig. 1C). The T helper subsets displayed transcriptional signatures similar to what has previously been described(30, 31), and Tr1 cells (~3% of CD45RA-CD4+ T cells) displayed characteristics of both Tregs (expression of PRDM1, CTLA4, PDCD1) and Th1 cells (expression of CXCR3, TBX21) (Fig. 1D). We observed distinct patterns of granzyme expression between Tr1 (GZMA, GZMK) and cytotoxic Th1 cells (GZMB, GZMH) (Fig. 1D). Tcm cells accounted for ~50% of CD45RA-CD4+ T cells, while the various effector memory populations accounted for 1–16%, with Th2 and Th17 most abundant (Fig. S6A). We did not observe a distinct population of circulating Tfh, and cells that expressed CXCR5 transcripts were mostly observed among Tcm cells (Fig. S6B). Additionally, Tr1 cells did not seem to express ITGA2 (encoding CD49b), even though prior studies have used CD49b (when co-expressed with LAG-3) to enrich for Tr1 cells in flow cytometry-based studies (Fig. S6C). In sum, we identify, characterize, and quantify the relative abundances of transcriptionally distinct CD4+ T cell subsets among Ugandan children exposed to malaria.

Seven CD4+ T cell subsets display clonal fidelity

In vivo studies employing single-cell genomics to study the plasticity and intraclonal heterogeneity of human CD4+ T cell clones over time are lacking. One prior study employing bulk-TCR beta-chain sequencing suggested a high degree of intraclonal diversity such that CD4+ T cells with the same TCR can adopt multiple different fates(32). In contrast to this, we observed very large clonotypes that seemed to display strong preferences for a single CD4+ T cell subset (Fig. 2A). To explore this further, we performed hierarchical clustering of expanded TCR clonotypes in our dataset (those that consisted of at least 5 clones) based on their distribution across the resting CD4+ T cell subsets (excluding cells in the “TCR-Activated”, “IFN-Stimulated”, and “Proliferating” clusters) (Fig. 2B). Strikingly, nearly all clonotypes displayed a strong preference for one of seven subsets, suggesting subset-specific fidelity among CD4+ T cell clones. This hierarchical clustering, thus, yielded 7 families, each displaying a certain degree of clonal fidelity. Clonotypes from family 6 displayed the highest clonal fidelity, as 98% of clones among these family 6 clonotypes on average were Tregs (Fig. 2C), consistent with naïve Treg clones diverging early from naïve T helper cells during thymic selection(33). Surprisingly, Tr1 cells also displayed high average clonal fidelity (0.93) and relative independence from Th1 cells (Fig. 2C), even though prior studies have suggested that the Tr1 state may represent a transient or terminal functional program of Th1 cells(34, 35). Distinct clonal communities were observed within Th1 cells, as clonotypes displayed some preference for either the cytotoxic (family 5) or non-cytotoxic (family 3) subsets of this lineage (Fig. 2B-C). Importantly, many clonotypes could be detected at multiple timepoints and retained their fidelity for hundreds of days (Fig. 2D). Only 8 of the 5,770 clonotypes assigned to a family were expanded in two different study participants. 7 of these 8 clonotypes adopted the same fate despite deriving from different individuals (Table S1), suggesting a role for TCR specificity in directing T cell differentiation.

Clonal fidelity as well as clonotype size varied across the seven families. Family 5 (cytotoxic Th1) clonotypes were the largest, and family 6 (Treg) clonotypes were the smallest (Fig. 2E). Thus, we wondered how clonal expansion influences fidelity. Assuming no relationship, the probability, P, that a clonotype displays a fidelity of 1 can be modelled as P = fⁿ, where f is the average clonal fidelity and n is the clonotype size. So, given that family 1 (Th2) clonotypes displayed an average fidelity of 0.93 (Fig. 2C), the probability that a family 1 clonotype with 64 observed clones would be entirely Th2 is: P = 0.93⁶⁴ = 0.01. However, we observed many phenotypically homogenous family 1 clonotypes of sizes exceeding 64 clones. In fact, our data demonstrated that the probability of this specific event (where n=64) is likely greater than 0.5 (Fig. 2F). A similar phenomenological relationship was observed among family 5 (cytotoxic Th1) clonotypes as well (Fig. 2G), but not family 2 (Tcm) clonotypes (Fig. 2H). Thus, clonal fidelity increases substantially with clonal expansion, suggesting that effector polarization and proliferation signals may be linked in vivo.

Clone tracking identifies subset specific activation programs

Given the apparent clonal fidelity of CD4+ T cells, TCR sequences can be used to link diverse activation programs with the well-defined resting states that precede them (Fig. 3A). More precisely, once clonotypes have been grouped into families based on their relative frequencies across different resting subsets, clones can be tracked from the resting populations to any activated clusters, facilitating the transfer of family-level, cell subset annotations (Fig. 3A). Comparing gene expression across activated and resting populations can reveal transcriptional programs of activation that are broadly conserved and subset specific. This approach improves upon manually annotating activated clusters based on transcriptomics because it is unbiased and does not rely on a priori knowledge of subset-specific activation genes (i.e., that Th1 cells should express IFNG). TCR sequences are intrinsically related to T cell ontogeny, and thus, our approach can identify subset specific activation programs with high confidence and minimal bias.

As expected, nearly all cells that were stimulated in vitro with αCD3/αCD28 coated beads clustered together in the “TCR-Activated” population (Fig. 3B), confirming that stimulation in this way induced activation of all CD4+ T cell subsets. On the other hand, many cells were not activated in response to in vitro stimulation with iRBCs, and the T cells that did become activated were less transcriptionally diverse than bead-activated cells (Fig. 3B). Using clone tracking as described above to follow clones from a resting to an activated state, we found that each clonotype family occupied a distinct region of the “TCR-Activated” cluster (Fig. 3C). Tracking clones in this way, we validated prior findings and uncovered new biology with respect to subset-specific/skewed TCR activated genes (TAGs) (Fig. 3D). For example, family 1 (Th2) clonotypes responded to TCR activation by expressing canonical Th2 lymphokines such as IL4, IL5, IL9, IL13, and IL24; however, these cells were also the primary expressors of IL3 and CSF2, important hematopoietic cytokines that have not previously been attributed specifically to Th2 cells (Fig. 3D). Additionally, family 7 (Tr1) clonotypes expressed IL10 and IL21 upon TCR activation, which is consistent with prior literature(36); however, these cells also upregulated many other genes that have not been previously attributed to Tr1 cells (Fig. 3D). For example, they upregulated TNFRSF1B, the gene encoding tumor necrosis factor receptor 2 (TNFR2), more than any other subset. Previously thought to be specific to FOXP3-expressing Tregs, this receptor binds TNF and transduces a signal that supports survival and regulatory function(37). Interestingly, family 3 (Th1) and family 5 (cytotoxic Th1) cells responded to TCR activation by upregulating CCL3 and CCL4 (Fig. 3D), chemotactic ligands for CCR5, which is most abundantly expressed by Tr1 cells(38) (Fig. 1D), suggesting that Tr1 cells may exhibit preferential chemotaxis towards activated Th1 cells. While murine Tregs express IL10 upon activation(39), we saw little evidence of this among human Tregs activated in vitro (Fig. 1D). Family 6 (Treg) did, however, specifically upregulate IL21R (which may help regulate their proliferation(40) and survival(41)) and BATF (which was recently shown to be indispensable for the activation of tumor associated Tregs(42)). In general T helper subsets upregulated TNF, IL2, and FASLG, while regulatory subsets expressed much less TNF and IL2 and instead upregulated FAS.

Clone tracking was also used to elucidate general and subset-specific interferon stimulated genes (ISGs) (Fig. 3D). Canonical ISGs such as OAS1, IFIT1, IRF7, ISG15, and MX1 were upregulated across all subsets, while expression of other genes (IRF3, CCR1, GSN, STAT1) were preferentially induced in some subsets. Ultimately, this dataset describing subset specific activation programs of CD4+ T cells provides a much-needed reference for the transcriptional programs defining diverse CD4+ T cell subsets and states.

Tr1 cells dominate the Pf-specific CD4+ T cell response

Knowledge of the Pf-specific human CD4+ T cell repertoire is extremely limited. Here, CD4+ T cells of Ugandan children and adults responded transcriptionally to in vitro stimulation with iRBCs in an antigen-specific and non-specific, interferon-mediated manner (Fig. S7). And, while in vitro stimulation with αCD3/αCD28 coated beads led to TCR activation of all CD4+ T cell subsets, stimulation with iRBCs specifically activated only family 7 (Tr1) clonotypes (Fig. 4A). We identified Pf-specific TCRs from clonotypes that saw an increase in frequency of TCR-activated clones (by at least 0.4, or 40%) following iRBC stimulation (Fig. 4B). Of the 39 Pf-specific clonotypes that we identified, 29 were observed at high enough frequencies to be assigned to a clonotype family. Strikingly, all of these 29 clonotypes belonged to the Tr1 family (family 7) (Table S2). So, while Tr1 cells accounted for ~3% of memory CD4+ T cells with any specificity in Ugandan children, these cells made up nearly 90% of CD4+ T cells specific to Pf blood-stage antigens (Fig. 4C). In contrast, malaria-naive adults from North America had very few Tr1 cells (as well as reduced Th2 frequencies) (Fig. S8A-B). Furthermore, when stimulated with iRBCs in vitro, very few CD4+ T cells of malaria-naïve adults became activated (Fig. S7A), and those that did appeared to be Tcm/naïve, as they upregulated TNF and IL2 but not other subset-specific cytokines (Fig. S8C-D). Compared to Ugandan children, the CD4+ T cell response of malaria-naïve adults to iRBCs was more inflammatory and did not include Tr1 cells (Fig. S8D).

Importantly, the Ugandan Tr1 response to iRBCs is transcriptionally very similar to the Tr1 response to αCD3/αCD28 stimulation. Under both stimulation conditions, family 7 (Tr1) clones upregulated key genes downstream of TCR activation, including IL2RA, TNFRSF9, TNFRSF4, BCL2A1, NFKB1, RELB, IRF4, and MYC (Fig. 4D-E). In response to in vitro stimulation with uninfected red blood cells (uRBCs), memory CD4+ T cells did not produce canonical Tr1 cytokines (IL-10, IFNγ) as assessed by flow cytometry (Fig. S9). Furthermore, blocking HLA class II molecules in PBMC cultures completely abrogated IL-10 and IFNγ production by memory CD4+ T cells stimulated with iRBCs (Fig. S10). Thus, the Tr1 response to Plasmodium is antigen-specific and not induced by mitogens.

We observed reliable Tr1 expansion and contraction in the setting of repeated symptomatic malaria across multiple children as measured using scRNAseq (Fig. 4F, Fig. S11) and flow cytometry (Fig. S12). Following some infections, Tr1 frequencies more than tripled, and these cells accounted for nearly 10% of circulating CD45RA-CD4+ T cells within a week of diagnosis (Fig. 4F). Tr1 abundance declined following malaria but may persist at higher frequencies among clinically immune individuals (Fig. 4G). Lastly, parasitemia (measured by qPCR) at the time of diagnosis correlated with changes in Tr1 cell frequencies (Fig. 4H), suggesting that antigen load is a likely determinant of Tr1 expansion. In sum, the Pf blood stage-specific CD4+ T cell response appears to be dominated by Tr1 cells and not IL-10-producing Th1 cells.

Tr1 cells demonstrate long-term memory with functional stability

To further study CD4+ T cell memory in vivo, we leveraged the longitudinal aspect of our data and tracked the dynamics of Pf-specific clonotypes (the N=39 clonotypes that responded to in vitro stimulation with iRBCs) over the course of repeated infections. The Pf-specific clonotypes from participants 3178, 3410, and 3354 could be observed at multiple timepoints (Fig. 5A). Following symptomatic malaria, these Pf-specific Tr1 clones reliably underwent clonal expansion. The same clonotypes that expanded in response to the first study infection also frequently expanded in response to subsequent repeated infections (Fig. 5A), further confirming their specificity for Pf antigens and demonstrating their memory potential. Moreover, repeat infections seemed to boost clonal expansion of Pf-specific Tr1 cells, especially in participants 3178 and 3410, where the second infection confirmed during this study occurred hundreds of days after the first. Although non-Tr1 family clonotypes did not respond to in vitro stimulation with Pf blood-stage antigens in Ugandan children, we wondered whether other clonotypes might reliably expand upon repeated infection, such as those clonotypes specific for liver-stage antigens. Thus, we selected clonotypes of the Tcm, Th1, cytotoxic Th1, Th17, Th2, and Treg families that expanded the most following the first study infection of participant 3354 and then asked whether these clonotypes would also expand in response to two subsequent infections. They did not (Fig. S13A-F), suggesting that the initial expansion of these non-Tr1 clones was unrelated to malaria and that they are not Pf-specific.

Our observations that (1) Tr1 cells are clonally distinct from Th1 cells and (2) Pf-specific Tr1 clonotypes can reliably respond to repeated natural infections together suggest that Tr1 cells display long-term memory. We, thus, sought to determine whether Tr1 function might be epigenetically encoded. Using multiome sequencing, we quantified chromatin accessibility and gene expression in individual memory CD4+ T cells collected from pediatric participant 3394 at three distinct timepoints before and after malaria. Cells were annotated by using our scRNA/TCRseq (5’ chemistry) dataset as a reference for mapping cell identity within our 3’ multiome dataset (3’ chemistry) (Fig. S14A-C). Examining the epigenetic landscape of Tr1 cells, we found that regions around the loci of key Tr1 genes such as IL10 and LAG3 were more accessible compared to other CD4+ T cell subsets (Fig. S14D,E). Additionally, peaks near the IFNG locus were accessible among Th1, cytotoxic Th1, and Tr1 cells (Fig. S14F). To explore the molecular regulators that are potentially relevant to the Tr1 versus Th1 fate decision, we estimated transcription factor activity using chromVAR. Several transcription factors appeared to be differentially active between the two populations, with Th1 cells displaying greater accessibility around the binding motifs of AP-1/ATF superfamily members (FOS, JUN, BATF) and T-box binding factors (TBX21, EOMES) and Tr1 cells displaying greater accessibility around binding motifs of CTCF (a zinc finger protein with diverse roles in gene regulation and chromosomal architecture)(43) and IRF family proteins (potentially consistent with prior suggestions for a role of type I interferon in Tr1 differentiation)(44, 45) (Fig. S14G). Similar to other CD4+ T cell subsets, TCR activation of Tr1 cells corresponded with increased activity of members of the AP-1/ATF superfamily as well as the NFAT and NF-κB families (Fig. S14H). Though these data should be considered preliminary as they are from three longitudinal samples taken from a single individual, they suggest an epigenetic basis for Tr1 function that should be validated in future studies.

We next explored Tr1 transcriptional heterogeneity to uncover how Tr1 memory responses might be sustained. Reclustering Tr1 cells, we identified heterogeneous populations (Fig. 5B). Recently activated Tr1 cells expressed IL10 and IFNG (Fig. 5B-D) and were primarily observed in samples that were stimulated in vitro with either αCD3/αCD28 coated beads or iRBCs (Fig. S15). Putatively defined effector Tr1 cells expressed IL10 but minimal IFNG (Fig. 5B-D) and were abundant in many unstimulated samples (Fig. S15). Crucially, activated Tr1 cells demonstrated a transcriptomic signature consistent with recent TCR activation—effector Tr1 cells did not appear recently activated, despite their production of many effector molecules. We also identified a putative population of memory Tr1 cells that expressed minimal IL10 and IFNG but expressed high levels of KLF2 (Fig. 5B-E), an important regulator of T cell survival and quiescence(46). Differential gene expression highlighted key differences between effector and memory Tr1 cells (Fig. 5F), and some of these differences similarly characterized effector versus memory CD8+ T cells (Fig. 5G), suggesting conserved transcriptional networks broadly mediating T cell memory. Importantly, memory Tr1 cells were distinct from Th1 cells (Fig. S16). We also identified a rare population of naïve-like Tr1 cells that may represent an actively differentiating population (Fig. 5B). Compared to effectors, these naïve-like Tr1 cells appeared more quiescent with reduced expression of many metabolic genes including those encoding ribosomal proteins, electron transport proteins, and ATP synthases (Fig. 5H). These differences observed between naïve-like and effector Tr1 cells, once again, mirrored the differences observed between naïve and effector CD8+ T cells (Fig. 5I).

Transcriptomics alone may not conclusively demonstrate that memory and effector Tr1 subsets are functionally distinct. Therefore, we modelled the dynamics of these subsets in the context of malaria (Fig. 5J). Naïve-like, effector, and memory Tr1 cells rapidly expanded during symptomatic infection, peaking around 14 days. Contraction of effector and memory Tr1 cells ensued in the months following diagnosis and antimalarial drug treatment. While Tr1 effectors disappeared after 150 days, Tr1 memory cells persisted in circulation for over 300 days, as long as we followed them (Fig. 5J). These dynamics were also assessed at the clonal level in participant 3481, since this individual was sampled extensively in the context of a single infection. Shortly before an infection, participant 3481 had a majority of circulating Tr1 memory cells (Fig. 5K). At diagnosis, there was an early expansion of activated and effector Tr1 cells, which was likely initiated from the preexisting memory pool (Fig. 5K). By day 6, there was the appearance of naïve-like Tr1 cells and the prolific clonal expansion of effectors (Fig. 5K). Clonal relationships began to emerge between the newly expanded effectors and memory Tr1 cells by day 28 (Fig. 5K). Finally, Tr1 contraction and memory consolidation began to be observed on days 56 and 84 (Fig. 5K). Together these data demonstrate Tr1 persistence following infection and the memory potential of individual Pf-specific Tr1 clones.

Study Design

This study included samples provided by Ugandan children and adults living in Tororo or Busia who were followed longitudinally through the East African International Centres of Excellence in Malaria Research (ICEMR) cohorts. In this setting, malaria transmission is year-round, with two seasonal peaks. Screening and enrollment of households in the parent ICEMR cohort, including inclusion and exclusion criteria, have been previously described(23, 60). Briefly, households in a contiguous study area across Busia and Tororo districts were enumerated to provide a sampling frame, and 80 households were randomly selected for screening and enrolled if they met the following criteria: 1) having at least two members aged 5 years or younger; 2) no more than 7 permanent residents currently residing; 3) no plans to move from the study catchment area in the next 2 years; and 4) willingness to participate in entomological surveillance studies. Cohort participants were followed for all care at dedicated study clinics. Participants were monitored for parasitemia by microscopy and qPCR by active (monthly) and passive surveillance, with peripheral blood samples taken. Those who presented with a fever (tympanic temperature >38.0 °C) or history of fever in the previous 24 hours had blood obtained by finger prick for a thick smear. If the thick smear was positive for Plasmodium parasites, the patient was diagnosed with malaria regardless of parasite density and was treated with artemether-lumefantrine. Participants with asymptomatic parasitemia were not treated with antimalarial drugs in accordance with local guidelines. In the MUSICAL subcohort, study participants diagnosed with symptomatic malaria were asked to return for repeat clinic visits and blood sampling on days 7, 14, and 28, post-treatment. Those with asymptomatic parasitemia were asked to return for repeat clinic visits and blood sampling on days 14 and 28 following diagnosis. Participant characteristics of samples analyzed are provided in Data File S1. As additional malaria-naïve controls, samples were obtained from anonymous adult blood donors from the Stanford Blood Center (California, USA).

Written informed consent was provided by all participants or by their parents/legal guardians. The study protocols were approved by the Uganda National Council of Science and Technology (HS-2700), the Makerere University School of Medicine Research and Ethics Committee (2019–134), the University of California, San Francisco Committee on Human Research (19–28606), and the Institutional Review Board at Stanford University (IRBs 41197, 52977).

PBMC Isolation

At select study visits, 3 to 10 mls of blood were obtained in acid citrate dextrose tubes. Plasma was removed and peripheral blood mononuclear cells (PBMC) were isolated by density gradient centrifugation; these samples were then cryopreserved and shipped by liquid nitrogen to Stanford University for further analysis. Analysis of cell viability consistently demonstrated >90% viability after thaw.

Parasite Culture and Isolation

Plasmodium falciparum blood-stage 3D7 parasites were grown by standard methods and harvested at 5–10% parasitemia. Red blood cells infected with mature asexual stages were cultured in RPMI 1640 with 25mM HEPES, 25mM sodium bicarbonate, 1% gentamycin, and enriched with 0.5% Albumax (pH 6.75) and 250uM hypoxanthine, under atmospheric conditions (5% oxygen, 5% carbon dioxide, and 95% nitrogen). To retain synchronous cultures, parasite growth was treated with 5% D-sorbitol. Mycoplasma contamination was assayed using the MycoAlert Mycoplasma Detection Kit (Lonza). Schizont isolation (iRBC) was performed magnetically using MACS cell separation LD columns (Miltenyi Biotec), and iRBC stored in liquid nitrogen prior to use. Uninfected RBCs (uRBC) were used as additional controls. For experiments, iRBC or uRBC were thawed, washed and resuspended in RPMI before addition to cells.

Kinetics of cytokine gene expression by qPCR

Frozen PBMC samples from three Ugandan donors were thawed and then either rested or stimulated in vitro with iRBCs for 0, 8, 16, 24, and 32 hours. Then, memory CD4+ T cells (CD3+ CD4+ CD45RA-) were sorted into TRIzol and RNA was isolated according to the manufacturer’s instructions. The SuperScript^™ III Platinum^™ SYBR^™ Green One-Step qRT-PCR Kit was used in conjunction with the StepOnePlus^™ Real-Time PCR System to quantify the abundance of IL10 and IFNG transcripts in the RNA samples. Rox dye was used as a passive reference. Three replicates were assayed for each gene target-sample combination as part of a singleplex experiment performed in a 96-well plate. ΔC_T was calculated as the difference between C_T for target genes (IL10 and IFNG) and the houskeeping gene RPL13A. ΔΔC_T for unstimulated and stimulated samples incubated for various durations was calculated with respect to the T₀ sample (no stimulation; 0 hours).

Single-cell RNA and TCR sequencing (scRNA/TCRseq)

Frozen PBMC samples were thawed and allowed to rest at 37°C for 5 hours. Cells from each sample were then split across three conditions. One third of the PBMCs were incubated with iRBCs at a ratio of 2:1 (PBMCs:iRBCs) for 9 hours at 37°C. Another third of the PBMCs were incubated concurrently in the absence of any stimulation. Finally, memory CD4+ T cells were isolated from the remaining fraction using MACS, and then, these cells were stimulated with αCD3/αCD28 beads at a ratio of 1:1 for 6 hours at 37°C. After incubation, memory CD4+ T cells were isolated from the unstimulated and iRBC-stimulated fractions using MACS, and αCD3/αCD28 beads were removed from the third fraction. MACS of memory CD4+ T cells yielded a relatively high purity population of CD45RA- CD4+ T cells for sequencing (Fig. S18). Cell hashing was performed so that cells from each of the three conditions could be identified using via distinct molecular barcodes after pooling stimulation conditions for single-cell capture.

Approximately 11,000 cells from each stimulation condition were loaded onto a chromium controller such that each well contained a total of 33,000 cells (targeting 20,000 cells). Then single-cell capture and the subsequent preparation of RNA, TCR, and hashtag libraries was performed according to 10X Genomics published protocol (CG000330: Chromium Next GEM Single Cell 5-v2 Cell Surface Protein UserGuide RevD). Paired-end sequencing of the prepared libraries was performed using either the NovaSeq X Plus or NovaSeq 6000 platforms, followed by demultiplexing of pooled libraries.

scRNA/TCRseq data processing

Raw FASTQ files were used to generate filtered feature-barcode matrices and TCR annotations following Cell Ranger Count v7.1.0 (reference Human GRCh38 2020-A) and Cell Ranger V(D)J v7.1.0. These data were further processed and analyzed using Seurat v4.2.1 built under R v4.1.2. Our full data processing and analysis pipeline is described in a flowchart (Fig. S19). Code used for all analyses and visualizations present in this manuscript is available on GitHub (jason-nideffer/clone-tracking-in-malaria). To summarize, the filtered feature-barcode matrices from each sample were loaded into Seurat and TCR annotations were integrated with these data using scRepertoire. Centered log ratio (CLR) transformations were applied to the hashtag data, and then hashtag expression cutoffs were used to exclude multiplets and assign cells to stimulation conditions. Data were additionally filtered such that all retained cells had a minimum of 500 features, a maximum of between 2,500 and 5,000 features (depending on sample sequencing depth), and a maximum of 8% mitochondrial reads.

Annotating scRNAseq data

Seurat objects from all samples (excluding those from participant 3394) were then combined; gene expression data were normalized; and, variable features were identified. Some hallmark genes of CD4+ T cell subsets were added to the list of variable genes, and TCR genes were removed prior to scaling these features. Then, to identify gene modules, correlation and distance matrices were generated using the scaled data for all variable genes, and hierarchical clustering was performed using Ward’s method. Gene modules were established as groups of genes with correlated patterns of expression that also included canonical CD4+ T cell subset-defining genes (i.e., GATA3 for Th2 cells) (Fig. S3). Then, the average expression values of the genes within each module were calculated to provide each individual cell with a score for each module. After, the scaled variable features data were used to perform principal component analysis and Louvain clustering, the mean module scores were calculated on a per-cluster basis and then used to assign cell identity to each cluster (Fig. S4). Following annotation, UMAP was performed for visualization.

Batch effects were observed between samples from participant 3394 and the other samples, and for this reason, samples from participant 3394 were excluded from initial clustering, annotation, and visualization. To overcome these batch effects, we used Seurat to transfer annotations and embeddings via reference mapping, where the query data from participant 3394 were mapped using the rest of our data as a reference.

Identifying clonotype families and assessing clonal fidelity

To analyze the relationship between phenotypes and clonotypes within our data, we developed the R package Phenoclone. Phenoclone was first used to calculate the abundance of each clonotype within samples and the average abundance of clonotype across samples from the same individual. Then clonotypes were hierarchically clustered based on their phenotypic composition—for this, Phenoclone computed a list of phenotypes for every clonotype where phenotypes were ranked based on their abundance within the clonotype. For this analysis, cells belonging to the TCR-activated, IFN-stimulated, or proliferating clusters were excluded (since these clusters represent highly heterogeneous populations). Then, Phenoclone computed distances between clonotypes using a rank-based distance metric that is weighted such that rank dissimilarities contribute to greater distances between clonotypes the earlier that these dissimilarities occur within the list. Clustering was performed using Ward’s method. As a result, clonotypes were grouped, first, according to the most abundant phenotype, then according to the second most abundant phenotype, and so forth.

Clonal fidelity was calculated as the maximum proportion of a phenotype within a clonotype. This was calculated either by aggregating clonotype frequencies from the same individual across multiple timepoints (Fig. 2B-C), or it was calculated independently for each timepoint (Fig. 2D). To assess the relationship between clonal expansion and fidelity, we first calculated the average fidelity of clonotypes belonging to a given family (Fig. 2C). Then, we used this value to model the theoretical probability that a clonotype of size N would have a fidelity of 1, assuming no relationship between clonotype size and fidelity. This function was then compared to a binomial regression model fit using our observed data, where clonotype size served as the independent variable and the binary outcome of whether or not the clonotype had a fidelity of 1 served as the dependent variable.

Identifying activation states using clone tracking

The subset identities of cells within heterogeneous clusters (TCR-activated, IFN-stimulated, and proliferating clusters) were inferred through clonal relationships between cells observed in these clusters and cells observed in the more homogeneous resting clusters. For example, cells belonging to the Th2 clonotype family that were observed in the TCR-activated cluster were assumed to be TCR-activated Th2 cells. We calculated the mean expression of all genes by cells from the 7 clonotype families when those cells were observed in either resting, TCR-activated, or IFN-stimulated clusters. This analysis excluded cells that could not be assigned to a clonotype family. Genes displayed (Fig. 4D) represent genes that were hand-selected because they appeared to demonstrate some subset-specific bias in their upregulation among cells that were TCR-activated or IFN-stimulated.

Identifying Pf-specific clonotypes

We identified Pf-specific clonotypes as those clonotypes that responded transcriptionally to in vitro stimulation with iRBCs. This was only assessed for clonotypes that had at least 3 clones observed in the unstimulated condition and at least 3 clones observed in the iRBC-stimulated condition. Additionally, a clonotype was only deemed malaria specific if the frequency of TCR-activated clones in the iRBC-stimulated condition was higher than the frequency of TCR-activated clones in the unstimulated condition by a value of at least 0.4. This cutoff was selected conservatively to minimize type 1 error—since no clonotype had a frequency of TCR-activated clones that exceeded 0.4 in the unstimulated condition.

Characterizing Tr1 cell heterogeneity

Tr1 cells used in the analyses assessing Tr1 heterogeneity (Fig. 5B-K) were the same cells initially annotated as Tr1 (Fig. 1C), with the addition of a sub-cluster (cluster 50, Fig. S4) within the TCR-activated population (Fig. 1C) that displayed a Tr1 signature and that was clonally related to the other Tr1 cells. This subset of cells was used to rediscover variable features (excluding TCR genes), which were rescaled and used to compute principal components for downstream clustering and UMAP visualization. Clusters were manually combined to yield the four annotated Tr1 populations—activated, effector, memory, and naïve-like. These annotations were supported by performing differential gene expression and gene set enrichment using Seurat and fgsea v.1.20.0, respectively.

Characterizing Tr1 responses via flow cytometry

To characterize Tr1 responses by flow cytometry, PBMCs were stimulated with media, P. falciparum infected RBCs, and/or uninfected RBCs and CD4 T cell production of IFNγ and IL10 were measured via intracellular staining as previously described(25). Briefly, PBMCs were thawed using standard methods, and a total of 10⁶ PBMCs were stimulated with media, intact purified trophozoite/schizont-stage P. falciparum (clone 3D7)–infected RBCs, and/or uninfected RBCs at an effector to target ratio of 1:2. For HLA class II blocking experiments PBMC were additionally stimulated for 24-hours in the presence of 1) HLA class II blocking antibodies (10 μg/ml clone Tü39 [BioLegend 361702], 10 μg/ml clone SPVL3 [Beckman Coulter IM0416], 10 μg/ml clone B7/21 [Thermo Fisher Scientific H260–1MG], and 10 μg/ml clone L243 [BioLegend 307602]), or 2) isotype control antibodies (10 μg/ml mouse IgG3 [Thermo Scientific MG300] and 30 μg/ml mouse IgG2a,κ [BioLegend 400202]). In all experiments, following 6 hours of stimulation, Brefeldin A and monensin (BD Pharmingen) were added (10 μg/mL). At 24 hours, cells were washed, and surface and intracellular staining was performed with antibodies from one of two panels. The first panel included: from BioLegend, anti-CD14-BV510 (M5E2), anti-CD19-BV510 (HIB19), anti-CD45RA-BV605 (HI100), anti-CD4-BV650 (RPA-T4), anti-CD8-BV785 (RPA-T8), anti-IFNγ-PerCP/Cy5.5 (4S.B3); from BD Pharmingen, anti-IL-10-PE (JES3–19F1); from BD, anti-CD3-APC-H7 (SK7); and from Invitrogen, Live/Dead aqua amine. The second panel included: from BD, anti-CD4-BUV395 (SK3), anti-CD45RA-BV785 (HI100), anti-FOXP3-BB700 (236A/E7); from BioLegend, anti-CD8-BV650 (SK1), anti-IL-10-PE/dazzle594 (JES2–19F1), anti-CD3-AF700 (OKT3); and from Thermo Fisher Scientific, Live/Dead Near IR. Samples were acquired on an Attune NxT Acoustic Focusing Cytometer or a BD FACSymphony A5 and analyzed by FlowJo (TreeStar). Cytokine production was quantified among memory CD4+ T cells derived from the various experimental conditions. Paired, two-sided T tests with multiple hypothesis test correction (Benjamini Hochberg method) were employed to assess statistically significant differences between conditions.

Multiome sequencing and data integration

Frozen PBMC samples were thawed and stimulated in vitro with iRBCs as performed for scRNA/TCRseq. After stimulation, single nuclei suspensions were prepared, GEMs were generated, and cDNA libraries were prepared following guidelines provided by 10X Genomics. Sequencing was performed as described for scRNA/TCRseq. Read alignment was performed, once again, using the Cell Ranger pipeline. Features and cells were filtered using Seurat. Then, the gene expression data was used to transfer embeddings and annotations from the scRNA/TCRseq dataset (5’ chemistry) to the multiome dataset (3’ chemistry). We were able to recapitulate the general structure of our UMAP (Fig. S14A). However, to ensure that annotations captured epigenetically distinct populations, we additionally clustered the multiome data based on chromatin accessibility (Fig. S14B) and then assigned identity to clusters based on the relative abundance of transferred annotations within each cluster (Fig. S14C). Finally, comparative epigenetics analyses were carried out using Seurat.

Data visualization and schematics

Data plots were made using ggplot2 v.3.4.3, and additional R packages including Seurat (Fig. 1D), EnhancedVolcano v.1.12.0 (Fig. 4D-E, Fig. 5F,H), and ComplexHeatmap v.2.10.0 (Fig. 3D). Alluvial plots displaying the frequency and phenotype of Pf-specific clonotypes over time were made using Phenoclone (Fig. 5A). Schematics (Fig. 1A, Fig. 3A) were made using BioRender.

Statistics

Comparisons of means were performed using paired, two-sided T tests of statistical significance (alpha = 0.05) (Fig. S7, Fig. S9, Fig. S10, Fig. S12). When multiple pairwise-comparisons were made, p-values were corrected using the Benjamini-Hochberg procedure. Differential gene expression was determined using a Wilcoxon Rank Sum test and p-values were adjusted using Bonferroni correction (Fig. 4D,E, Fig. 5F,H, Fig. S8D, Fig S14G,H, Fig. S16). Binomial regressions (Fig. 2F-H) and simple linear regression (Fig. 4H) were fit with 95% confidence intervals.