Research ArticleIMMUNOLOGY

Weighting tumor-specific TCR repertoires as a classifier to stratify the immunotherapy delivery in non–small cell lung cancers

See allHide authors and affiliations

Science Advances  19 May 2021:
Vol. 7, no. 21, eabd6971
DOI: 10.1126/sciadv.abd6971


Analysis of T cell receptor (TCR) repertoires may contribute to better understanding of the response to immunotherapy. By deep sequencing of the TCR β chain complementarity-determining regions in the paired biopsies and peripheral blood specimens of 31 patients with non–small cell lung cancer (NSCLC) treated with anti–programmed death 1 (PD-1) or PD-ligand 1 (PD-L1) therapy, we developed a previously unidentified index, the TCR-based immunotherapy response index (TIR index), that estimated the degree of overlap of the TCR repertoire between tumor-infiltrating lymphocytes and circulating PD-1+CD8+T cells (shared TCR clones). This index correlated with response and survival outcomes of anti–PD-(L)1 treatment. All the TCR sequences of neoantigen-stimulated T cells were included in the shared TCR clones, indicating that TCR clones involved in TIR index estimation represented tumor-specific T cells. Therefore, the TIR index is a feasible approach for assessing tumor-specific TCR and stratifying patients with NSCLC for anti–PD-(L)1 therapy.


Immune checkpoint inhibitors (ICIs) targeting programmed death 1 (PD-1) or PD-ligand 1 (PD-L1) have been extensively administered and have resulted in substantial survival benefits in patients with non–small cell lung cancers (NSCLCs) (15). However, patients who do not benefit from this treatment remain a major challenge despite the recent identification of several potential, yet still clinically inadequate, biomarkers, including PD-L1 expression (6), tumor mutational burden (TMB) (7), and gene expression profile (GEP) score (8), highlighting the urgent need to explore novel predictors.

Elucidation of the contributions of various immune cells in the tumor microenvironment (TME) to antitumor immunity is likely to be a fundamental step in better matching these therapies to specific patients likely to benefit. One central feature driving the antitumoral response is the activation of tumor neoantigen–specific tumor-infiltrating lymphocytes (TILs), more specifically, PD-1–expressing CD8+ T cells (PD-1+CD8+TILs), the main subgroup previously identified (9). A superior response to anti–PD-(L)1 therapy may be indicated by a high density of PD-1+CD8+TILs, which can be roughly evaluated by semiquantitative approaches such as immunohistochemical staining or flow cytometry. The efficacy of ICIs could also be better predicted by the identification of specific neoepitopes that can effectively stimulate TILs in vivo (10, 11). However, the complexity of this process limits the clinical application of these techniques. As a pivotal surrogate for TIL expansion upon recognition of neoantigens, the T cell receptor (TCR) repertoire allows direct and precise assessment of the T cell immune activation, representing a quantitative estimation of neoantigen recognition. However, most targeted tumor antigens remain unknown and are highly variable between patients. The diversity and clonality of the TCR repertoire have previously been reported to be associated with the response to ICIs, but these results remain controversial (12, 13). Our previous study demonstrated that the baseline TCR diversity of PD-1+CD8+TILs but not peripheral blood mononuclear cells (PBMCs) could predict the response to ICIs (14), suggesting the importance of isolating tumor-infiltrating T cells for assessment of the TCR repertoire. Theoretically, a high proportion of PD-1+CD8+T cells in the overall TILs should reflect an increased possibility of the enrichment of functional antigen-specific TILs, but how to accurately and conveniently measure this enrichment through TCR repertoire assessment for practical clinical usage is unclear.

The neoantigen-specific TCR repertoires of the circulating and tumor-infiltrating PD-1+ CD8+ T cells have also been reported to be similar (15), implying that expanded intratumoral TCR clones can be identified in matched blood samples, and thus, circulating PD-1+CD8+TCR repertoires could be used as substitutes for those in the tumor. In this context, we aimed to elucidate the relationship between the intratumoral TILs and the circulating PD1+CD8+ T cell subpopulations via detailed genomic profiling and TCR repertoire analyses, followed by functional assays of autologous T cell activation in paired samples from ICI-treated patients with advanced NSCLCs. Subsequently, a previously unidentified indicator of the survival outcomes of anti–PD-(L)1 therapy as well as the corresponding TME features was developed, and its potential clinical applications were explored.


Exploration and clinical validation of a previously unidentified TCR-based immunotherapy response index

The aggregation of “bystander” (non–tumor-specific) TILs activated by ICIs in the TME might result in a high TIL that may inadequately distinguish patients who will benefit from ICIs. PD-1+CD8+T cells represent a subset of T cells containing reinvigorated tumor-specific TILs (16), which could compensate for the deficiencies of TILs alone as a predictive factor. Pre-ICI biopsies and matched peripheral blood specimens from 31 consecutive ICI-treated patients with NSCLC were assessed by next-generation deep sequencing (NGS) for TCR β chain complementarity-determining regions (CDR3s) (table S1). The overall study algorithm is presented in Fig. 1A. A median of 247 TCR repertoire sequences (range, 54 to 805) were identified in all individual tumor tissues. There was no significant difference in the TCR clone numbers between patients who benefited [confirmed complete response (CR) + partial response (PR) + stable disease (SD)] and those with no benefit (progressive disease at first tumor evaluation) from immunotherapy. The diversity and clonality of the TCR repertoire were evaluated by the Shannon metric (see Materials and Methods). Neither TCR diversity nor clonality could differentiate the patients who benefited from the patients who did not benefit from immunotherapy (table S2). Besides, no significant correlation was found between the ratio of PD-1+CD8+ T cells to total PBMCs and response to immunotherapy (P = 0.195; fig. S1).

Fig. 1 Correlations of TIR index with clinical outcomes upon anti–PD-1/PD-L1 treatment.

(A) Overview of study design. (B) Definition of the TIR index. TIR index was calculated with the overlapping degree of TCR repertoire between the tumor tissue [part 1 in (A)] and the isolated peripheral PD-1+CD8+T cells [part 2 in (A)]. (C) Comparison of TIR index distributions between nonresponders (red points, n = 25) and responders (gray points, n = 6) to anti–PD-(L)1 treatment. (D) Identification of cutoff value of TIR index. An ROC curve was used to distinguish responders from nonresponders. (E) Association of TIR index with PFS. (F) Association of TIR index with OS.

