Comutations in DNA Damage Response Pathways Serve as Potential Biomarkers for Immune Checkpoint Blockade

Immune checkpoint blockades (ICB), including antibodies to programmed death receptor 1 (PD-1) or its ligand (PD-L1) and CTL-associated protein 4 (CTLA-4), are now standard therapies for a range of solid tumors after approval by the FDA (1, 2). However, only a subset of patients benefit from ICB monotherapy. Biomarkers such as PD-L1 expression (3, 4), tumor mutational burden (TMB; ref. 5), neoantigen load (NAL; ref. 6), tumor-infiltrating lymphocytes (TIL; ref. 7), and immune-regulatory mRNA expression signatures (8) are potentially applicable to the clinical selection of patients for ICBs, but each has limited utility.

Challenges in defining cut-off values, intratumoral heterogeneity distributions, test platform uniformities, and dynamic changes have limited the clinical application of PD-L1 expression (9). High TMB and NAL are associated with improved response to ICB treatment (5, 10), but the lack of a validated cut-off value has limited the use of this method (9). Mismatch repair deficiency (dMMR) may contribute to genome instability and lead to increases in TMB and NAL, which can be indicated by microsatellite instability-high (MSI-H) status (11), and thus, pembrolizumab was approved by the FDA for dMMR/MSI-H solid tumors (1). However, less than 5% of solid tumors are dMMR/MSI-H (12), which is the limiting factor for the clinical utilization of dMMR/MSI-H.

The DNA damage response (DDR) system comprises eight pathways: mismatch repair (MMR); base excision repair (BER); check point factors; Fanconi anemia; homologous recombination repair (HRR); nucleotide excision repair (NER); nonhomologous end-joining; and DNA translesion synthesis (13). The DDR system is essential for the preservation of genomic integrity (14), and thus, mutations in this system may induce a hypermutational phenotype. Previous studies have revealed that gene deficiencies in two DDR pathways of extraordinary fidelity, MMR (MLH1, MSH2, MSH6, and PMS2) and BER (POLE), resulted in a durable clinical benefit from ICBs (1, 15). Patients with melanoma responding to ICBs commonly harbored mutations in BRCA2, a major gene in the HRR pathway (6). Another study revealed that patients with deleterious alterations in 34 genes involved in DDR pathways exhibited high TMB levels and had improved clinical outcomes of ICB therapy in urothelial cancer (16). However, the number of DDR genes involved in these studies is limited, and most other DDR genes are poorly studied. In addition, due to the interactions among different DDR pathways, the contribution weights of different DDR pathways to TMB should be examined but have not been well established to date.

Defects in one DDR pathway lead to a greater dependency on the remaining DDR pathways, implying higher genomic instability when multiple deficiencies exist in different DDR pathways (17). Therefore, we hypothesized that mutations in multiple DDR pathways are associated with higher genomic instability; as a result, we would expect higher NAL and TMB. By analyzing whole-exome sequencing (WES) data from The Cancer Genome Atlas (TCGA) and validating the results using the International Cancer Genome Consortium (ICGC) database, we explored the optimized pattern of mutation combinations across distinct DDR pathways for TMB and NAL estimation. We herein report the relationships between comutations in DDR pathways and TMB, NAL, immune-regulatory mRNA expression signatures and the clinical benefit of ICBs.

Supplementary Table S1 summarizes the detailed information regarding the data used in this study. In brief, we obtained the WES data of 8,552 solid tumors from TCGA, the mRNA expression of 8,482 solid tumors across 29 tumor types from cBioPortal (www.cbioportal.org), and whole-genome sequencing data of 2,638 solid tumors from ICGC (dcc.icgc.org). The experimental procedures regarding DNA and RNA extraction from tumors, library preparation, sequencing, quality control, and subsequent data processing were published previously by TCGA (18). Neoantigens of 5,935 solid tumors estimated by TCGA data were downloaded from https://tcia.at/neoantigens (19). MSI of 5,930 solid tumors across 16 tumor types was estimated by Hause and colleagues with TCGA WES data (20).

To further explore the association between DDR pathways and the clinical benefit of ICBs, we included genomic and clinical data from four clinical cohorts treated with ICBs. The first cohort consisted of 34 patients with non–small-cell lung cancer (NSCLC) treated with anti-PD-1 therapy (Rizvi cohort; ref. 21). The second and third cohorts consisted of 110 and 64 patients, respectively, with advanced-stage melanoma treated with anti-CTLA-4 therapy (Allen cohort and Snyder cohort, respectively; refs. 10, 22). The last cohort comprised 38 patients with melanoma treated with anti-PD-1 therapy (Hugo cohort; ref. 6).

