Early-life thymectomy leads to an increase of granzyme-producing γδ T cells in children with congenital heart disease

Abstract

Congenital heart disease (CHD) is the most common birth defect in newborns, often requiring cardiac surgery with concomitant thymectomy that is known to increase disease susceptibility later in life. Studies of γδ T cells, which are one of the dominant T cells in the early fetal human thymus, are rare. Here, we provide a comprehensive analysis of the γδ T cell compartment via flow cytometry and next-generation sequencing in children and infants with CHD, who underwent cardiac surgery shortly after birth. A perturbation of the γδ T cell repertoire is evident, and Vδ1 T cell numbers are reduced. However, those cells that are present, do retain cytotoxicity. In contrast, GZMA+CD28+CD161hi innate effector Vγ9Vδ2 T cells are found in higher proportions. TCR-seq identifies an increase in TRDJ3+ γδ T cell clones in children with CHD, but not in a confirmatory group of neonates prior to CHD surgery, which overall points to a persistence of fetal-derived effector γδ T cells in children with CHD.

Fig. 1. Early-life thymectomy affects the αβ and γδ T cell populations.

Flow cytometric and TCR-seq analyses of peripheral mononuclear blood cells from 5- to 12-year-old children with CHD (n = 18, CHD) after cardiac surgery involving total or partial thymectomy, and age-matched non-CHD controls (n = 11, ctrl). a Percentage of CD31⁺ cells within naive CD4 T cells, indicative of Recent Thymic Emigrants (RTE). b Dot plot shows CD3 T cell frequency of lymphocytes (defined as FSC-A^low/SSC-A^low). c Percentage of CD4 T cells (left) and CD8 T cells (right) among lymphocytes. d Stacked bar plot showing the percentage of a naïve (CD45RA⁺CCR7⁺), memory (CD45RA⁻, either CCR7⁺ (T_CM) or CCR7⁻ (T_EM)) or terminally differentiated effector memory phenotype (T_EMRA; CD45RA⁺CCR7⁻) among CD4 T cells. e Dot plot showing the percentage of PD-1 (left) or CD57 (right) expression among CD4 T cells and CD8 T cells, respectively. f TRB repertoire analysis of FACS-sorted CD4 or CD8 T cells. Inverse Simpson diversity indices of CD4 T cells (left) and CD8 T cells (right). g Dot plots present the percentage of CD3⁺ γδ T cells among lymphocytes. h Percentage of CD3⁺γδTCR⁺ Vδ1 or CD3⁺γδTCR⁺ Vγ9Vδ2 T cells within lymphocytes. i Stacked bar plots present the abundance of respective γδ T cell subsets as median values within ctrl and children with CHD. Statistical analyses were performed using the two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). No adjustments for multiple comparisons were made. Horizontal bars indicate median values. Each dot represents one donor. Source data are provided as a Source Data file.

Fig. 2. Enrichment of TRDV2⁺TRDJ3⁺ γδ T cell clones in children with CHD.

TRD repertoires analyses of FACS-sorted peripheral blood γδ T cells from 5- to 12-year-old children with CHD post-surgery (n = 16, CHD) and age-matched controls (n = 11, ctrl). a Inverse Simpson indices of TRDV1⁺ (left) and TRDV2⁺ (right) TRD repertoires. b Distribution of paired TRDV and TRDJ gene elements calculated from all ctrl and all CHD samples are depicted in the pie chart. The size of pie segments represents the median value of respective gene pairs calculated from all ctrl or CHD samples. c J-gene usage in TRDV1⁺ (left) and TRDV2⁺ (right) TRD sequences in percent, differentiated according to respective J-gene (TRDJ1, TRDJ2, TRDJ3, TRDJ4). d Percentage of RTE (CD31⁺ cells of naïve CD4⁺ T cells) plotted against the proportion of TRDJ3⁺TRDV2⁺ sequences for children with CHD, showing negative correlation (Spearman’s rho = −0.69, p-value = 0.004), and the age-matched controls (Spearman’s rho = 0.36, p-value = 0.33). The error band describes the standard deviation of the displayed data points. e N-insertion counts in TRDV1⁺ (left) and TRDV2⁺ (right) sequences, calculated as weighted mean, weighted by the total number of clones per sample. Statistical analyses were performed using the two-sided Wilcoxon–Mann–Whitney U test. No adjustments for multiple comparisons were made; not significant (ns). Bars indicate median values; each dot represents one donor. Source data are provided as a Source Data file.

Fig. 3. Increase of CD28^hiCD161^hi Vγ9Vδ2 T cells in children with CHD.