We speculated that the bystander-expanded TCRs may have diluted the tumor antigen–specific TCR repertoires, which were previously reported to be similar between the circulating and tumor-infiltrating PD-1+CD8+T cells (15). We therefore calculated the degree of overlap of the TCR repertoire between the tumor tissues and the isolated peripheral PD-1+CD8+T cells (called shared TCR clones; Fig. 1B and table S3). We hypothesized that the degree of shared TCR clones in general TIL TCRs [namely, the TCR-based immunotherapy response index (TIR index)] before therapy indicated a high proportion of the tumor antigen–specific T cells preexisting in the TME and would indicate an improved response to immunotherapy (Fig. 1B). Consistent with this hypothesis, a positive correlation was observed between the TIR index and the response to anti–PD-(L)1 treatment {responder versus nonresponder; median, 2.34% [95% confidence interval (CI), 1.78 to 4.09] versus 0.64% (95% CI, 0.57 to 1.17), P < 0.001; Fig. 1C}. The cutoff value is obtained by maximizing the Youden index. On the basis of an optimal threshold of 1.91 with the maximum receiver operating characteristic (ROC) curve value (area under curve, AUC = 0.974), the TIR index was used to dichotomize patients into “high” and “low” TIR index subsets (Fig. 1D) and exhibited a Youden’s index of 0.96 (100% sensitivity, 95% CI, 54.07 to 100% and 96.2% specificity, 95% CI, 80.36 to 99.9%) for stratifying the clinical response to ICIs (responder versus nonresponder). Patients with a high TIR index demonstrated a longer progression-free survival (PFS) and overall survival (OS) than those with a low TIR index [median PFS, 10.7 versus 2.3 months; hazard ratio (HR), 0.40; 95% CI, 0.19 to 0.85, log-rank P = 0.022; Fig. 1E; median OS, not reached versus 6.9 months; HR, 0.21; 95% CI, 0.08 to 0.53, log-rank P = 0.018; Fig. 1F]. These results indicated that the overlapping degree of TCRs between the TILs and the peripheral PD-1+CD8+T cells correlated well with the efficacy of anti–PD-(L)1 treatment.

Combined use of the TIR index with other biomarkers to predict response to immunotherapy

To further evaluate the predictive performance of the TIR index for ICI therapy, we measured the PD-L1 expression on tumor cells and the tissue-based TMB (tTMB) in patients with NSCLC receiving ICIs by using multiplex immunofluorescence (mIF) (n = 31) and targeted NGS with a 1021-gene panel (n = 26) (table S4), respectively. A high correlation for the TMB estimation between targeted NGS and whole-exome sequencing based on The Cancer Genome Atlas (TCGA) database was obtained (r2 = 0.973; fig. S2). In addition, PD-L1 expression on tumor cells (1% as the cutoff) (2, 17) was positively associated with the efficacy of anti–PD-(L)1 treatment (responder versus nonresponder, P = 0.043; Fig. 2A), which was not observed for the tTMB (14 as the cutoff; response versus nonresponse, P = 0.76; Fig. 2B). Although all the patients with a high TIR index (n = 7) had positive PD-L1 expression (≥1%), the distribution of tTMB and PD-L1 expression did not correlate with TIR index (Fig. 2C). These observations suggested that the TIR index might be a biomarker independent of PD-L1 expression and tTMB. Moreover, by using multivariate Cox regression models including the TIR index, PD-L1 expression, tTMB, pathological type, gender, and age, we identified the TIR index as an independent factor for both PFS (HR, 0.47; 95% CI, 0.28 to 0.89, P = 0.013; Fig. 2D, top) and OS (HR, 0.62; 95% CI, 0.33 to 0.91, P = 0.030; Fig. 2D, bottom). When the TIR index was combined with PD-L1 expression or tTMB, the objective response rates (ORRs; Fig. 2E) and duration of response of ≥6 months (Fig. 2F) were further improved in patients with a high TIR index/PD-L1 positivity (≥1%) or a high TIR index/high tTMB, respectively, which is superior to the predictive power of the individual assays, although no significantly statistical difference was observed between sole TIR index and combined models (including TIR index, PD-L1 expression, and tTMB) (P > 0.05; fig. S3).

Fig. 2 Utility of TIR index alone or in combination with previously established predicting indexes.

(A) Comparison of the distribution of PD-L1 expression between responders (gray points, n = 6) and nonresponders (red points, n = 25) to anti–PD-(L)1 treatment. (B) Comparison of tTMB distribution between responders (red points, n = 6) and nonresponders (gray points, n = 25) to anti–PD-(L)1 treatment. (C) Distribution of TMB and PD-L1 expression in different TIR index groups. (D) Multivariate Cox regression model assessing the correlation of clinical variables with PFS (top) and OS (bottom). (E) ORRs in different subgroups divided by individual or combined predicting indexes to anti–PD-(L)1 treatment. (F) Proportions of patients with PFS ≥ 6 months in different subgroups divided by individual or combined predicting indexes to anti–PD-(L)1 treatment. (G) Pie chart showing the distribution of EGFR mutation subtypes in 10 patients harboring EGFR mutations. (H) ORRs obtained according to different predicting indexes to anti–PD-(L)1 treatment in patients harboring EGFR mutation (n = 10). (I) Association of TIR index with PFS in EGFR mutated patients. (J) Comparison of TIR index distribution between patients with (mutant type, gray points, n = 19) and without (WT, red points, n = 12) oncogenic driver mutations (including EGFR exon 19 deletion/L858R/T790M, KRAS G12C/G12D, and MET exon 14 skipping) or potential resistant gene aberrations (including DNMT3A V560M and JAK 1/2 loss) to anti–PD-(L)1 therapy.

The TIR index as a classifier for immunotherapy response in a specific patient population with gene alterations potentially resistant to anti–PD-(L)1 treatment