The detailed profiles of genes involved in eight DDR pathways are listed in Supplementary Table S2 (13). As most DDR genes have not yet been studied, we defined mutations in DDR pathways as any nonsilent mutations in gene-coding regions in the corresponding pathways, including missense, nonsense, insertion, deletion, and splice mutations. Comutations of HRR-BER or HRR-MMR were defined as “co-mut⁺.”

TMB was defined as the total nonsilent somatic mutation counts in coding regions. To explore how many DDR pathway mutations predicted higher TMB or NAL, TMB or NAL was classified into high/low based on the top quartile, and a receiver operating characteristic (ROC) curve was applied to determine the cutoff of pathway mutation counts with the highest Youden index.

The response to ICBs has been reported to be related to immune tumor microenvironments, including T-effector, IFNγ-associated genes, T-cell receptors, and immune factors (23, 24). The associations between DDR pathway comutation status and relevant immune-related genes were analyzed in 7,570 patients from TCGA for whom both RNAseq and DNAseq data were available. The immune gene list was mainly based on a published article that summarized the genes related to activated T cells, immune cytolytic activity, and IFNγ release (25). Other immune genes were added according to two relevant clinical trials (26, 27). A list of 47 immune genes is provided in Supplementary Table S3. The mRNA expression from cBioPortal was quantified by RSEM (RNAseq by expectation-maximization; ref. 28). The data were log₂-transformed before analysis.

