T cell dynamics and response of the microbiota after gene therapy to treat X-linked severe combined immunodeficiency

Human subjects

Patients were recruited as described [19, 33]. We collected the same sample types from six healthy children between the ages of 21–43 months under IRB 13-010072. We obtained sorted CD3+ T cells from three anonymous healthy adult donors above the age of 18 from the Human Immunology Core at the University of Pennsylvania. All samples were stored at − 80 °C.

Integration site analytical methods

Integration site sequences were determined using two different methods due to changes in technology over the period of patient monitoring (details in Additional file 1: Table S2). In the first, 454/Roche pyrosequencing was used to determine integration site placement [24]. In the second, Illumina paired end sequencing was used [21, 26, 27, 34]. In both, DNA was broken using shearing or restriction enzyme cleavage, then DNA adaptors ligated on to the broken DNA ends. Nested PCR was then used to amplify from the linker to the integrated vector, and the intervening segment of human DNA sequenced. All integration site sequence analysis was be carried out in quadruplicate to minimize PCR jackpotting. All sample sets were worked up together with human DNA lacking integrated lentiviral sequences to monitor for PCR contamination, which was typically undetectable. Different linkers were used for ligation-mediated PCR for each sample in a set to block PCR cross over. All samples were bar coded on both ends of the molecule, and only those with correct bar code pairs analyzed, thereby suppressing artefactual molecules resulting from PCR recombination. A total of 31 samples were analyzed, yielding a total of 24,170 integration sites.

TCR sequence analysis

TCR sequencing was performed on whole blood samples that had been fractionated to yield T cell (CD3+) or PBMC fractions. Genomic DNA was isolated and sequenced at Adaptive Biotechnologies using their immunoSeq protocol to determine CDR3 region sequences of the TCR-beta locus [35]. The protocol uses proprietary primer design and post-processing steps to implement controls, adjust for bias, and detect gene segments and sequence productivity. Data were analyzed using the immunoSeq Analyzer version 3.0. A total of 40 samples were analyzed, yielding a total of 32 million TCRB sequences. All statistics were performed only on the part of the clonal population for which a productive rearrangement was detected.

Microbiome sequencing

DNA was isolated from fecal, oral, and nare samples using the following procedures: Small aliquots of fecal material (< = 1 ml) and the tips of each swab (for oral and nares swabs) were deposited into a PowerSoil bead tube. DNA was extracted using the standard MoBio PowerSoil DNA extraction protocol, with one or more blank extraction controls worked up simultaneously with each set of samples. Work spaces were decontaminated using bleach and UV irradiation. The resulting DNA from all samples and blank controls was sequenced on an Illumina HiSeq 2500 using NextSeq chemistry and standard Illumina dual barcoding for each sample.

Due to artifacts in genomic sequences for members of the Apicomplexa family, many reads from a common water contaminant (Bradyrhizobium) were cross-annotated as belonging to Apicomplexa. Consequently, we removed all Apicomplexa reads before further analysis due to their uncertain provenance.

Virome analytical methods

Viral particles were isolated from a subset of the fecal samples using a protocol adapted from [36]. Fecal samples were homogenized and filtered through a 0.4-μm filter. The filtered samples were then treated with DNase and RNAseq to remove exogenous nucleic acids. Combined nucleic acids were then purified from the sample using the Qiagen UltraSens Virus Kit.

To obtain DNA viruses, we amplified viral genomes in an aliquot of the combined nucleic acids using the Illustra Genomiphi V2 DNA amplification kit. Resulting amplified DNA was quantified using PicoGreen and stored at − 20 °C. To obtain RNA viruses, we treated a separate aliquot of the combined nucleic acids with DNAse+ and then performed reverse transcription of the RNA to cDNA using the SuperScriptIII First-Strand Synthesis System from Life Technologies and second strand synthesis using Sequenase. The resulting cDNA was quantified with PicoGreen and stored at − 20 °C.