EGFR mutant NSCLCs are considered representative tumors with a “cold” immune microenvironment, and the efficacy of ICIs is highly limited (18). Several gene alterations, such as MDM2/4 amplification, DNMT3A mutation (19), and JAK 1/2 mutation (20), have also been shown to be associated with resistance to anti–PD-(L)1 therapy. To validate the predictive value of the TIR index, we first analyzed the correlation of this index with the efficacy of anti–PD-(L)1 treatment in 10 epidermal growth factor receptor (EGFR) mutant patients (Fig. 2G), and all patients with a high TIR index obtained an objective response (Fig. 2H). The patients with a high TIR index demonstrated prolonged PFS compared with those with a low TIR index (median, 10.7 versus 3.1; Fig. 2I). Similarly, in 19 patients harboring driver gene mutations (including EGFR exon 19 deletion/L858R/T790M, KRAS G12C/G12D, and MET exon 14 skipping) and those potentially associated with disease progression (including DNMT3A V560M and JAK 1/2 loss) (fig. S4A), a high TIR index still yielded superior ORRs (75%; fig. S4B) and PFS (HR, 0.30; 95% CI, 0.12 to 0.74, P = 0.006; fig. S4C) with ICIs in this population, which was not obtained when stratified by PD-L1 expression (1 or 25% as the cutoff; fig. S4, D and E) or tTMB (≥14 as the cutoff; fig. S4F). Furthermore, no correlation was observed between the distribution of the TIR index and the aforementioned gene alterations that were associated with poor efficacy of anti–PD-(L)1 treatment (Fig. 2J).

In vitro validation of the TIR index as an indicator of functional tumor-specific cytotoxic T cells

We next sought to verify whether the shared TCR clones used in the TIR index calculations represented functional mutation-associated neoantigen (MANA)–specific T cells. Peptides were designed and synthesized according to the candidate MANAs and their corresponding wild-type (WT) sequences predicted via a neoantigen prediction algorithm in four patients with PR and two patients with progressive disease (range, 1 to 6 MANAs per patient) (table S5). A majority of the predicted MANAs (6 of 11) in patients obtaining PR led to significantly up-regulated expression of 4-1BB on CD8+T cells compared to the WT peptide controls (Fig. 3, A and B). However, no significant difference in 4-1BB expression on CD8+T cells was observed with the three predicted MANAs of the patients with progressive disease (fig. S5A). Besides, interferon-γ (IFN-γ) production detected by cell cytometry and enzyme-linked immunospot (ELISpot) showed elevated level in MANA-stimulating T cells (fig. S6). In the in vitro assays of MANA stimulation, the isolated PBMCs exhibited notable amplification and aggregation after incubation with the MANAs, suggesting the recognition of neoantigens by MANA-specific T cells (Fig. 3C). To further determine whether the in vitro amplified MANA-activated TCR clones were covered in the shared TCR repertoire used for the TIR index calculation, we performed comparative analyses of the TCR sequences among the TILs, the peripheral PD-1+CD8+T cells, and the MANA-stimulated 4-1BB–expressing CD8+T cells in patient nos. 2 and 3. In this analysis, all the detected TCR sequences in the in vitro stimulated T cells (MANAs 1, 2, and 4) and MANA4-tetramer–specific T cells were found in the overlapping TCR sequences between the TILs and the peripheral PD-1+CD8+T cells, which was not observed in patients with progressive disease and healthy donors (Fig. 3D and figs. S5B and S7). The gene mutations corresponding to MANAs 1, 2, and 4 referred to CDKN2A p.D74Afs*72, c.2TP53 p.G244V, and KRAS p.G12D, respectively. We further predicted and listed all the candidate neoantigen peptides with median inhibitory concentration (IC50) values of <500 nM. The numbers were significantly higher in the patients with a high TIR index than in those with a low TIR index (fig. S5C). These results indicated that the shared TCR repertoire contained tumor neoantigen–specific TCR clonotypes, and thus, the corresponding TIR index could be a functional index correlated with ICI therapy.

Fig. 3 In vitro validation of functional T cell stimulated by candidate MANAs.

(A) Representative flow cytometry analysis for the expression of 4-1BB on CD8+ lymphocytes after coculture with the MANA (left), WT (middle), and no peptide (right). The percentage of the cells expressing different combinations of receptors is shown. (B) Candidate MANA stimulation in patients obtaining PR. Adjusted P values are shown for pairwise comparisons. Blue bars indicate mutant peptides, orange bars indicate WT peptides, and gray bars indicate no peptides. (C) The algorism of in vitro coculture of immune cells with MANAs/WT peptide/no peptide. T cells were isolated from PBMCs (top left) and cocultured with DCs carrying MANA peptides or control peptides (middle). Notable T cell expansion was observed after 7 days of coculture with MANA (right) comparing to WT peptide and no peptide. The representative microphotographs are shown. (D) Comparative analyses of the TCR sequences among TILs, peripheral PD-1+CD8+T cells, and MANA-stimulated T cells. In patient nos. 2 and 3, all the detected TCR sequences in the in vitro stimulated T cells (MANAs 1, 2, and 4) were found in the overlapping TCR sequences between TILs and peripheral PD-1+CD8+T cells, which was not observed in healthy donors. *P < 0.05 and **P 0.01; n.s., not significant.

The TIR index was associated with the degree of similarity of the TCRs

A total of 1471 kinds of shared TCR clonotypes between the TILs and the peripheral PD-1+CD8+T cells were identified across all 31 patients. Theoretically, epitope-specific amino acid residues could be identified by a series of TCR sequences with a high degree of similarity (21). This finding indicated that the fraction of TCRs with few CDR3 differences within a patient would be relatively high due to similar genomic backgrounds. By calculating the edit distance of individual TCRs, we found that 2.6% of the shared TCRs had relatively close distances (0 to 3 defined as low) (Fig. 4A). More concisely, the patients obtaining a response or with a high TIR index showed a high fraction of TCRs with a close edit distance in the total shared TCRs compared with the nonresponders (3.8 versus 2.2%, P = 0.048; Fig. 4B) or those with a low TIR index (3.6 versus 2.1%, P = 0.034; Fig. 4C). As expected, the shared TCR sequences were more similar in the patients with a high TIR index than those with a low TIR index (Fig. 4, D and E). Similarly, the shared TCR clonotypes showed significantly higher clonality than the nonshared TCR clonotypes (P = 0.034; Fig. 4F). As validation, a lower fraction of remote edit distance (defined as up to 30%, ≥10) was obtained in the patients with a high TIR index than in those with a low TIR index (22.08 versus 31.45%, P = 0.041; fig. S8). Together with the positive correlation of the TIR index with response, these observations indicated that a high TIR index represented the enrichment of similar TCRs that were more likely to effectively recognize neoantigens.

Fig. 4 Association of TIR index with the degree of TCR similarity.

