Patients

The patients included in this study were followed at the Hospital Sant Joan de Déu and Hospital de la Santa Creu i Sant Pau (Barcelona, Catalonia, Spain). Written informed consent for the collected samples and analyses reported here, and for publication of the article, was obtained according to the procedures of the local Ethics Committee [code: PR(AG)202/2021]. This study was conducted in full compliance with the Declaration of Helsinki and all relevant ethical regulations for human subjects research.

Sanger sequencing

Genomic DNA was obtained from different samples: peripheral blood (collected in EDTA tubes), saliva and purified cell populations (CD3⁺, CD4⁺ T cells, CD8⁺ T cells, CD19⁺ B cells and CD15⁺ polymorphonuclear cells). DNA was isolated using the QIAsymphony SP purification system (Qiagen).

Polymerase chain reaction (PCR) was performed to amplify the regions of interest (primers and PCR conditions are available upon request), and purified PCR products were sequenced on an ABI 3500 DNA Sequencer using the BigDye Terminator sequencing kit 3.1 (Applied Biosystems).

Whole exome sequencing (WES)

Genomic DNA from the patient, obtained prior to HSCT, was used for WES. SeqCap EZ MedExome Target Enrichment Kit (Roche Nimblegen) was used to perform whole-exome enrichment in combination with SeqCap EZ Mitochondrial Genome Design (Roche Nimblegen). The short-insert paired-end libraries for the whole exome capture were prepared with KAPA HyperPrep kit (Roche). In short, 1 microgram of genomic DNA was sheared on a Covaris™ LE220-Plus (Covaris). The DNA fragments were end-repaired, adenylated, and xGen Dual Index UMI Adapters (Integrated DNA Technologies) were ligated. The adaptor-modified end library was enriched by 7 pre-capture PCR and with 14 cycles of post-capture PCR. The size, concentration and quality of the captured material were determined using an Agilent Bioanalyzer DNA 7500 chip.

The WES libraries were sequenced on HiSeq4000 (Illumina) in paired-end mode with a read length of 2 x 101 bp following the manufacturer’s protocol for dual indexing. Image analysis, base calling and quality scoring of the run were processed using the manufacturer’s software Real-Time Analysis (v2.7.7) and followed by the generation of FASTQ sequence files. The data was mapped and aligned against the human reference genome (GRCg38/hg38) using BWA¹⁰. Following the GATK best practices¹¹, we performed base quality score recalibration and local indel realignment before variant calling. The resulting VCF file was annotated using ANNOVAR¹².

NGS-based deep amplicon sequencing (NGS-DAS)

Amplimers containing the candidate region were PCR-generated using 40 ng of DNA. PCR products were quantified with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific). A 300- to 900-ng amount of DNA was fragmented using NEBNext dsDNA Fragmentase (New England Biolabs) to obtain fragments of approximately 200 bp in size. End repair, ligation of the adapter for Illumina, and sample indexing by a small amplification were then carried out using the NEBNext Ultra DNA library prep kit for Illumina and the NEBNext MultiplexOligos for Illumina Dual Index Primers Set 1 (New England Biolabs). All necessary purifications and size selections to recover only the fragments of interest were done using AMPure XP beads (Beckman Coulter).

Libraries were quality-evaluated using QIAxcel (Qiagen) and quantified with a Qubit 2.0 Fluorometer. The DNA sample libraries were then mixed in equimolecular amounts and sequenced in a MiSeq instrument (Illumina) using the 500-cycle MiSeq reagent kit v2 with a paired-end run of 2 × 250 bp reads. All procedures were performed according to the manufacturer’s instructions.