Gene set enrichment analysis (GSEA) was performed using the javaGSEA 3.0 Desktop Application (http://software.broadinstitute.org/gsea/index.jsp) to identify whether immune-related gene signatures were associated with DDR pathway comutation status. The genes identified to be on the leading edge of the enrichment profile were subjected to pathway analyses. Genes with expression equal to 0 in more than 80% of samples were excluded from the GSEA. The normalized enrichment score (NES) is the primary statistic for examining gene set enrichment results.

Across four ICB-treated cohorts, we analyzed comutation status, TMB, and PD-L1 expression for association with patient prognosis. The TMB-high group was defined as the top quartile, and the TMB-low group was defined as the bottom three quartiles (29). PD-L1⁺ was defined as PD-L1⁺ in ≥ 50% of tumor cells, and PD-L1⁻ was defined as PD-L1⁺ in <50% of tumor cells (3, 4). We further tested the association of comutation status with survival outcomes stratified by TMB and PD-L1. Measures of patient survival and response were based on definitions consistent with how they were evaluated in the abovementioned trials. For the Rizvi cohort and the Allen cohort, response was defined by RECIST criteria. For the Hugo cohort, responding tumors were defined as exhibiting a complete response or a partial response or as stable disease, and nonresponding tumors were defined as exhibiting disease progression by irRECIST. For the Snyder cohort, response data were not available. We also analyzed the association between comutation status and patient prognosis in NSCLC and melanoma from TCGA, using propensity-score matching with a ratio of 1:2. Propensity score was estimated by age, sex, pathologic tumor stage, other history of malignancy, and histologic classification.

All the analyses were performed by R 3.4.2. Correlations between number of mutated DDR pathways and TMB or NAL were determined by Spearman rank correlation coefficient. Multiple linear regression with a stepwise method was used to examine the correlations among eight pathways and TMB or NAL. If NAL, TMB, or mRNA were normally distributed, Student t test was used to determine the differences between two groups; otherwise, the Mann–Whitney U test was used. Differences among three or more groups were determined by the Kruskal–Wallis test with Dunn posttest. Heatmap was used to depict the mean difference of immune-related genes between co-mut⁺ and co-mut⁻ subgroups by pheatmap package in R. Correlation-based group average hierarchical clustering was used for gene clustering and Minkowski distance for samples clustering. Survival analysis was performed using Kaplan–Meier curves, with a P value determined by a log-rank test. HR was determined through Cox regression. Proportional hazards assumption was tested before the Cox regression. To explore the modification effects of PD-L1 or TMB on the association between the comutation status and prognosis, stratified analyses were conducted, and interactions were evaluated with the log likelihood ratio statistic. Fisher exact test or the χ² test was used to identify the association of comutation status and objective response rate. All reported P values were two-tailed, and P < 0.05 was considered significant unless otherwise specified. A false discovery rate (FDR) < 0.05 was considered significant in the pan-cancer mRNA analysis and GSEA.

By analyzing the sequencing data from TCGA and validating the findings in ICGC, we identified that both TMB and NAL were significantly higher in patients with mutations in one or more DDR pathways than in those with wild-type status (Fig. 1A). TMB and NAL increased significantly with the number of mutated DDR pathways (Supplementary Fig. S1A–S1C) with Spearman correlation coefficients of 0.66 (for TMB in TCGA), 0.68 (for TMB in ICGC), and 0.59 (for NAL in TCGA). The TMB levels of the top quartile are commonly considered TMB-high (TMB-H; ref. 29), which were 128 and 88 counts of TMB in the TCGA and ICGC databases, regardless of cancer type. However, the median mutation counts for patients with mutations in a single DDR pathway were 52 and 53 in the TCGA and ICGC databases, respectively (Supplementary Fig. S1A and S1B), suggesting that mutation of a single DDR pathway may not account for TMB-H. The ROC curve showed that mutations covering ≥ 2 DDR pathways demonstrated an optimal Youden index, with 79.6% sensitivity and 80.9% specificity, compared with other numbers of pathway combinations in TMB-H estimation (Fig. 1B). In addition, similar results were verified in the ICGC TMB data, with 65.7% sensitivity and 84.8% specificity (Fig. 1C), and in TCGA NAL data, with 75.3% sensitivity and 78.8% specificity (Fig. 1D). These results confirmed our hypothesis that comutations in DDR pathways could predict higher TMB and NAL.

Given the potential interactions between different DDR pathways, to identify the most effective patterns of pathway combinations, we studied the associations of eight DDR pathways with TMB or NAL to select DDR pathways that contributed more to TMB or NAL, estimated by multiple linear regression analyses in patients with mutated DDR pathways limited to five. Together with the consistency between TMB and NAL data, three pathways, BER, HRR, and MMR, were selected with satisfactory standardized β coefficients (Fig. 1E). Although we could not fully exclude the functions of other comutation combinations, more specifically, patients with HRR-BER or HRR-MMR comutations illustrated significantly or borderline higher TMB and NAL than those with BER-MMR and other comutations (Fig. 1F). Thus, we defined comutations in the HRR-BER and HRR-MMR pathways as the “co-mut⁺” subgroup.

A total of 740 among 8,552 cases (8.7%) were classified as co-mut⁺ in TCGA dataset. Incidences of co-mut⁺ varied across different cancer types, ranging from 0% in pheochromocytoma and paraganglioma to 29.2% in bladder urothelial carcinoma (Fig. 1G). The frequencies of global DDR-related gene mutations in each cancer type were similar among the HRR, BER, and MMR pathways (Supplementary Fig. S2). We then evaluated whether co-mut⁺ status was associated with a higher level of TMB or NAL. The results showed significantly higher levels of TMB (Fig. 1H) and NAL (Fig. 1I) in the co-mut⁺ subgroups than in the co-mut⁻ subgroups across the cancer types with indications for ICBs delivery according to the FDA. Similar results were also obtained in other cancer types (Supplementary Tables S4 and S5). In addition, co-mut⁺ was significantly associated with MSI-H (P < 0.0001) and predicted MSI-H with 76.1% sensitivity and 98.8% specificity (Fig. 1J). In tumors with high frequency of MSI-H, such as colorectal adenocarcinoma, stomach adenocarcinoma, and uterine corpus endometrial carcinoma, most MSI-H samples were co-mut⁺, and importantly, co-mut⁺ covered more cancer types than MSI-H (Supplementary Fig. S3). No significant difference in TMB was observed between the co-mut⁺/MS stability (MSS) and co-mut⁻/MSI-H groups (Fig. 1K). These findings suggested that co-mut⁺ was a favorable surrogate for TMB estimation covering more percentage of patients than MSI-H.

To investigate the correlations between the co-mut⁺ subgroup and immune-regulatory gene expression signatures, we analyzed 7,570 samples from TCGA with both RNAseq and WES data. The GSEA revealed prominent enrichment of mRNA signatures involved in IFNγ pathway, inflammatory response, and allograft rejection in the co-mut⁺ subgroup over a pan-cancer analysis (Fig. 2A). Among the 47 selected immune-related genes, the mRNA expression of 31 genes was significantly higher in the co-mut⁺ subgroup than in the co-mut⁻ subgroup (FDR < 0.05) across all cancer types. The mean differences of these 31 genes between the co-mut⁺ and co-mut⁻ subgroups are shown in the heatmap (Fig. 2B). Specifically, the co-mut⁺ subgroup demonstrated higher levels of mRNA expression than did the co-mut⁻ subgroup in the following gene clusters: immune checkpoint (PD-L1 and LAG3; Fig. 2C); tumor immune microenvironment (IL1B, IL6, and PTGS2; Fig. 2D); T-effector and IFNγ-associated signatures (INFG, CXCL9, CXCL10, CXCL11, EOMES, GZMB, GBP1, and STAT1; Fig. 2E); and T-cell receptor–related genes (CD8A, CD3G, CD3D, IKZF3, and TIGIT; Fig. 2F). In addition, in 14 cancer types, at least one immune-regulatory gene had an expression level significantly higher in the co-mut⁺ subgroup than in the co-mut⁻ subgroup (Supplementary Fig. S4).

We further investigated the relationship between co-mut status and the clinical outcomes of ICBs. Four publicly available datasets were utilized for this analysis, including one NSCLC cohort and three melanoma cohorts. As in the above pan-cancer analyses, co-mut⁺ status was significantly positively associated with TMB and NAL in these four cohorts (Supplementary Fig. S5A–S5C).

The Rizvi cohort (21) included 34 patients with NSCLC treated with pembrolizumab, among whom 11 patients were classified as co-mut⁺; their progression-free survival (PFS) was superior to that of those without co-mut⁺ [median, not reached (NR) vs. 4.1 months; HR, 0.21; 95% CI, 0.06–0.71; P = 0.006, Fig. 3A]. The objective response rates (ORR) in co-mut⁺ and co-mut⁻ were 63.6% and 21.7%, respectively (P = 0.02, Fig. 3B).

We performed a pooled analysis on two metastatic melanoma cohorts [the Allen cohort (10) and the Snyder cohort (22) treated with CTLA-4 antibody, n = 174]. A total of 50 patients were classified as co-mut⁺, and their overall survival (OS) was significantly longer (median, 32.4 vs. 10.8 months; HR adjusted for cohort, 0.64; 95% CI, 0.42–0.98; P = 0.04; Fig. 3C). However, no difference was observed in tumor response rate (43.8% vs. 29.5%, P = 0.15, Fig. 3D). The Hugo cohort (6) enrolled 38 patients with advanced-stage melanoma who were treated with anti-PD-1 therapy. Co-mut⁺ patients demonstrated numerically superior OS (median, NR vs. 27.9 months; HR, 0.52; 95% CI, 0.18–1.44; P = 0.20; Fig. 3E) and ORR (73.3% vs. 43.4%, P = 0.10, Fig. 3F) compared with those of co-mut⁻ patients. In addition, no significant difference was found in survival outcomes between the single DDR pathway mutant and DDR-wild-type groups in these four public cohorts (Supplementary Fig. S6A–S6C), further supporting the hypothesis that a single mutated DDR pathway was not sufficient to predict greater benefit from ICBs.

To verify whether co-mut⁺ was a predictive or prognostic factor, we performed survival analyses according to comutation status in TCGA using propensity-score matching. As a result, no significant difference was found in OS between the co-mut⁺ and co-mut⁻ subgroups in either the patients with NSCLC or melanoma (Fig. 4A and B), suggesting that co-mut was a predictor of the clinical benefit of ICBs instead of a prognostic marker.

The predictive power of comutation status was further compared with PD-1 expression or TMB in the clinical cohorts. TMB calculated by WES with the top quartile or median as a cut-off value failed to consistently predict survival improvement by ICBs and thus was apparently inferior to comutation status as a predictive marker (Supplementary Table S6). Furthermore, in TMB-low (defined as less than the top quartile) NSCLC patients of the Rizvi cohort, both the median PFS and ORR of co-mut⁺ patients were significantly superior to those of co-mut⁻ patients (median NR vs. 3.8 months, P = 0.04, Fig. 5A left; ORR: 100% vs. 18.2%, Fig. 5B left). There was no difference between the co-mut⁺ and co-mut⁻ subgroups in terms of median PFS or ORR in TMB-high patients with NSCLC (Fig. 5A and B, right). No interaction between co-mut and TMB was observed (P = 0.16). In the melanoma cohorts, co-mut⁺ status also predicted improved OS and ORR in the TMB stratification analysis, although significant P values were not reached (Supplementary Fig. S7A–S7D).

In patients with NSCLC with PD-L1 < 50%, co-mut⁺ patients had significantly prolonged PFS (median, 8.3 vs. 2.7 months; HR, 0.26; 95% CI, 0.07–0.93; P = 0.03; Fig. 5C, left) and higher ORR (33.3% vs. 7.1%, P = 0.02; Fig. 5D, left) compared with those without co-mut⁺. In patients with NSCLC with PD-L1 ≥ 50%, co-mut⁺ status also predicted a trend of longer PFS (median, NR vs. 14.5; HR, NA; P = 0.08; Fig. 5C, right) and higher ORR (100% vs. 50%; Fig. 5D, right) than co-mut⁻. No significant interaction was observed between comutation status and PD-L1 expression (P = 0.99).

In this study, we demonstrated that the presence of comutations in the HRR-MMR and HRR-BER pathways was associated with higher TMB, NAL, and immune-regulatory gene expression and is a potential predictive biomarker for ICBs. The co-mut⁺ status identifies a patient population that may potentially benefit from ICBs irrespective of PD-L1 expression and TMB status.

The DDR system is required to preserve genomic integrity and has multiple functional links with innate immune signaling to protect against pathogens such as damaged DNAs (11). Functional mutations in the DDR system reduce genomic stability (17), although it remains unclear which DDR pathway contributes more to genomic stability and how their interactions may impact clinical outcomes. We found that comutations involving two pathways, HRR-BER and HRR-MMR, are the key reasons for the elevation of TMB and NAL. The associations between these comutations and increased TMB and NAL are consistent with previous reports based on WES data regarding genomic and transcriptomic features related to ICB response (6). The HRR and BER pathways monitor and repair DNA double-strand breaks (DSB) and single-strand breaks (SSB), respectively (30), whereas the MMR pathway is responsible for DNA insertion/deletion corrections (31). We speculated that high TMB and NAL were the result of increasing DDR mutation accumulation. Defects in the HRR pathway may lead to an inability to repair the DNA damage caused by BER or MMR deficiency. The underlying mechanism of coordination between HRR-BER and HRR-MMR should be further investigated in the future.

Comutation is a potential predictive biomarker for the use of PD-1 blockades in patients with advanced NSCLC. Both the tumor response rate and median PFS were improved in the co-mut⁺ subgroup. In melanoma, the median OS was longer in the co-mut⁺ subgroup, while the difference in ORR was not significant. Approximately 36.8% (14/38) of the patients in the Hugo cohort had received MAPK inhibitors, which may mediate cross-resistance to ICBs (6). By comparing co-mut⁺ with PD-L1 expression and TMB in clinical NSCLC and melanoma cohorts, we found that co-mut⁺ might be a clinically practicable biomarker predictive of response to ICBs and survival. Currently, PD-L1 expression is a standard predictive biomarker for first-line pembrolizumab (3). To date, the treatment for patients with PD-L1 expression < 50% is platinum-based chemotherapy, either alone (32) or in combination with pembrolizumab (33). Biomarker selection for the use of ICBs in the population remains debatable. We have demonstrated that co-mut⁺ patients may benefit from anti-PD-1 therapy even if the PD-L1 expression is less than 50%. Similarly, in the TMB-low subset, the comutation status may enrich a patient population with a higher response rate to ICBs. Comutation is a qualitative biomarker that does not require a cut-off value. Standard next-generation sequencing (NGS) may be able to determine comutation status accurately. This method would be relatively simple and cost-effective compared with measurement of TMB. Prevalence of co-mut⁺ status in the most common cancer types is estimated at 10% to 20%, which is higher than the prevalence of dMMR/MSI-H (12). Furthermore, the gene panel for the determination of comutation status is detectable in plasma cell–free DNA (cfDNA). These features may make comutations in the DDR pathways an attractive biomarker for ICBs. Future prospective studies are warranted.

A recent study suggested that deleterious mutations in 34 DDR genes designed in a MSK-IMPACT NGS gene panel was associated with better clinical outcome of ICBs in urothelial cancer (16). However, we extracted the 34 DDR genes included in the MSK-IMPACT gene panel and failed to obtain positive results of these deleterious DDR genes in predicting a clinical benefit in NSCLC or melanoma treated with anti-PD-1 (Supplementary Fig. S8A and S8B). One important potential reason was that the TMB and the corresponding NAL caused by a single DDR mutation were not high enough to lead to significantly improved immunogenicity. According to our results, patients with mutations in a single DDR pathway did not show a superior response to ICBs compared with that of those with wild-type DDR. The median TMB levels of patients with single mutations in DDR genes were 52 (TCGA) and 53 (ICGC), which were very different from the previously reported predictive TMB cut-off values (6, 10, 21, 22). On the basis of these data, comutations in DDR pathways could better predict ICB efficacy than a single mutated DDR pathway could.

This study is first limited by the retrospective profiling of the analysis. However, TCGA is the most reliable public dataset, and validation using the ICGC database supports the reliability of our results. Moreover, we here did not differentiate whether the DDR gene mutations were functional. Our attempts to recruit functional DDR mutations into our co-mut⁺ pattern was handicapped by the limited information available regarding the functions of different mutations and the lack of hotspots in DDR gene mutations. In addition, some benign mutations may also influence protein function (34). Therefore, it is very difficult to clearly differentiate functional DDR mutations, and thus, missense DDR mutations were utilized in this study. Finally, we could not fully exclude the functions of comutations in other DDR pathways. However, the necessity of at least 2 pathways and the co-mut⁺ pattern provided optimal accuracy for TMB estimation and efficacy prediction for ICB delivery.

In conclusion, comutations in DDR pathways are associated with an immune phenotypic profile of increased TMB, NAL, and immune response–related mRNA expression. Preliminary data from four clinical cohorts strongly suggested better treatment outcomes of ICBs in co-mut⁺ patients. Comutations of HRR-MMR or HRR-BER are a potential predictive biomarker for ICB therapy that warrants future prospective investigation.

T.S.K. Mok. reports receiving commercial research grants from AstraZeneca, Boehringer Ingelheim, Bristol-Myers Squibb, Clovis Oncology, Merck Sharp & Dohme, Novartis, Pfizer, SFJ Pharmaceuticals, and XCovery; has ownership interest (including stocks and patents) in Sanomics and Hutchison Chi-Med; is a consultant/advisory board member for AstraZeneca, Boehringer Ingelheim, Roche/Genentech, Pfizer, Eli Lilly, Merck Serono, Merck Sharp & Dohme, Novartis, SFJ Pharmaceuticals, ACEA Biosciences, Vertex Pharmaceuticals, geneDecode, OncoGenex Technologies, Celgene, Ignyta, Fishawack Facilitate, Cirina, Janssen, Takeda, Hutchison Chi-med, OrigiMed, Hengrui Therapeutics, Sanofi-Aventis, Yuhan Corporation, and Amoy Diagnostics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: Z. Wang, S. Cai, T.S.K. Mok, J. Wang

Development of methodology: Z. Wang, S. Cai, T.S.K. Mok, J. Wang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Z. Wang, Z. Zhang, F. Zhang (Peking University), H. Dong, J. Duan, H. Bai, R. Wan, L. Liu, T.S.K. Mok, J. Wang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Z. Wang, J. Zhao, G. Wang, Y. Zhang, X. Zhao, J. Duan, Y. Tian, S. Wang, S. Cai, T.S.K. Mok, J. Wang

Writing, review, and/or revision of the manuscript: Z. Wang, J. Zhao, G. Wang, F. Zhang (Chinese PLA General Hospital), Y. Zhang, X. Zhao, H. Bai, M. Han, Y. Cao, L. Xiong, S. Wang, S. Cai, T.S.K. Mok, J. Wang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): H. Dong, X. Zhao, J. Duan, S. Wang, T.S.K. Mok, J. Wang

Study supervision: T.S.K. Mok, J. Wang

We thank the financial support from the National Natural Sciences Foundation Key Program, National Key R&D Program of China, National High Technology Research and Development Program 863, CAMS Innovation Fund for Medical Sciences, China National Natural Sciences Foundation, Beijing Natural Science Foundation, and Beijing Novel Program Grants. This work was supported by grants from the National Natural Sciences Foundation Key Program (81630071 and 81330062 to J Wang), National Key R&D Program of China (2016YFC0902300 to J Wang), National High Technology Research and Development Program 863 (SS2015AA020403 to J Wang), CAMS Innovation Fund for Medical Sciences (CIFMS 2016-I2M-3-008 to J Wang), China National Natural Sciences Foundation (81472206 to J Wang and 81871889 to Z.J. Wang), Beijing Natural Science Foundation (7172045 to Z.J. Wang), and Beijing Novel Program Grants (Z141107001814051 to J. Wang).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.