(A) Heatmap based on the edit distances of TCRs. The edit distances were calculated between every two TCRs extracted from 1471 shared TCR clonotypes in 31 enrolled patents. The blue color indicates lower values of the edit distance, which means that the two TCR clones are more similar in CDR3 regions. (B) Proportions of TCRs with close edit distance (0 to 3) between responders (red points, n = 6) and nonresponders (gray points, n = 25). (C) Proportion of TCRs with close edit distances (0 to 3) between high–TIR index group (red points, n = 7) and low–TIR index group (gray points, n = 24). (D and E) Representative examples of the edit distance distribution (top) and alignment of CDR3 sequences from a single cluster (bottom) in patient no. 7 with high TIR index (D) and patient no. 21 with low TIR index (E). (F) Comparison of the clonality between shared and nonshared TCR clonotypes in the enrolled patients (n = 31).

The TIR index was associated with immune activation in the TME

To better demonstrate that the TIR index is a functional index that can predict the efficacy of ICIs, we further analyzed the immune microenvironment characteristics across different TIR index subgroups. Cytokines reflect the immune conditions in mediating immune cell trafficking and lymphoid tissue development (22). The levels of 10 immunoregulation–related cytokines were assessed in our 31-patient cohort: those secreted mainly by T cells that can promote T cell growth and activation [interleukin-3 (IL-3), IL-4, IL-5, and IL-6], those secreted by macrophages and fibroblasts that promote T cell survival and chemotaxis (IL-7, IL-8, IL-9, and IL-12), and those that regulate the trafficking of monocytes, T cells, and neutrophils for immune surveillance (SDF-1 and CTACK). On the basis of these results, higher levels of the abovementioned cytokines were observed in the high–TIR index subgroup than in the low–TIR index subgroup (Fig. 5A).

Fig. 5 The association of TME with TIR index.

(A) Distribution of indicated cytokines in plasma between high–TIR index subgroup (red points, n = 7) and low–TIR index subgroup (gray points, n = 24). *P < 0.05 and **P < 0.01. (B) Heatmap clustered by the mRNA expression of 395 immune-related genes in 14 tumor samples (high TIR index, n = 4 and low TIR index, n = 10). The column and row represent individual patient and gene respectively, which have been grouped using unsupervised clustering. The expression levels were normalized by z score and standardized (centered and scaled) within rows for visualization. (C) Expression differentiation validated by an 18-gene panel involved in T cell inflammation in 14 tumor samples (high TIR index, n = 4 and low TIR index, n = 10). Expression levels have been normalized by z score within columns for visualization. The four columns on the left side indicate patients with high TIR index. The rows and columns have been grouped using unsupervised clustering. (D) Distribution of T cell–inflamed GEP score in the high–TIR index subgroup (red points, n = 4) and low–TIR index subgroup (gray points, n = 10). (E) Comparisons of fractions of M1 and M2 between high–TIR index and low–TIR index subgroups. The fractions of M1 and M2 were assessed by CIBERSORT based on the RNA expression data. The red points and gray points indicate patients with high TIR index (n = 4) and low TIR index (n = 10), respectively. (F) Proportions of CD8+T cells between the high–TIR index (red points, n = 4) and low–TIR index (gray points, n = 10) subgroup. (G) Representative microphotographs of mIF from patient no.3 with high TIR index (top) and patient no. 23 with low TIR index (bottom). An enlarged image is shown on the right. (H) Proportions of PD-1+CD8+T cells between high–TIR index (red points, n = 7) and low–TIR index (gray points, n = 24) subgroup. (I) Association of the fraction of TIR index with PFS in all enrolled patients (n = 31).

The GEPs for 14 enrolled patients, consisting of 4 patients in the high–TIR index subgroup and 10 patients in the low–TIR index subgroup, were then determined on the Oncomi-Seq platform. The heatmap of unsupervised clustering of 395 genes (table S6) clearly distinguished the TIR index subgroups (Fig. 5B). We further analyzed the GEP scores of 18 genes involved in T cell–inflamed gene expression, including IFN-γ–responsive genes related to antigen presentation, chemokine expression, cytotoxic activity, and adaptive immune resistance (8). Up-regulated expression of these genes was observed in the high–TIR index subgroup (Fig. 5C), while the GEP score was substantially higher in the high–TIR index subgroup than in the low–TIR index subgroup (median, 1.2 versus −0.57) (Fig. 5D). The abundance of 22 immune cell types within the TIL population was estimated according to the normalized read counts of related genes by using a deconvolution approach, CIBERSORT, which uses a validated leukocyte gene signature matrix (LM22) of 547 genes (fig. S9A). The fraction of macrophage 1 (M1) was positively correlated with the TIR index (median, 10.6 versus 5.3%), while the percentage of macrophage 2 (M2) was lower in the high–TIR index subgroup than the low–TIR index subgroup (median, 6.3 versus 17.6%) (Fig. 5E). The percentage of CD8+T cells was relatively higher in the high–TIR index subgroup than in the low–TIR index subgroup (23.4 versus 11.6%; Fig. 5F).

In addition, mIF was performed to detect the expression status of CD4, CD8, PD-1, PD-L1, LAG3, and TIM3 on immune cells in 31 enrolled patients (Fig. 5G). The coexpression of PD-1 and CD8, which was used to represent the fraction of PD-1+CD8+T cells, displayed a significant correlation with the TIR index level (P = 0.027; Fig. 5H); however, it failed to distinguish patients that showed an increase in PFS, although a trend was observed (HR = 0.63, 95% CI, 0.30 to 1.32, P = 0.22; Fig. 5I). The percentage of CD8 T cells with coexpression of TIM3 or LAG3 showed no difference in the high–TIR index and low–TIR index groups (P = 0.24 and 0.59, respectively; fig. S9B).

Dynamic evolution of the shared TCR repertoire during anti–PD-(L)1 treatment