High-dimensional flow cytometry data analysis of children with thymectomy (CHD, n = 18) and non-CHD controls (ctrl, n = 11) reveals an increase CD28^hiCD161^hi Vγ9Vδ2 T cells in 5- to 12-year-old CHD children after early-life cardiac surgery. a UMAP of obtained clusters (c1–c8) from pre-gated live CD3⁺TCRγδ⁺ cells, clustered by an unsupervised clustering approach based on surface marker expressions, and respective cluster proportions per group presented as stacked bar plot. b Expression UMAPs of obtained clusters colored by Vδ2 and Vδ1 surface marker expression. c Ridge plots show the expression of selected markers (CD161, NKG2A, CD31, CD8, CD57 and CD28) in each cluster. The vertical lines indicate the median marker expression within the respective cluster. d Heatmap visualizing the median scaled expression of surface markers within each cluster c1–c8. e Cluster frequencies among γδ T cells, compared between pediatric CHD patients and non-CHD controls. The clusters c1–c8 were annotated based on their surface marker expression characteristics. f Expression UMAPs of obtained clusters colored by surface marker expression of cluster defining markers. g Manual gating of Vγ9Vδ2⁺CD28⁺CD161⁺ and Vδ1⁺CD31⁺ γδ T cell populations in representative samples. h Percentage of Vγ9Vδ2⁺CD28⁺CD161⁺ and Vδ1⁺CD31⁺ cells among all γδ T cells. Statistical analyses were performed using the two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). No adjustments for multiple comparisons were made. Bars indicate median values, with each dot representing a single donor. Source data are provided as a Source Data file.

Fig. 4. Elevated expression of cytotoxicity gene module scores in thymectomized children with CHD.

a UMAP visualization of identified clusters from scRNA-seq of FACS-sorted γδ T cells from 5- to 12-year-old children who received thymectomy early after birth (CHD, n = 4) and control children (ctrl, n = 4), colored by cluster. b The expression of TRDV1 and TRDV2 transcripts per each cluster. c The bar plot displays fractions (%) of absolute cell numbers from ctrl and CHD γδ T cells that contributed to clusters c1–c6. d Dot plot of the average gene expression (columns) per cluster (rows). Dots are colored by average logFC (Scaled Exp) and sized by the percentage of cells per cluster that expressed this gene (% Exp). e–g Box plots of the single-cell gene signature module score for recent thymic emigrants (RTE) (e), TCF7 (f), and cytotoxicity (g) of each cluster in both conditions. The RTE module score was computed based on KLF2, CCR9, PECAM1, S1PR1, LEF1, TCF7, SOX4, NT5E, and SELL. The cytotoxicity score was computed based on GZMA, GZMB, GZMH, GZMK, NKG7, GNLY, KLRK1, KLRD1, KLRF1, KLRC1, KLRB1, KLRG1, CCL5, CCL4, and CXCR6. h Left: protein expression of GZMA and GZMB in one representative ctrl and one CHD child on alive CD3⁺γδTCR⁺ Vδ1 T cells. Right: the frequency of GZMA^posGZMB^pos cells among Vδ1 T cells and total γδ T cells shown as mean ± SEM. n_ctrl = 3, n_CHD = 9. e–g The cell numbers of each cluster from 4 ctrls and 4 CHDs are as follows: c1, nctrl = 1653, nCHD = 91; c2, nctrl = 405, nCHD = 32; c3, nctrl = 389, nCHD = 42; c4, nctrl = 1989, nCHD = 696; c5, nctrl = 1851, nCHD = 1511; c6, nctrl = 77, nCHD = 5816. The minima of boxplot: the lowest data point within 1.5 times the first quartile’s Interquartile Range (IQR). Maxima: the highest data point within 1.5 times the IQR of the third quartile. Center: median of the data. e–h Statistical significance was determined by unpaired two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). Source data are provided as a Source Data file.

Fig. 5. Increase of GZMA⁺ innate effector Vδ2 T cells in children with CHD.

a Volcano plot shows the average expression of the differentially expressed genes between Vδ2 cluster c6 and other Vδ2 clusters (c2-c5). b Box plots of the single-cell gene signature module score for KLRB1 and CD28. The cell numbers of each cluster from 4 ctrls and 4 CHDs are as follows: c1, nctrl = 1653, nCHD = 91; c2, nctrl = 405, nCHD = 32; c3, nctrl = 389, nCHD = 42; c4, nctrl = 1989, nCHD = 696; c5, nctrl = 1851, nCHD = 1511; c6, nctrl = 77, nCHD = 5816. The minima of boxplot: the lowest data point within 1.5 times the first quartile’s Interquartile Range (IQR). Maxima: the highest data point within 1.5 times the IQR of the third quartile. Center: median of the data. c Flow cytometric analyses depict GZMA, GZMB, CD28 and CD8A expression of unstimulated alive CD3⁺γδTCR⁺Vδ2⁺ T cells in one representative control and one child with CHD. d, e The frequency of CD8A^posCD28^posGZMA^posGZMB^neg (d) and GZMB^pos (e) among Vδ2 T cells, shown as mean ± SEM. n_ctrl = 6, n_CHD = 9. f Representative FACS plot of CD28 and CD8B expression on Vδ2 T cells in one control (1 out of 1) and one child with CHD (1 out of 4). g, h In vitro HMBPP + IL-2 stimulation of peripheral blood samples of 5- to 12-year-old CHD children (n = 7) and controls (n = 5), or IL-2 alone was performed. g Intracellular GZMB expression in both unstimulated (only IL-2) and stimulated (IL-2 plus HMBPP) PBMCs of alive/CD3⁺/γδTCR⁺/Vγ9Vδ2 T cells in representative samples, and as determined by frequency at day 7 post-stimulation. h TNFα-producing Vγ9Vδ2 T cells in one representative ctrl and CHD sample, and determined by frequency at 7 days post-stimulation. c–g Statistical analyses were performed using by unpaired two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). Bars indicate median values; each dot is representative of one donor. The error band describes the standard deviation of the displayed data points. Source data are provided as a Source Data file.