Resulting DNA was prepared and sequenced using NextSeq for library preparation and an Illumina HiSeq 2500 for sequence acquisition. The same post-processing pipeline was used, including all quality-control, host removal, and read annotation steps. For analysis, we considered only reads that fell under the virus classification and removed the following viral annotations as being reagent contamination [37, 38]: Enterobacteria phage M13, Enterobacteria phage T7, Enterobacteria phage phiX-174 sensu lato, Bacillus phage phi29, and Pseudomonas phage phi6, human herpesvirus 6 and 7, and Shamonda virus. Influenza sequences were detected, but only the HA segment was found. The HA sequence matched sequences from a previous run of the sequencing instrument, and so were excluded as carry over.

Bioinformatic methods

To estimate population sizes of T cell progenitors from integration sites, we used Chao1, and of TCR sequence data, we used Chao2.

Sequence reads for all metagenomic samples were processed using the Sunbeam pipeline (https://github.com/sunbeam-labs/sunbeam). Reads were quality-controlled by trimming low-quality bases and adapter sequences using Trimmomatic [39], and host reads were removed using BWA [40]. The remaining reads were filtered for low-complexity sequences using dustmasker [41] and Komplexity (https://github.com/eclarke/komplexity). The reads that remained were assigned taxonomy via Kraken [42] and a custom database built on all microbial genomic sequences in RefSeq release 79 [43].

Antibiotic resistance gene levels were assessed using ShortBRED [44] and a marker gene database built from CARD [45] available on the ShortBRED website (https://bitbucket.org/biobakery/shortbred/downloads/ShortBRED_CARD_2017_markers.faa.gz). The same reads used as input to the Kraken taxonomic classifier were used as input to ShortBRED.

T cell receptor antigen recognition was queried using VDJdb [46], a database and software package that compares CDR3 regions from TCR sequencing experiments to public epitope databases. The software compared our TCRB sequences to human TCRB sequences in the database with allowances for low numbers of indels and substitutions (default parameters), and only matches with high confidence (VDJdb score > 2) were kept.

Final analysis and figure generation were performed using the R statistical software [47]. Bray-Curtis dissimilarity and other ecological metrics were calculated using the R package vegan [48]. The full code listing and post-processed data used for analysis and figures are available online at https://github.com/eclarke/scid-multiomics-paper.

Experimental strategy

Our comparative analysis of gene-corrected progenitors and daughter T cells focused on nine patients from the SCIDn2 trial [19] and five patients from the SCIDn1 trial [8] for whom samples were available (Fig. 1a). Patient characteristics are summarized in Additional file 1: Table S1. Two of the subjects from the SCIDn1 trial developed leukemia requiring chemotherapy: F107 and F110 at 68 months and 33 months, respectively, after infusion of gene-modified progenitor cells [8, 19, 25]. The samples for these patients analyzed here are well after these adverse events, allowing assessment of the effects of leukemia and chemotherapy on progenitor cell and T cell populations. One subject in the SCIDn2 trial, U205, was transplanted twice—T cell numbers failed to reach a clinical protocol end point after the first infusion with modified cells, but succeeded the second time. Samples from the two treatment periods are designated U205 and U205b. Patient U204 was not successfully reconstituted but was rescued after allogeneic stem cell transplantation.

Vector integration sites and the T cell repertoire were monitored through regular per-protocol blood draws, followed by isolation of CD3+ populations. Sorted cells were subject to DNA extraction and then amplification of either integration sites (mostly previously reported in [8, 19, 25]; Additional file 1: Table S2) or mature TCR-beta loci (Additional file 1: Table S3; new data here), followed by sequencing. Vector integration sites mark progenitor cells capable of differentiating into mature peripheral T cells. Rearranged TCR-beta loci identify mature T cells present in the blood. Time points from the integration site analysis were chosen to match those used in the TCR analysis—characterization of additional time points for integration site distributions in SCID-X1 gene therapy was published previously [8, 18, 19, 33, 49, 50, 51].

Microbiome samples were available for six SCIDn2 patients (U201, U203, U204, U205, U207, and F201). Oral, nasal, and gut microbiota were sampled via collection of oropharyngeal swabs, nasopharyngeal swabs, and stool samples. Sampling times ranged from 4 to 181 months post infusion of corrected cells. Sample acquisition was at times limited by clinical and practical considerations.

As controls, we analyzed TCR-beta sequences from CD3+ cells from five healthy children and three healthy adults after informed consent. Healthy subject demographics are in Additional file 1: Table S1. Cross-sectional microbiome samples were collected from the same subjects and analyzed along with the SCID gene-corrected samples.

Integration site analysis

To characterize the hematopoietic progenitor cells capable of differentiating into peripheral blood T cells, we determined the sites of vector integration in patient chromosomes. Because of the large size of the human genome, each integration site uniquely marks a precursor cell capable of additional division after integration of the vector. Summaries of genes near integration sites in the most expanded clones at each time point are in Additional file 2, Figure S1. The numbers of integration sites detected in purified T cell samples ranged widely, from as few as 62 to as many as 2009. Reconstruction of population sizes using the Chao1 estimator suggested minimal sizes of 144 to 6018 active progenitors.

For the SCIDn2 subjects, we found that the estimated total number of vector integration sites was relatively constant over the time intervals analyzed (Fig. 1b), in the range of ~ 1000 predominant clones yielding circulating cells (Additional file 1: Table S2). Numbers varied by subject, with subject U203 showing consistently higher levels, while U205, a subject in which a suboptimal number of gene-modified cells were delivered, showing consistently lower levels.

For the SCIDn1 subjects, numbers of unique sites identified varied from 263 to 682; reconstructed minimal population sizes ranged from 628 to 1287. Comparison of mean values shows no difference in the numbers of unique integration sites in SCIDn2 subjects compared to SCIDn1 subjects (p = 0.17, Wilcoxon rank-sum test).

Population sizes of progenitors were compared between SCIDn1 subjects who suffered adverse events and were treated with chemotherapy (F107, F110) versus those who did not (F102, F106), or versus SCIDn2 subjects. Treating time points as independent tests, no systematic differences were detected comparing within SCIDn1 subjects; a modest difference in median values was detected with a comparison of adverse events versus no adverse events (adverse event samples compared to pooled SCIDn1 and 2; p = 0.041, Wilcoxon rank-sum test). Thus, the data suggest that adverse events and chemotherapy may have diminished the pool size, though we note that these subjects were also at relatively long times after gene therapy, so slow loss of diversity over time is a potential alternative explanation.

TCR-beta CDR3 analysis

Clonal structure of T cell populations was investigated by analyzing TCR-beta rearrangements in genomic DNA from peripheral blood CD3+ cells (Fig. 2). CDR3 region nucleotide sequences were conceptually translated into amino acid sequences, and numbers of productive rearrangements quantified (Fig. 2a). All the following analysis focused on T cell clones with productive rearrangements. For healthy adults, productive rearrangements ranged from 18,000 to 27,000 per sample. For healthy children, numbers ranged from 18,000 to 22,000 per sample. For SCIDn2, most samples were close to this range (12,000 to 33,000 per sample; no significant difference in medians). For SCIDn1, results were slightly lower, with four out of five samples between 15,000 and 23,000. The exceptional sample was from subject F107 who only showed 5,000 unique CDR3 sequences. F107 both suffered an adverse event at month 68 and was treated by chemotherapy, and demonstrated a low number of engrafted gene-modified cells at baseline. Patient F110 was also treated by chemotherapy following a serious adverse event at month 33, but was initially treated with a higher dose of gene-modified cells, potentially explaining the higher diversity in this patient.

For all the CDR3 samples studied, we could only sample a small fraction of the total CDR3 clonotypes in the subjects. To account for this, we used replicate sampling and the Chao2 population estimator to estimate a minimum bound for the true number of CDR3 clonotypes (Additional file 1: Table S3). Estimators can be sensitive to sampling effort, so our estimates represent minimum values. Focusing on samples analyzed with multiple replicates, we estimate population sizes of 26,000 to 2,600,000 CDR3 variants. The median for SCIDn2 samples was not different from that of healthy children, but the median for SCIDn1 was significantly lower than for the healthy children (p = 0.026; Wilcoxon rank-sum test), possibly due to the much longer follow-up time available for the SCIDn1 patients. Analysis of trends over time suggested that subjects U201, U203, and U207 repertoires became more rich over time (Fig. 2b). The one exception was the case of suboptimal reconstitution (U205), which did not show longitudinal increase.

To begin to characterize repertoire composition, we used Bray-Curtis dissimilarity and t-SNE to cluster samples (Fig. 2c). Replicates from each subject closely resembled those from the same subject at different time points, but differed from other subjects. No systematic differences were observed comparing SCIDn1 and SCIDn2, or comparing either data set to healthy controls (PERMANOVA test on Bray-Curtis dissimilarity, n.s.).

Recombination involving the most frequently used V (Fig. 2d) and J (Fig. 2e) gene segments was next quantified and compared. Usage of gene segments was quantified for healthy children and averaged, then the profile was compared to gene-corrected SCID subjects. The great majority of gene-corrected samples did not show significant differences from the distribution in healthy children, in the analysis of either V or J gene segment usage (chi-squared test). The only exception was V gene usage in early time points for SCIDn2 subject U201 (chi-squared test, p < 0.05). This subject was successfully corrected, so the unusual heavy use of specific V gene segments is not indicative of clinical failure. Of possible significance, U201 had an ulcer on his palate of unknown origin that became inflamed as T cells were produced, possibly skewing the repertoire. Fig. 2f shows the usage of V-J pairs within subjects, emphasizing the occasional outgrowth of expanded clones, potentially in response to antigen, followed by decrease of abundance.

TCR sequences were scanned for evidence of public epitopes in order to assess whether the subjects had formed typical responses to common antigens. However, no clear-cut link was found between putative epitopes and the microbial and viral taxa described below.

Tracking T cell ontogeny

Comparison of the estimated lower bounds for the number of unique integration sites per sample with the number of unique TCR sequences allows estimation of the minimum number of cell divisions required to generate the TCR-beta cell population from gene-corrected precursors (Fig. 3). We calculated the minimum number of cell divisions as the base-2 logarithm of the difference between the estimated population size of unique T cells and the estimated population size of integration sites:

CellDivisions=log2(TCRs−IntSites)CellDivisions=log2(TCRs−IntSites)

We found a relatively consistent range of minimum cell division values across patients and timepoints, with a median of 9.1 and a minimum and maximum of 5.6 and 15.4. The fraction of progenitor cells that die in the thymus is unknown, so our estimates are lower bounds for the required number of cell divisions.

Response of the microbiome to T cell reconstitution

Microbiota community structure was analyzed for six gene-corrected subjects (Fig. 4a) by extraction of DNA from swabs (oral and nares samples) or stool, followed by shotgun metagenomic sequencing (samples analyzed are in Table S4; summaries of clinical infections are in Table S5). Numbers of samples available per subject ranged from one to seven. The numbers of sequences acquired per sample averaged 72,000 (nasal; Fig. 4b), 8,746,000 (fecal; Fig. 4c), and 3,625,000 (oral; Fig. 4d). Sequencing reads were quality filtered as described in the “Methods” section and then assigned to microbial taxa using Kraken [42] .

The analysis of gut microbiota (Fig. 4c) emphasized the differences between gut microbiota of healthy children (right columns) versus SCID gene-corrected subjects (left columns). The healthy subjects were colonized predominantly with Bacteroides, which is typical of healthy gut, whereas the SCID subjects showed a range of major colonists. U201 was colonized mainly with Bifidobacteria, which is characteristic of healthy breast-fed babies [52]. However, the guts of subjects U203 and U205 were colonized with Veillonella at early timepoints. As Veillonella is commonly an oral bacteria, this is suggestive of abnormal colonization. Subject U204 showed high-level colonization with Enterobacteriaceae, typical of dysbiotic states. Early samples from subjects U207 and F201 showed abundant representation of Adenovirus and Bocavirus in the gut (Fig. 4c, bottom rows), comprising a significant portion of the classified taxa even though whole stool (and not virus-enriched material) was sequenced.

Oral samples (Fig. 4d) from the healthy controls were dominated by typical oral bacteria, including Prevotella, Streptococcus, Neisseria, and Haemophilus. In contrast, U204 and U205 showed high-level domination by Streptococcus, and U207 was dominated by Rothia. Viral sequences were rare in oral sample.

In nasopharyngeal samples (Fig. 4b), healthy subjects were dominated by Moraxella, Staphylococcus, and Propionibacterium. Bradyrhizobium was detected in B207 and healthy subjects, but this is a soil bacteria and a common contaminant, so we infer that detection of this lineage is artifactual. SCID subjects contained normal lineages to varying degrees, but also showed high-level colonization by Streptococcus, Corynebacterium, and others. In addition, nasopharyngeal samples from U203 were dominated by Bocavirus.

We compared the full microbial compositions of the samples using Bray-Curtis dissimilarities and found that sample types clustered by body site of origin (Fig. 5a, PERMANOVA p value < 0.001), as expected from many studies. Nasopharyngeal swabs also clustered near the reagent controls, suggesting that low authentic microbial biomass was present in these samples and that a significant proportion of the reads were likely derived from environmental contamination. We also found that each subject clustered with themselves throughout all timepoints, reflecting consistently higher inter-subject differences than intra-subject differences (PERMANOVA p value < 0.001). In several patients (U201, U205, and U207), the gut microbial communities began resembling healthy children more at late times after cell infusion (Fig. 5b), potentially reflecting a combination of improving immune function and reduction in usage of antibiotics and other medications. U203 and F201, whose microbial communities remained abnormal throughout the sampling period, were on prolonged antibiotics for treatment of disseminated Bacille Calmette-Guerin after gene therapy, potentially contributing to dysbiosis.

Analysis of taxonomic richness provides insight into microbial community health, since low richness is often associated with abnormal outgrowth of pathogenic organisms. For stool, richness increased for 4/6 subjects over time, potentially indicative of improving gut health (Fig. 5c). However, of the four subjects for whom oral samples were available, all showed abnormally low richness at every time point. This finding was associated with particularly high Streptococcus colonization in multiple subjects. Nasopharyngeal microbiota showed richness comparable to healthy controls.

SCID patients were treated with a wide variety of antibiotics before and after therapy to mitigate opportunistic infections. To assess whether this treatment was associated with an increase in antibiotic resistance gene representation, we quantified antibiotic resistance genes in the stool sequence samples using ShortBRED and the CARD antibiotic resistant factor database (Fig. 5d). We saw a greater quantity of antibiotic resistance genes in SCID subjects than in healthy controls in stool (pooled across timepoints, p = 0.048). Particularly high levels were seen in U204, for whom gene therapy was unsuccessful. Lesser frequencies of antibiotic resistance genes were evident in the oropharyngeal and nasopharyngeal samples. For patient U201, who had a timepoint 24 months after therapy, the level of antibiotic resistance genes decreased substantially between the 5- and 24-month timepoint. At this timepoint, the patient was no longer on any medications, which suggests that antibiotic resistance gene load may have decreased after the selective antibiotic was withdrawn. This suggests that the antibiotic administration, and potentially hospitalization, increased antibiotic resistance gene proportions.

Response of the virome

To enrich for virus detection beyond what was identified in whole stool sequencing (as in Fig. 4), viral particles were purified from stool, and DNA and RNA purified from particles and sequenced. Only patients with sufficient material were able to be virome-sequenced (patients U207 and F201 were excluded). Reads were assigned using Kraken, then extensive filtering was carried out to remove contaminants and artifacts. Figure 6 shows the remaining attributions. For RNA viruses, high-level colonization was detected with Astrovirus (U204 and U205) and Sapovirus (U203). For DNA viruses, numerous bacteriophage lineages were seen. Among DNA viruses infecting animal cells, high viral loads were seen for Adenovirus in U201, Bocavirus in U203, and Betatorquetenovirus (Anellovirus) in U205. The Adenovirus and Bocavirus infections diminished over time after successful gene correction. Betatorquetenovirus is a normal commensal in healthy individuals and persisted.