To examine the dynamic evolution of the shared TCR repertoire clones during the anti–PD-(L)1 treatment, we analyzed the change in the frequency of the shared TCR clones and its correlation with the change in the maximum variant allele frequency (mVAF) and the efficacy of ICI treatment. The genomic profiles were assessed by NGS using a 1021-gene panel (table S4) in both the tumor tissues and the peripheral blood specimens (fig. S10, A and B). In 31 patients with matched pre- and on-ICI peripheral PD-1+ CD8+ TCR sequencing data, the frequencies of the shared TCR clones were well correlated between the peripheral PD-1+ CD8+ TCR clones (namely, the circulating shared TCR clone frequency) and the tumor TILs (r2 = 0.268, P = 0.002; Fig. 6A and fig. S10C). The frequency changes in the circulating shared TCR clones displayed a consistent trend (increased versus decreased) with the response and PFS (superior versus inferior; Fig. 5B) for anti–PD-(L)1 treatment (Fig. 6B, top). However, these changes were negatively associated with the mVAF of somatic gene mutation clones (Fig. 6B, bottom). In addition, patients with increased frequencies of circulating shared TCR clones demonstrated a longer PFS than those with decreased frequencies (median, 8.9 versus 2.0 months; HR, 0.28; 95% CI, 0.09 to 0.86, log-rank P < 0.001) (Fig. 6C). The frequencies of the shared TCR repertoire clones in the peripheral PD-1+ CD8+ TCR clones were significantly increased after 2 cycles of ICI treatment in the response group comparing to the nonresponse group (median, 0.51 versus −0.12%, P = 0.002) (Fig. 6D). These results reflected the expansion of specific clonal populations in ICI responders.

Fig. 6 Dynamic evolution of the shared TCR repertoire during anti–PD-(L)1 treatment.

(A) Pearson correlation of the frequencies of the shared TCR clones between peripheral PD-1+CD8+T cells and TILs. Shaded area, 95% CI for the correlation. (B) Association of the efficacy of anti–PD-(L) treatment with the change of ctDNA abundance and the proportion of shared TCR clones in peripheral PD-1+CD8+T cells. A total of 19 patients were divided into two subgroups with increased (n = 10) and decreased (n = 9) proportions of shared TCR clones in peripheral PD-1+CD8+ T cells. The top waterfall plot shows the observed best response and corresponding PFS of anti–PD-(L)1 treatment. The bottom shows the corresponding change of mVAF in dynamic ctDNA. (C) PFS stratified by dynamic frequencies of the shared clone in peripheral PD-1+ CD8+ T cells. (D) Change of the proportion of the shared TCRs in peripheral PD-1+CD8+ T cells before and after anti–PD-(L)1 treatment in the responders and nonresponders. (E) The relative proportion of TCRs (top) and frequencies of ctDNA subclones (bottom) in peripheral PD-1+CD8+ T cells collected at six longitudinal time points in patient no.4 obtaining initial PR with a 9-month PFS. (F) Dynamic changes of the targeted lesion at the left lower lobe from baseline to drug resistance in patient no. 4.

To further elucidate whether the changes observed in the shared TCR repertoires during anti–PD-(L)1 therapy were primarily a functional reflection over the selection and expansion of T cell clonotypes, we then conducted longitudinal peripheral blood–based dynamic tracking of the shared TCR clones in patient no. 4, who achieved PR upon atezolizumab treatment. The frequency of the circulating shared TCR clone initially increased when the patient showed PR at the first time and remained stable along with tumor shrinkage and continuously decreasing mVAF (Fig. 6, E and F, and fig. S10D). The tumor began to progress after 9 months of atezolizumab treatment; however, the shared TCR clone frequency did not correspondingly increase (Fig. 6, E and F, and fig. S10D). Moreover, analysis of the circulating tumor (ct)DNA-based subclone showed that the subclonal clusters presented different evolution patterns during atezolizumab treatment, and the enlarged subclones (subclones 6 and 7) at tumor progression that were considered major clones were actually minor clones before atezolizumab treatment (Fig. 6F), suggesting that tumor neoantigens due to these subclones could not be recognized by the initially activated tumor-specific TCRs, leading to subsequent tumor progression. This observation may explain why the seemingly “gained” mutations in resistance to anti–PD-(L)1 treatment were less likely to be associated with an antitumor immune response than other mutations.


The clinical and economic burden of the large proportion of patients who do not benefit from extensive administration of anti–PD-(L)1 drugs have led to a search for more effective matching of patients to these therapies. We herein developed a previously unidentified TCR repertoire–based predictor, which was indicative of TME status and could effectively identify patients with NSCLC who benefit from anti–PD-(L)1 treatment.

The PD-1+ CD8+ phenotype is a crucial marker of exhausted T cells with the highest possibility of containing neoantigen-specific cytotoxic T cells, which are the main targets of anti–PD-(L)1 therapy (9). A high proportion of PD-1+CD8+T cells in the overall TILs, reflected by TCR repertoire assessment more precisely, could represent an increased possibility of the enrichment of functional antigen-specific TILs. However, it is difficult to sequence the specific TCR repertoires of tumor-infiltrating PD-1+ CD8+ T cells directly. Circulating reinvigorated PD-1+ CD8+ T cells after inhibition of the PD-1 pathway contained both tumor-specific and bystander T cell subpopulations (15, 23). The shared TCR sequences between circulating PD-1+ CD8+ T cells and tumor regional TILs can be used to represent the tumor-specific PD-1+ CD8+ TCR clones, which was partially supported by a previous report showing that higher overlapping degrees of TCR in tumors and PBMCs were correlated with an improved immune response (24). In comparison, the TIR index established in our study adopted PD-1+ CD8+ T cells rather than PBMCs for TCR repertoire assessment, which could reduce the influence of non–PD-1+ CD8+ T cells in estimating antitumor immune effects. A high TIR index was found to be associated with a superior efficacy of ICIs, which was possibly due to the enhanced chance of encountering neoantigen-harboring tumor cells upon the massive release of tumor-specific cytotoxic T cells into peripheral blood. In addition, the in vitro stimulated TCR clones by candidate MANAs were successfully sequenced in both the shared TCR clones and the tumor regional TILs, strongly indicating that tumor-specific cytotoxic T cells are well represented by the shared TCR clones that were involved in the TIR index.

The cold immune microenvironment is considered an important correlate of poor efficacy to ICIs (25), especially for patients with lung cancer harboring specific gene alterations, such as EGFR mutations (18), MDM2/4 amplification, DNMT3A mutations (19), and JAK 1/2 mutations (20), which represent a large portion of lung malignancies. Several other biomarkers have been assessed for clinical utility, including PD-L1 ≥ 25% (26) and contaminant KRAS and TP53 mutations (27). However, very limited improvement has been achieved to date. Our preliminary evidence showed that the TIR index could predict survival benefits of patients with the above gene mutations, including activating EGFR mutations. Patients with a high TIR index had an ORR of 75%, along with elevated levels of immune regulation–related cytokines, T cell–inflamed GEP scores, and fractions of tumor regional PD-1+CD8+T cells. These observations confirmed that a higher TIR index could be a functional predictor for ICI delivery and an indicator of an activated TME.