Fig. 6. γδTCR repertoires in neonates with CHD but prior cardiac surgery are highly similar to age-matched neonatal controls.

TRG and TRD repertoire analyses from neonates with CHD prior cardiac surgery (nCHD, n = 15) and neonatal controls (neonatal ctrl, n = 16). a Inverse Simpson indices of TRGV9⁺ TRG sequences. b Proportion of TRGV9⁺ clones within all TRG sequences. c Inverse Simpson indices of TRDV2⁺ TRD clones. d Pie charts, color-coded by TRDV-TRDJ pairing within TRD repertoires at group-level. e J-gene usage in TRDV2⁺ sequences in percent. f N-insertion counts in TRDV2⁺ sequences, calculated as weighted mean, weighted by the total number of clones per sample. Statistical analyses were performed using the two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). No adjustments for multiple comparisons were made. Bars indicate median values; each dot is representative of one donor.

Fig. 7. Longitudinal single-cell transcriptomic analysis reveals a proliferative advantage of γδ T cells in infants with CHD at 6 months post-thymectomy.

a UMAP visualization of the γδ T cell scRNA-seq dataset from two longitudinally followed CHD infants before (Baseline, BSL, at age day 8 and 9) and 6 months post-surgery (Follow-up, FU, at age day 188 and 189), colored by the identified clusters c1–c8. Control samples represent transcriptional profiles of age-matched preterm neonates. b Dot plot of the average gene expression (columns) per cluster (rows). Dots are colored by average logFC (Scaled Exp) and sized by the percentage of cells per cluster that expressed this gene (% Exp). c The bar plot reveals fractions of absolute cell numbers of indicated subpopulations from each group. d The volcano plot demonstrates DEGs among the Ctrl FU c6–c7 and CHD FU c6–c7. Upregulated DEGs are identified with log2FC > 0.25 and p_val_adj <0.05. e Gene ontology enrichment analysis of the upregulated differential expressed genes from the indicated subpopulations. f Box plots of the single-cell gene signature module score for proliferation, based on gene list: CBLB, PTPN6, IRF1, HLA-DPB1, HLA-DPA1, HLA-DRB1, CCL5, TNFRSF1B, SH2D2A, FYN, BTN3A1, IL12RB1, SPN, CD81, MSN, ANXA1, SASH3, PYCARD, PTPN22, RASAL3, SELENOK, TNFSF14, ITCH, HLA-A, PSMB10, PTPRC, HLA-E, TNFRSF14, IL2RA, CD40LG, TNFRSF4, CR1, ANXA1, XCL1, SOS1, IL6ST, CD55, PRNP, and PPP3CA. The cell numbers of each cluster from each group are as follows: c1, nCtrl BSL = 235, nCtrl FU = 384, nCHD BSL = 215, nCHD FU = 55; c2, nCtrl BSL = 399, nCtrl FU = 2296, nCHD BSL = 1432, nCHD FU = 129; c3, nCtrl BSL = 300, nCtrl FU = 1320, nCHD BSL = 754, nCHD FU = 138; c4, nCtrl BSL = 788, nCtrl FU = 282, nCHD BSL = 229, nCHD FU = 138; c5, nCtrl BSL = 89, nCtrl FU = 31, nCHD BSL = 159, nCHD FU = 127; c6, nCtrl BSL = 2825, nCtrl FU = 1589, nCHD BSL = 141, nCHD FU = 1401; c7, nCtrl BSL = 662, nCtrl FU = 671, nCHD BSL = 27, nCHD FU = 923; c8, nCtrl BSL = 58, nCtrl FU = 56, nCHD BSL = 5, nCHD FU = 3. The minima of boxplot: the lowest data point within 1.5 times the first quartile’s Interquartile Range (IQR). Maxima: the highest data point within 1.5 times the IQR of the third quartile. Center: median of the data. Statistical significance was determined by unpaired two-sided Wilcoxon–Mann–Whitney U test; not significant (ns). Source data are provided as a Source Data file.