For the bioinformatic analysis, NGS sequencing data were aligned to the reference genome, hg38, using BWA (https://github.com/lh3/bwa). Duplicate reads were marked using Picard MarkDuplicates (version 2.18.6, https://github.com/broadinstitute/picard). Next, GATK 4.1.8.1 (https://github.com/broadinstitute/gatk) was used to calibrate the data, calculate BQSR scores (BaseRecalibrator, ApplyBQSR), and perform variant calling (Mutect2). The resulting VCF files were annotated using Annovar (https://annovar.openbioinformatics.org/en/latest/). Finally, Integrative Genomic Viewer was used to count the number of reads supporting each allele.

Single cell RNA sequencing (scRNAseq): sample preparation and sequencing

Sorted CD8+ cells were centrifuged at 400 x g for 5 min at 4 °C, and final concentration and viability were assessed using a TC20™ Automated Cell Counter (Bio Rad) upon staining with Trypan Blue. The sample was encapsulated for a Target Cell Recovery of 6000 cells on the Chromium Controller instrument (10X Genomics), using standard throughput Chromium Next GEM Single Cell 5’ Reagent Kit v2 (10X Genomics, PN-1000263). Libraries were prepared following the manufacturer’s instructions of protocol CG000331. Briefly, after GEM-RT cleanup, cDNA was amplified during 13 cycles, purified, and quantified on an Agilent Bioanalyzer High Sensitivity chip (Agilent Technologies). To construct the Gene Expression (GEX) library, 10 ng of cDNA were fragmented, end-repaired, A-tailed, and sample-indexed using the Chromium Single Cell 5’ Library Construction Kit (10X Genomics PN-1000190) and the Dual Index Kit TT Set A (10X Genomics PN-1000215). Human T cell receptor (TCR) sequences were enriched from the amplified full-length cDNA with the Chromium Single Cell Human TCR Amplification Kit (PN-1000252). Fragmentation, end repair, A-tailing and library indexing of the enriched cDNA were carried out using the aforementioned kits. Finally, the size distribution and concentration of 5’ GEX and TCR libraries were verified on an Agilent Bioanalyzer High Sensitivity chip. Sequencing was carried out on a NovaSeq6000 system (Illumina), aiming for approximately 40,000 and 10,000 read pairs per cell for the GEX and the TCR libraries, respectively. The following sequencing conditions were used: 28 bp (Read 1) + 10 bp (i7 index) + 10 bp (i5 index) + 90 bp (Read 2).

Single cell RNA sequencing (scRNAseq): data analysis

The 10x Genomics Cell Ranger software was used for demultiplexing, alignment to GRCh38 human reference genome, filtering, and generation of feature-barcode matrices through unique molecular identifiers (UMIs). Seurat v5.1.0¹³, and R v4.4.1 were used for downstream analysis and visualization. Cells containing fewer than 600 features and/or with more than 7% mitochondrial content were considered of low quality and subsequently removed. Estimated doublets were removed using DoubletFinder v2.0.4¹⁴. The raw count data were normalized with SCTransform, regressing for mitochondrial content and cell cycle phase scoring¹⁵. Principal Component Analysis was applied to the SCTransform-normalized data to identify significant principal components. These components were then used for constructing a Shared Nearest Neighbor graph. Louvain clustering was applied to this graph to determine cell clusters. Uniform Manifold Approximation and Projection was utilized for the visualization of clusters in a two-dimensional space. Automatic annotation of cells was performed using SingleR v2.6.0¹⁶, and cell types different from CD8 + T cells were removed. After filtering, SCTransform transformation was applied a second time.

To distinguish between cells from the patient and her brother, cells were classified based on specific gene expression. TNFRSF9 expression was observed in 269 CD8+ T cells. Since scRNAseq does not typically sequence the entire transcript, the somatic reversion could be detected only in a few of those cells, but all cells carrying R1 (8 cells) and R2 (3 cells) variants were classified as male cells (XY) based on the expression of genes located outside the pseudoautosomal region of the Y chromosome. Cells expressing genes located outside the Pseudoautosomal Region (PAR) of the Y chromosome (such as USP9Y, EIF1AY, TTTY10, TTTY14, KDM5D, RPS4Y1, ZFY, DDX3Y, UTY, TMSB4Y, NLGN4Y, LINC00278, RPS4Y2, and PCDH11Y) were classified as the brother’s cells. Conversely, cells expressing the gene XIST were classified as the patient’s cells.

The combineTCR function from scRepertoire v2.0.0¹⁷, was used to concatenate clonotype information into an integrated Seurat object. The repertoire landscape was visualized using the built-in functions of scRepertoire.

Cells expressing either one of the genetic variants of interest in TNFRSF9 were selected using samtools v1.20. To infer the presence of the three specific variants in cells where the target gene (TNFRSF9) could not be read directly, a TCR-based proximity approach was employed. Leveraging the clonal nature of TCRs, cells exhibiting the same TCRs—yet lacking direct gene readout—were predicted to also harbor the mutations of interest. This assumption was based on the rationale that cells with identical TCRs have undergone clonal expansion from a common progenitor cell and, therefore, are likely to share genomic mutations present in that progenitor.