Furthermore, although TCR recognition is the fundamental and critical step in the response to immunotherapy, the complexities of the TME and the immune system have suggested the need for combinations of multiple predictors. Previous studies have shown that high levels of PD-L1 expression and tTMB represented independent populations (28), and the two biomarkers could be complementary to each other for better predictions of ICI delivery. In the present study, excellent response rates from ICIs were observed when the TIR index was combined with PD-L1 expression or tTMB. All the patients with a high TIR index had positive PD-L1 expression, and only one-third of the patients with a high TIR index exhibited high tTMB levels. These results demonstrated that the utilization of the TIR index may appropriately broaden the population of patients likely to benefit with immunotherapy and highlight combined detection of multiple predictive biomarkers as a promising strategy, which warranted the further validation in future study design and can be potentially used in clinical practice.

ICI responders have been shown to undergo a process of T cell clone expansion, along with the elimination of neoantigenic mutation-defined tumor subclones (29). Acquired resistance for initial ICI responders would inferentially lead to tumor subclone growth without further activation of tumor-specific T cells. In our analysis, an opposite trend was observed for the change in the shared TCR clone frequencies against the mVAF of somatic gene mutation clones through longitudinal ctDNA detection. Notably, disease progression due to the expansion of initial minor tumor genetic subclones did not accompany synchronous elevation of the shared TCR clone frequency. This phenomenon further demonstrated that the shared TCR clones could be representative of tumor neoantigen–specific cytotoxic T cells. These observations together with the association of the frequency change of the shared TCR clones and the efficacy of anti–PD-(L)1 treatment suggested that the early detection and dynamic monitoring of the shared TCR clones may be complementary in the identification of patients benefiting from anti–PD-(L)1 therapy.

Despite the multiple-aspect validations followed by in vitro functional experiments, the TCR sequences identified by the current assay still did not accurately represent the functional T cells that are specific to an individualized tumor-specific neoantigen. Thus, TCR repertoire analysis of the TIR index provides a practical, precise, and quantitative approach to assess neoantigen recognition beyond T cell quantities. For strategies using an NGS panel for neoantigen prediction, there is a risk of missing some key candidate MANAs and somatic mutations, the corresponding peptides of which could activate autologous T cells. In addition, our TIR index can serve as an indicator for TME status, which could also differentiate the patients benefiting from anti–PD-(L)1 treatment to a different degree. However, the causality and the predominant effects of these two factors remain unclear, and further investigations of the underlying mechanisms are needed. Another limitation of our study was the relatively small sample size, which might lead to a potential bias. Nevertheless, multiple platforms were used for the assessments and analyses as cross validations, supporting the reliability of our results. Thus, future investigations are warranted to validate these exploratory findings in larger prospective cohorts. In summary, this study systematically demonstrated that weighting the tumor-specific TCR repertoire could be a clinically feasible classifier and provided a previously unidentified promising approach to stratify patients with NSCLC treated with anti–PD-(L)1 treatment, providing new insights for the development of future TCR-based investigations.


Study design and patient population

A total of 31 consecutive histologically confirmed patients with NSCLC diagnosed with stage IV disease were recruited and treated with anti–PD-(L)1 monotherapy from 14 March 2017 to 2 May 2018. All patients provided tumor tissue and paired peripheral blood specimens within 2 weeks before delivery of PD-1/PD-L1 inhibitors. Dynamic peripheral blood samples were prospectively collected every 6 to 8 weeks at the same time of image evaluation, and PD-1+ CD8+ T cells were sorted through flow cytometry for TCR sequencing (fig. S11). Cell-free DNA (cfDNA) was extracted for genetic analysis. Partial spare tissue specimens were used for mIF and gene expression profiling analyses. Blood samples used for MANA/WT peptide stimulation were obtained from two healthy donors after written informed consent was provided. This study was approved by the ethics committees of the National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences and the Peking Union Medical College (NCC2016JZ-03 and NCC2018-092), and all enrolled patients provided written informed consent.

Isolation of PD-1+ CD8+ T cells and identification of T cell activation markers

PBMCs were isolated from 10 ml of fresh peripheral blood containing an anticoagulant by density gradient centrifugation using LymphoPrep (Progen, Heidelberg, Germany). On the basis of fluorescence-activated cell sorting (FACS) analysis (BD FACSAria II), PD-1+ CD8+ T cells were purified from PBMCs. The following human protein-specific flow cytometric antibodies were purchased from eBioscience: CD3-efluor 450 (clone: OKT3), CD8-APC (clone: RPA-T8), IFN-γ-percp-cy5.5 (4S.B3), and CD279 (PD-1)-PE (MIH4). MANA from patient no. 3 was loaded at 10 mM onto the QuickSwich Quant HLA-A*11:01 Tetramers (PE) (MBL International) according to the manufacturer’s instructions. CD8+T cell activation was assessed by staining for 4-1BB (CD137), MANA-specific tetramer, and IFN-γ.

Next-generation sequencing

Circulating cfDNA was isolated from plasma samples using a QIAamp Circulating Nucleic Acid Kit (Qiagen, Hilden, Germany). Genomic DNA from PBMCs and tissue samples was extracted using the QIAamp DNA DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). All DNA extractions were performed according to the manufacturer’s instructions. For detection of genomic alterations (single-nucleotide variants, insertions, deletions, and amplifications) in plasma and tissue biopsies, an NGS-based 1021-gene panel with potential clinical relevance was used to capture target regions. For ctDNA analysis, molecular barcoding and proprietary bioinformatic algorithms were performed. DNA sequencing was carried out with paired-end reads on an Illumina HiSeq 3000 platform in a CLIA/CAP-accredited laboratory (Geneplus) (30). The mVAF was used to monitor dynamic changes in ctDNA in each patient.

Neoantigen predictions

According to the germline sequencing data from the PBMCs, the four-digit HLA type was first inferred using OptiType. Then, somatic mutations, including nonsynonymous single base substitutions, insertions, and deletions, were used to perform a comprehensive assessment of peptide 8 to 11 amino acids in length at every position surrounding a somatic mutation. The major histocompatibility complex (MHC) class I binding potential of each somatic and WT peptide was predicted by NetMHCpan 4.0. Neoantigen candidates with a high MHC class I binding affinity (IC50 value < 500 nM) were further characterized for tumor-associated expression levels derived from TCGA to generate a final ranking of peptides for experimental follow-up.

Peptide-MHC tetramer production

The peptide exchange experiment was performed with QuickSwitch Quant HLA-A*11:01 Tetramer Kit-PE (MBL) according to the manufacturer’s instructions. Briefly, dissolve each lyophilized peptide in dimethyl sulfoxide at a stock concentration of 10 mM. Fifty microliters of QuickSwitch Tetramer was pipetted into a microtube, and 1 μl of target peptide was added and mixed gently with pipetting. Then, 1 μl of peptide exchange factor was added and mixed gently with pipetting. The samples were incubated for at least 4 hours at room temperature protected from light.

Dendritic cell and T cell preparation and in vitro expansion of MANA-stimulated T cells

T cell functional assays were performed for four patients with NSCLC with PR and two patients with progressive disease treated with PD-1/PD-L1 antibody. Fifteen pairs of MANAs and WT peptides were synthesized (GL Biochem, Shanghai Ltd.) and used to stimulate T cells in vitro.

PBMCs were collected from the patient peripheral blood samples and isolated using Ficoll-Paque PLUS (GE Healthcare, USA). The total cells were cultured in AIM V serum-free medium (Gibco, Grand Island, NY) and allowed to adhere for 4 hours. The adherent cells were cultured using DC serum-free medium (CellGenix, Germany) containing IL-4 (1000 U/ml) and granulocyte-macrophage colony-stimulating factor (GM-CSF; 500 U/ml). Following next 6 days for facilitating DC growth, another 10 ng/ml of tumor necrosis factor–α (TNF-α) was added to the DCs to induce maturation on the sixth day. The nonadherent cells were collected at the first day and induced to become T cells followed by stimulating with anti-CD3 antibody (50 ng/ml) and recombinant human IFN-γ (1000 U/ml; Peprotech, USA) at 37°C with 5% CO2 for 24 hours. Then, recombinant human IL-2 (1000 U/ml; Peprotech) was added to the medium. IL-2–containing medium was added to the culture system every 2 days.

For T cell activation, after 7 days of T cell culturing, the T cells were mixed with DCs loaded with 10 mM peptide at ratio of 20:1 and cocultured using conditioned medium supplemented with IL-2 (1000 U/ml), IL-4 (500 U/ml), GM-CSF (250 U/ml), and TNF-α (10 ng/ml) (Peprotech) for another 7 days. As mature DC cells have a short life, most DCs in the coculture system died when the T cells were collected. The T cells were harvested and analyzed by TCR sequencing and flow cytometric analysis. In terms of TMB, the top quarter (14 muts/MB) is used as the cutoff value.


IFN-γ–secreting T cells stimulated by MANA/WT peptide/scramble peptide/no peptide were detected by Human IFN-γ ELISpot assay kits (3512-4APW-2, Mabtech) according to the manufacturer’s protocol. PBMCs were plated in triplicate at 3 × 105 per well and then incubated for 24 hours. Spots were then counted using an S6 ultra immunoscan reader (Cellular Technology Ltd.), and the number of IFN-γ–positive T cells was calculated by ImmunoSpot 5.1.34 software (Cellular Technology Ltd.).

TCR-Vβ deep sequencing

DNA was extracted from longitudinally sorted PD-1+CD8+ T cells, baseline tumor tissues, and peptide-stimulated T cells using a Qiagen DNA blood kit or DNA FFPE (formalin-fixed paraffin-embedded) kit. TCR-Vβ sequencing was performed in a CLIA/CAP-accredited laboratory (Geneplus) as described previously (14). The sequences of the V, D, and J genes were compared with the ImMunoGeneTics (IMGT) database by MIXCR software, and the CDR3 sequence of TCR was identified. Then, 106 qualified reads were randomly selected for downstream analysis. Bioinformatic and biostatistical analyses were performed to identify antigen-specific expansions using the following criteria: (i) significant expansion (Fisher’s exact test with Benjamini-Hochberg false discovery rate, P < 0.05) compared to T cells cultured without peptide and (ii) no significant expansion in any other peptide-stimulated culture.

The diversity of the TCR repertoire was calculated on the basis of the Shannon-Wiener index (Shannon index), which scaled from 0 to 10: Shannon index H=i=0npilnpi, where pi refers to the frequency of clonotype i for the sample with n unique clonotypes. Clonality can be considered as a normalized Shannon index over n unique clones: clonality = 1 − H/ln(n).

To assess the similarity of TCR repertoires between samples, we used the metric of TCR repertoire overlap. Specifically, if n TCR clones were identified in both TILs (T) and baseline peripheral PD-1+CD8+ T cells (P), the sum frequency of these TCR clones was calculated, and the TIR index was defined as the average frequency of overlapping clones, with the formula TIR index=i=1nFreq(TP)iTn ranging from 0 to 1, i indicates the frequency of clonetype i among n overlapped clones. The follow-up tracking is based on the dynamic frequency of the overlapping clones in the peripheral blood PD-1+CD8+TCR during treatment.

RNA extraction and gene expression analysis

Total RNA was extracted from 5-μm-thick FFPE sections of tumors fixed on positively charged slides using the Ambion RecoverAll kit for RNA isolation from FFPE tissue (Thermo Fisher Scientific) following the manufacturer’s protocols. Library preparation was conducted with the Oncomine Immune Response Research Assay, and sequencing was performed by a Thermo Fisher Ion S5 System. A custom code set consisting of a 395-gene panel related to T cell biology, immune regulation, and cellular markers of TILs and tumor-associated macrophages was used. Preliminary bioinformatic analysis was performed by using the Thermo Fisher Torrent Suite Immune Response RNA plugin.

The analytical tool CIBERSORT was used in the study; this tool can accurately quantify the percentage of different types of tumor infiltrating cells (TICs) under the complex “gene signature matrix” based on 547 genes (31). In the current study, we determined the immune infiltration of each sample with the LM22 signature file, which can define 22 subtypes of immune cells—including naïve B cells, memory B cells, plasma cells, CD8+ T cells, naïve CD4+ T cells, resting memory CD4+ T cells, activated memory CD4+ T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, resting natural killer (NK) cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting DCs, activated DCs, resting mast cells, activated mast cells, eosinophils, and neutrophils—with the preset signature matrix at 1000 permutations. After analysis with the CIBERSORT program, the distribution of 22 subtypes in different TIR index groups was presented.

The gene expression data were also used to evaluate the T cell–inflamed signature, an 18-gene signature defined by Ayers et al. (8). Starting from a tumor dataset spanning nine cancer types (KEYNOTE-012 and KEYNOTE-028 studies), the authors defined the T cell–inflamed GEP as an 18-gene signature, and based on this improved set of genes, they determined the T cell–inflamed signature score as a linear function defined as follows: TIS = Σi = 118xiwi, where xi is the ith gene’s expression value obtained as described above and wi is a predefined weight for the ith gene derived by Ayers et al. (8), which were derived from the logistic regression model.

Multiplex immunofluorescence

FFPE tissues were cut into 4-μm-thick sections and placed on plus-charged slides. Slides were heated at 57°C overnight; then, residual paraffin was removed using xylene, and the tissue was rehydrated in a series of graded alcohols in distilled water. Antigen retrieval was performed using tris-EDTA buffer and microwave treatment. The slides were washed, and blocking was performed with 3% H2O2 blocking solution followed by Dako antibody diluent. The first primary antibody, i.e., position 1 in the table above, was then applied and allowed to incubate. Opal polymer horseradish peroxidase Ms + Rb (PerkinElmer, Hopkinton, MA) was used as the secondary antibody. The slides were washed, and 4 tyramide signal amplification dye (Opal 7 color kit, PerkinElmer, Hopkinton, MA) for position 1 was applied. The slides were then microwaved to strip the primary and secondary antibodies, washed, and blocked again using blocking solution. The second primary antibody, i.e., position 2, was applied, and the process was repeated through position 6, where 4′,6-diamidino-2-phenylindole (DAPI) was applied, rather than another primary antibody. After unbound DAPI was washed off, the slides were coverslipped using Vectashield (Vector Laboratories, Burlingame, CA). In addition to the multiplex assay described above, a single-color slide was generated for each antibody using an archival Merkel cell carcinoma case. Each single-color control slide was compared with standard immunohistochemical methods for CD8, CD68, PD-L1, PD-1, and NSE (neuron-specific enolase).

Multiplexed and single-color control slides were loaded onto the PerkinElmer Vectra automated multispectral microscope. Representative fields from the single-color slides were imaged, and InForm Image Analysis software (version 2.1) was used to generate a spectral library for unmixing. Index cases stained using the multiplex method were then imaged. Channels were unmixed using the spectral library, and tissues and cells were segmented and scored.

Cytokine analysis

Plasma was collected for cytokine protein level detection via Meso Scale Discovery (MSD) technology using a U-PLEX Biomarker Group 1 (Human) Multiplex Assay (K15049D-1). A MESO QuickPlex SQ 120 machine and MSD Discovery Workbench software version 4.0.12 were used for the detection and data analyses, respectively.

Treatment evaluation

Radiographic imaging was acquired by the indicated approaches, such as CT or MRI, for tumor response assessment, which was evaluated by both the investigator and an independent radiologist. Baseline tumor assessments were performed within 1 to 28 days before the initiation of anti–PD-1/PD-L1 treatment, with subsequent assessments performed every 4 to 8 weeks until objective disease progression. The efficacy was determined by the Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1, including CR, PR, SD, and progressive disease. Disease control indicated the patients who achieved disease control (i.e., CR, PR, or SD according to RECIST version 1.1). PFS was defined as the time from the start of anti–PD-1/PD-L1 treatment until either objective disease progression (assessed by an investigator using RECIST version 1.1) or death from any cause. OS was defined as the time from the start of anti–PD-1/PD-L1 treatment until death from any cause. We censored patients who had not progressed at the time of statistical analysis or who were lost to follow-up before death at the time of their last evaluation.

Statistical analysis

The optimal thresholds for the TIR index were calculated by ROC analysis using a nonparametric bootstrap. The Youden index (J = sensitivity + specificity − 1), as a commonly used measure of overall diagnostic effectiveness, was used in conjunction with ROC analysis. Independent predictive variables of PFS were identified by a multivariate Cox regression analysis using a backward selection procedure. The association between response and the high versus low TIR index was assessed by Fisher’s exact test. Kaplan-Meier curve analysis of PFS and OS was compared using the log-rank test. Statistical analyses were performed using Student’s t test and a Mann-Whitney U test for unpaired data, and the Wilcoxon rank sum test was used for pairwise analyses. All reported P values are two-tailed, and for all analyses, P < 0.05 was considered statistically significant, unless otherwise specified. IBM SPSS software (version 23.0), GraphPad Prism (version 8.4), and R 3.4.2 were used in the statistical analysis.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank all the patients who were involved in this study. We also thank the Geneplus-Beijing Institute, Beijing, China, for their support. Funding: This work was supported by the China National Natural Science Foundation Key Program (81630071 and 81330062) to J.W., National Key R&D Program of China (2016YFC0902300 and 2019YFC1315700 to J.W.; 2019YFC1315704 to Z.W.), CAMS Innovation Fund for Medical Sciences (CIFMS 2016-I2M-3-008 to J.W.; 2017-I2M-1-005 to Z.W.), Aiyou Foundation (KY201701 to J.W.), CAMS Key Laboratory of Translational Research on Lung Cancer (2018PT31035 to J.W.), National Natural Science Foundation of China (81871889 and 82072586 to Z.W.; 81972905 to J.D.; 82003275 to J.H.), and Beijing Natural Science Foundation (7212084) to Z.W. Author contributions: Concept and design: Z.W. and J.W. Development of methodology: J.H., W.Z., J.L., H.B., and Y.W. Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.H., R.Y., J.D., H.B., Y.W., X.Z., J.X., R.W., X.W., K.F., Z.Y., Z.W., and J.W. Analysis and interpretation of data (e.g., statistical analysis, biostatistics, and computational analysis): J.H., R.Y., J.D., H.B., G.F., Y.W., Y.G., X.X., K.F., Z.W., and J.W. Writing, review, and/or revision of the manuscript: Z.W., J.W., J.H., R.Y., and D.P.C. Administrative, technical, or material support (i.e., reporting or organizing data and constructing databases): J.H., Z.W., and J.W. Study supervision: Z.W. and J.W. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article