Geospatial Immune Heterogeneity Reflects the Diverse Tumor–Immune Interactions in Intrahepatic Cholangiocarcinoma

Intrahepatic cholangiocarcinoma (iCCA) is the second most common primary liver cancer, with increasing incidence worldwide (1). Characterized by high invasiveness and frequent postoperative recurrence, iCCA exhibits one of the highest mortality rates among human cancers. Despite the recent progress with novel therapies targeting IDH1 mutation (2) and FGFR2 fusion (3–5), iCCA is still incurable in most cases. Immune-checkpoint inhibitors (ICI) have revolutionized the standard treatment for patients with cancer, but the objective response rates range from only 3% to 22% in biliary cancers (6). One major challenge is the poor understanding of the tumor–immune interaction in iCCA, which impedes the identification of patients for effective immunotherapy.

The tumor–immune interaction is a dynamic and continuous process. Infiltrated immune cells hold the potential of killing malignant cells, representing the major external selection pressure to orchestrate tumor evolution (7). We have previously illustrated the spatiotemporal genomic evolution of iCCA through multiregion sampling–derived cancer cell cultures, linking with drug resistance, but we did not focus on the immune microenvironment (8). Single-cell RNA sequencing (RNA-seq) analyses have emphasized the impact of transcriptomic diversity of cancer cells on patient survival and treatment response in iCCA (9, 10). Accumulating studies have uncovered the effect of the tumor–immune interaction on tumor progression, drug resistance, and immune tolerance in many other cancers (11–13). Likewise, a recent study has also characterized the high intratumoral heterogeneity (ITH) and “cold” tumor microenvironment in patients with IDH-mutant iCCA (14). An in-depth exploration of such a relationship in iCCA may have the promise to enable more precise development for targeted therapy and immunotherapy (15).

Herein, we collected 207 tumor samples from 45 patients with treatment-naïve iCCA and performed whole-exome sequencing (WES), transcriptome sequencing, T-cell receptor sequencing (TCR-seq), and multiplex staining to portray the immuno­genomic landscape of iCCA. Our results demonstrate the diversity of immune infiltration in iCCA, and those with sparse, heterogeneous, or high immune infiltration have different genetic characteristics and escape mechanisms. The integration of immunogenomic analyses provides a novel insight into how iCCA's genetic makeup affects immune composition and function and vice versa.

To elucidate the immunogenomic landscape of iCCA, we prospectively assembled 207 tumor samples, together with paired nontumor liver tissues and peripheral blood, from 45 patients (four to six tumor regions per patient), and performed high-depth WES, RNA-seq, TCR-seq, and multiplex immunofluorescence (Supplementary Fig. S1A and S1B; Supplementary Table S1). None of the patients had received preoperative anticancer treatments. Two tumor samples were excluded for further analysis due to low tumor purity (<0.1). The average depth of WES was 284.7× for tumor samples (range, 196.5–584.7) and 247.5× for normal controls (range, 150.8–527.9), respectively. In total, we identified 30,348 somatic mutations (range, 47–509, averaging 148 per tumor sample), including 27,749 (91.4%) point mutations and 2,599 (8.6%) indels (Supplementary Table S1).

We observed considerable intratumoral genomic heterogeneity, with a median of 65.1% (range, 43.6%–89.3%) of single-nucleotide variations (SNV) and 66.0% (range, 9.1%–100%) of copy-number variations (CNV) identified as subclonal (Fig. 1A). The median pairwise SNV ITH and CNV ITH were 42.8% (range, 23.8%–71.4%) and 41.1% (range, 3.6%–81.8%), respectively, suggesting that over 40% of the SNVs and CNVs were distinct among tumor subregions within the same patients (Fig. 1A; Supplementary Table S1). There was a moderately positive correlation between SNV ITH and CNV ITH (Spearman correlation = 0.402, P = 0.006), indicating that tumor evolution coexisted at the mutational and chromosomal levels (Fig. 1B). Loss-of-function mutations in BAP1 can damage genomic stability and DNA repair (16), and iCCA patients with BAP1 mutations showed significantly higher levels of SNV ITH (Fig. 1C).

We identified a total of 510 driver mutations (range, 4–30 per patient, averaging 11.4), of which 237 (46.5%) were clonal and 273 (53.5%) were subclonal. We further investigated the timing of recurrent driver mutations by integrating three recent multiregion sampling studies of iCCA (refs. 8, 14, 17; Fig. 1D; Supplementary Fig. S1C). Generally, TP53, KRAS, BAP1, and PBRM1 mutations were clonal, highlighting their essential roles in tumorigenesis. IDH1/2, ARID1A, SMAD4, ATM, and PIK3CA mutations were primarily clonal but also subclonal in some patients, implying their dual roles in driving tumorigenesis and progression. We found that subclonal driver mutations could also induce pathway alterations of tumor cells in corresponding subregions (Supplementary Fig. S2). For example, P41–R3 with ARID1A mutation displayed downregulated DNA double-strand break repair compared with other subregions in P41, consistent with ARID1A's role in mismatch repair (18).

FGFR2 fusions were identified in four patients of our cohort, with the same breaking point but different partner genes (Supplementary Fig. S3A). The occurrence of the FGFR2 fusion was validated in all tumor subregions and nearly all tumor cells using Sanger sequencing and fluorescent in situ hybridization (FISH), respectively, suggesting that it occurred early during tumor evolution (Fig. 1E; Supplementary Fig. S3A). Meanwhile, clonal FGFR2 mutations were detected in another two patients, and the six cases with FGFR2 alterations exhibited significantly higher expression of FGFR2 than wild-type (Fig. 1F). We also obtained multiomics data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) iCCA cohort (n = 249; ref. 19) and analyzed the CNV of FGFR2. The frequency of FGFR2 loss was similar between our cohort and the CPTAC cohort (15.6% vs. 13.7%), whereas a lower frequency of FGFR2 gain was observed in ours (8.3% vs. 17.3%; Supplementary Fig. S3B). Focal FGFR2 amplification was only found in one patient in our cohort, co-occurred with FGFR2 fusion, and existed in all tumor subregions on FISH images (Supplementary Fig. S3C). There was no significant difference in FGFR2 expression between samples with FGFR2 copy-number neutral and gain in both cohorts, in contrast to its mutation/fusion-induced higher expression (Supplementary Fig. S3D).

Chromosome 6q loss (where HLA class I molecules are located) occurred early, implying an immediate decline of antigen presentation after tumor initiation (Fig. 1G). Loss of TP53 (17p) and CDKN2A (9p) was also an early event, whereas gain of EGFR (7p) and BRAF (7q) was a late event. We then analyzed mutation signatures and found different distribution between clonal and subclonal mutations (Supplementary Fig. S3E). Signature 16 (special signature in liver cancer) and signature 1 (spontaneous deamination of 5-methylcytosine correlating with age) were dominant in clonal mutations, whereas the signatures of subclonal mutations were more diverse, including signatures 1, 5, 8, 9, 15, and 25. Taken together, iCCA showed evolutionary branches at the early stage, with both driver gene alterations (TP53, KRAS, BAP1, FGFR2, PBRM1) and CNVs (6q loss, 8p loss, 17p loss, and 1q gain) in tumor initiation, together with persistent events contributing to the complex landscape of genomic ITH (Fig. 1H).

The clinical success of ivosidenib (2) and pemigatinib (3) may be attributed to the clonality of IDH1 mutation and FGFR2 fusion, respectively (Fig. 1A). Considering the prevalence of driver mutations in three other iCCA cohorts (19–21), we identified KRAS, PIK3CA, ARID1A, ATM, BRCA1, and BRCA2 mutations as the most promising therapeutic targets (Fig. 1I). KRAS mutations occurred in 7.5% to 26.7% of patients with iCCA, presenting as hotspot mutations such as KRAS^G12D and KRAS^G12V. Recently, multiple small molecules, including AMG510, MRTX849, and JNJ-74699157, have been developed to target KRAS^G12C specifically (22–24). However, KRAS^G12D, instead of KRAS^G12C, was the dominant mutant in iCCA but still lacked specific inhibitors. Gain-of-function mutations in PIK3CA occurred in 3.2% to 8.9% of patients, and those patients were ideal candidates for alpelisib treatment (25). Loss-of-function alterations in BRCA1, BRCA2, ATM, and ARID1A, occurring in 15.0% to 37.8% of patients, damaged homologous recombination repair and conferred sensitivity to PARP inhibitors such as olaparib (26).

To explore the influence of genomic ITH on the immune microenvironment, we used the immune signature by Danaher and colleagues (27) to characterize the immune context (Supplementary Fig. S4A and S4B; Supplementary Table S2). At the tumor region level, unsupervised hierarchical clustering stratified 205 tumor regions into two distinct clusters: high (85 regions from 33 patients) and low (120 regions from 38 patients) immune infiltration (Fig. 2A). Nearly all immune cell subgroups exhibited positive correlations, especially for lymphocyte subsets (Supplementary Fig. S4C). Strong positive correlations of transcriptomic immune deconvolution with multiplexed immunostaining results and histologic tumor-infiltrating lymphocyte (TIL) estimates were observed (Supplementary Fig. S4D and S4E). As expected, the high-infiltration cluster showed a significantly higher level of histologic TILs, intense immunostaining densities of CD3⁺ T cells, CD8⁺ T cells, and CD68⁺ macrophages (P < 0.001; Supplementary Fig. S4F and S4G), and elevated cytolytic activity (CYT; ref. 28) and IFNγ signature (29), demonstrating stronger immune responses in these regions (P < 0.001; Supplementary Fig. S4H). Intriguingly, the high-infiltration cluster displayed comparable tumor mutation burden (TMB) but significantly higher tumor neoantigen burden (TNB), implying that the immunogenicity of mutations might direct immune infiltration (Supplementary Fig. S4I).

At the patient level, the 45 patients were stratified into three immune groups: sparsely infiltrated (uniformly low immune infiltration across subregions, 12 patients containing 52 regions), heterogeneously infiltrated (heterogeneous immune infiltration among subregions, 26 patients containing 120 regions), and highly infiltrated (uniformly high immune infiltration across subregions, 7 patients containing 33 regions; Fig. 2A). Heterogeneous infiltration among different tumor regions was further confirmed by multiplex immunostaining and histology (Fig. 2B; Supplementary Fig. S4J). We applied four published transcriptome-based immunophenotypic classifiers (30–33) to stratify our multiregion samples, revealing that the heterogeneously infiltrated group harbored higher proportions of tumor subregions distributed in different immunophenotypes within the same patients (Fig. 2C; Supplementary Fig. S5A). Likewise, the IFNγ signature was also more variant in this group (Supplementary Fig. S5B). The results highlighted the challenges in the immune stratification and immunotherapy prediction of iCCA patients with a single biopsy.

As evaluated by multiplexed immunostaining on tumor microarrays (TMA; including matched tumor subregions, tumor margin, and nontumor liver; Supplementary Fig. S1B), the three immune groups displayed comparable densities of CD3⁺ T cells, CD8⁺ T cells, and CD68⁺ macrophages in tumor margins and nontumor livers but significantly differed in tumor subregions (Fig. 2D). T-cell exclusion was sequentially aggravated from highly infiltrated to sparsely infiltrated tumors (Fig. 2E). In contrast, macrophages were more abundant in tumor regions than in surrounding tissues (Supplementary Fig. S5C), suggesting that T cells were the targets of exclusion. Immunogenomic analysis showed that highly infiltrated tumors displayed comparable cancer testis antigen (CTA) load, chromosomal instability, and ploidy with heterogeneously and sparsely infiltrated tumors but had significantly higher TNB and the lowest tumor purity (Fig. 2F). Sparsely infiltrated tumors had comparable TMB and TNB with heterogeneously infiltrated tumors but exhibited universally fewer immune and stromal cells. Meanwhile, multiple chemokines, including CXCL9, CXCL10, and CXCL11, were inadequate in sparsely infiltrated tumors compared with highly infiltrated tumors, and heterogeneously infiltrated tumors displayed the distinct distribution of chemokines between tumor subregions with low and high immune infiltration (Fig. 2F; Supplementary Fig. S5D). We then compared transcriptomic features among the three immune groups using hallmarks from the Molecular Signatures Database (MSigDB; ref. 34). Immune-related pathways, angiogenesis, and classical cancer pathways were upregulated in highly infiltrated tumors, while sparsely infiltrated tumors upregulated cell-cycle signaling and several metabolic pathways such as oxidative phosphorylation and cholesterol homeostasis (Supplementary Fig. S5E). These results revealed distinct immune landscape, genomic, and transcriptomic characteristics among the three immune groups. We further calculated an immune ITH index for each patient and observed significant positive correlations between ITH of genomic features and immune contexts, indicating their potential interaction (Supplementary Fig. S5F).

A specific genetic profile could affect the infiltration and composition of immune cells in the tumor microenvironment (35). We combined the multiomics data of our cohort and the CPTAC iCCA cohort to investigate the correlations between recurrent gene alterations and immune groups (Fig. 3A; Supplementary Table S3). The result showed that iCCA with KRAS mutations had a tendency for more myeloid infiltration but fewer lymphocytes (Fig. 3A), reflecting its role in reprogramming the local milieu to support tumor progression (36, 37). Contrarily, negative correlations were found between BAP1 or ARID1A mutations and neutrophil infiltration (Fig. 3A). Recent studies reported that IDH1/2 mutations could shape a “cold” immune microenvironment in iCCA (14, 38), but patients with IDH1/2 mutations displayed more T-cell infiltration in our cohort. We found that three of the eight cases with IDH1/2 mutations were highly infiltrated tumors, with two of them having hepatitis B virus (HBV) infection and the rest showing a relatively high TNB (ranking 2/45 in our cohort). HBV infection positively correlated with lymphocyte infiltration in iCCA, similar to the finding in hepatocellular carcinoma (39), irrespective of IDH1/2 mutation or wild-type (Supplementary Fig. S6A and S6B). Thus, the immune impact of IDH1/2 mutations could be influenced by coexisting factors such as hepatitis and neoantigenicity, which may explain the unexpected genotype–immune correlations in our cohort. Indeed, in patients with HBV infection (n = 69) and liver cirrhosis (n = 22) from the CPTAC cohort, IDH1/2 mutations also showed a trend of positive correlations with immune infiltration (Supplementary Fig. S6C).

FGFR2 fusion is a fascinating drug target in iCCA (6), but its influence on the immune microenvironment remains unclear. Notably, FGFR2 alterations (fusion and mutation) were consistently associated with fewer immune cells, including total TILs, CD4⁺ T cells, and macrophages, together with lower TMB and TNB, in both our cohort and the CPTAC cohort (Fig. 3A and B). Downregulation of immune-related pathways and genes of cytotoxicity, chemotaxis, immune exhaustion, and inhibition was noted in iCCA with FGFR2 alterations (Supplementary Fig. S7A and S7B). Immunostaining on the TMA (n = 177) of the CPTAC cohort verified that iCCA with FGFR2 alterations had fewer CD3⁺ T cells (P = 0.033), CD20⁺ B cells (P = 0.003), and CD15⁺ neutrophils (P = 0.030) than wild-type tumors (Fig. 3C; Supplementary Fig. S7C). Furthermore, we constructed murine iCCA models with or without an FGFR2::BICC1 fusion by a hydrodynamic tail-vein injection. KRAS/p19/FGFR2 (KPF) and KRAS/p19 (KP) models led to macroscopic isolated iCCA tumors, whereas AKT/YAP/FGFR2 (AYF) and AKT/YAP (AY) models resulted in diffuse iCCA nodules in the liver (Fig. 3D; Supplementary Fig. S7D). KPF and AYF tumors had fewer CD4⁺ T cells, CD8⁺ T cells, and F4/80⁺ macrophages than KP and AY tumors, respectively, highlighting the repressed immune microenvironment by the FGFR2::BICC1 fusion protein (Fig. 3E; Supplementary Fig. S7E). RNA-seq analysis identified 3,055 differentially expressed genes between KPF and KP tumors, including downregulation of Cxcl2, Ccl6, Cd163, Cd24a, Cd177, and Ctla4 in KPF tumors (Fig. 3F). Downregulated genes in KPF tumors were enriched in immune responses like Th1/Th2 differentiation, TCR signaling, and immune checkpoint pathways, whereas upregulated genes were involved in pathways of multiple metabolic processes (Fig. 3G), consistent with the above findings in human samples.

Neoantigen depletion is a common immune escape mechanism of tumors (40). Like mutations, neoantigens also showed significant ITH where only a median of 36.1% (range, 9%–60%) were identified as clonal (Supplementary Fig. S8A; Supplementary Table S4). Highly infiltrated tumors tended to have more clonal and subclonal neoantigens, but the ratios of clonal neoantigens were comparable among the three immune groups (Supplementary Fig. S8B and S8C). Considering the complex process from gene mutations to neoantigen presentation, we only analyzed clonal neoantigen depletion at DNA (through copy-number loss), mRNA (through expression silence), and presentation levels [through loss of heterozygosity (LOH) of HLA]. At least half of clonal neoantigens could not be presented, where copy-number loss, expression silence, and HLA LOH led to the depletion of 2.2%, 47.1%, and 3.8% of clonal neoantigens, respectively (Fig. 4A).

Copy-number loss of clonal neoantigens occurred in 17.6% (36/205) of tumor regions, leading to 2.3% (46/2,034), 2.9% (152/5,164), and 0.9% (29/3,239) of clonal neoantigen depletion in sparsely, heterogeneously, and highly infiltrated tumors, respectively (Supplementary Fig. S8D). We compared clonal neoantigens with nonneoantigenic, nonsynonymous clonal mutations to determine whether clonal neoantigen loss is an active event as previously described (40). Clonal neoantigen loss occurred more frequently than their nonneoantigenic counterparts in sparsely (P = 0.002) and heterogeneously infiltrated tumors (P = 0.046), indicating active clonal neoantigen loss in those tumors but not in highly infiltrated tumors (P = 1.00; Fig. 4B). Further dissection of heterogeneously infiltrated tumors showed that tumor regions with low infiltration, instead of high infiltration, actively abandoned clonal neoantigens through copy-number loss, authenticating the above assumption (Supplementary Fig. S8E). We then quantified the DNA immunoediting score in each tumor region, where the lower score indicates stronger DNA immunoediting (28). Intriguingly, tumor regions with clonal neoantigen loss exhibited weaker DNA immunoediting than those without clonal neoantigen loss in sparsely and heterogeneously infiltrated tumors (Fig. 4C), suggesting that DNA immunoediting was a supplement to clonal neoantigen loss.

Mutated gene silence could be directed by mutation itself or transcriptional regulation. A median of only 52.0% (range, 12.5%–100%) of clonal neoantigens were expressed in iCCA, and the proportion was significantly lower in heterogeneously infiltrated tumors (median, 48.8%) than in sparsely infiltrated tumors (median, 55.6%; P < 0.001; Fig. 4D), partly reflecting the selection of neoantigen expression under immune pressure. However, this proportion rebounded to 53.8% in highly infiltrated tumors, possibly due to the extensive defects in antigen presentation machinery and T-cell exhaustion in this group (described below). Alternatively, other potential mechanisms, such as epigenetic silence by DNA methylation (41), may also exist, which is beyond the scope of this study. We indeed observed a higher frequency of HLA LOH in highly infiltrated tumors (71.4%, 5/7) than heterogeneously (34.6%, 9/26) or sparsely infiltrated (33.3%, 4/12; Fig. 4E) tumors. For highly infiltrated tumors, clonal neoantigens tended to be expressed more frequently than their nonneoantigenic counterparts in tumor regions with HLA LOH (P < 0.001; Fig. 4F), suggesting that HLA LOH might liberate clonal neoantigen from expression silence. Meanwhile, increased copy-number loss of B2M, TAP1/TAP2, and HSPA1A/HSPA1B in highly infiltrated tumors further damaged antigen presentation to avoid expression silence of clonal neoantigens compared with heterogeneously infiltrated tumors (Supplementary Fig. S8F and S8G).

Parallel evolution is a common phenomenon as a result of positive selection in tumors (42). We observed parallel clonal neoantigen depletions through distinct mechanisms in 17 of 45 tumors, more so in highly infiltrated tumors (5/7, 71.4%) than in sparsely (3/12, 25.0%) and heterogeneously infiltrated (9/26, 34.6%; P = 0.086) tumors. For example, the copy-number loss of ZFP90^C145T occurred in P38-R3, whereas the subclonal loss of HLA-B 39:01 damaged the presentation of the neopeptide ENYSYLVSL derived from ZFP90^C145T in P38-R2 and R4, resulting in the depletion of ZFP90^C145T-derived neoantigens in P38-R2, R3, and R4 (a highly infiltrated tumor; Fig. 4G).

We also compared the clonal neoantigen depletion derived from driver and nondriver mutations. Copy-number loss of clonal nondriver mutation–derived neoantigens (225/9,872, 2.3%) occurred more frequently than those from driver mutations (2/565, 0.4%; Supplementary Fig. S8H). As such, significantly higher proportions of clonal driver mutation–derived neoantigens (499/565, 88.3%) were expressed than nondriver mutation–derived neoantigens (4,874/9,872, 49.4%) in all three immune groups (Fig. 4H). In terms of immunogenicity, driver mutations were predicted to produce significantly fewer and also less immunogenic neoantigens than nondriver mutations, which may represent innate immune evasion (Supplementary Fig. S8I).

The assessment of subclonal neoantigen depletion may be complicated by sequencing depth and tumor purity. Alternatively, we calculated subclonal TNB for each tumor, defined as the sum of subclonal neoantigen numbers multiplied by corresponding cancer cell fractions (Supplementary Fig. S8J). Overall, within each patient, tumor regions with higher TIL scores exhibited significantly lower subclonal TNB compared with other regions (P = 0.009, linear mixed model). Such negative correlation was mainly contributed by heterogeneously infiltrated tumors (P = 0.007; Fig. 4I), as it disappeared in highly (P = 0.114) and sparsely infiltrated tumors (P = 0.555). To dissect the role of tumor purity, we calculated the subclonal burden of nonneoantigenic, nonsynonymous mutations using the same method. Indeed, subclonal TNB, rather than subclonal nonneoantigenic mutation burden, was significantly lower in regions with high versus low infiltration within heterogeneously infiltrated tumors (Fig. 4J). These findings highlighted the influence of heterogeneous immune infiltration on iCCA subclonal evolution.

To further evaluate T-cell quality, we performed TCR-seq on all tumor samples and peripheral blood (one patient using nontumor liver tissue) and obtained an average of 18,776 (range, 52–115,446) unique TCRβ rearrangements per tumor region. Compared with sparsely and heterogeneously infiltrated tumors, highly infiltrated tumors displayed more unique TCR clones but comparable TCR clonality (Supplementary Fig. S9A and S9B). The peripheral blood of patients with highly and heterogeneously infiltrated tumors displayed more unique TCR clones but also comparable TCR clonality compared with the peripheral blood of patients with sparsely infiltrated tumors (Supplementary Fig. S9A and S9B). Then, we investigated the overlap of TCR clones across different tumor regions. Consistent with previous studies (43, 44), most TCR clones were unique to individual tumor regions, ranging from 38.7% to 79.9%, whereas there were also TCR clones ubiquitous among all regions of the tumor, with as little as 0.9% and at most 30.6% (Fig. 5A). The percentages of shared TCR clones across tumor regions and between tumor and peripheral blood were both significantly higher in highly infiltrated tumors, indicating that extensive intratumoral immune responses could also cause perturbations in the peripheral TCR repertoire (Supplementary Fig. S9C). The Morisita–Horn index (44) was calculated to measure the similarity of the TCR repertoires between two tumor regions within a tumor, and highly infiltrated tumors exhibited more similar TCR profiles across regions as expected (Supplementary Fig. S9D).

Of particular interest was P25 (a highly infiltrated tumor), who exhibited clonal copy-number loss of B2M and HLA-B 58:01, along with wide parallel subclonal mutations in the remaining B2M allele, including B2M p.L13 frameshift deletions in R1, R3, and R5, and a B2M 72–78 nonframeshift deletion in R2 (Supplementary Fig. S9E). Within P25, R4 exhibited the strongest subclonal immunoediting with the lowest immune infiltration and TCR clonality, whereas the other four regions developed B2M mutations to further destroy antigen presentation (Supplementary Fig. S9E). Interestingly, TRBV4-1-TRBJ1-4 (CDR3aa: CASSQGQGREKLFF) and TRBV15-TRBJ1-2 (CDR3aa: CATSSGLSGADYGYTF) were dominant in P25-R1, R3, and R5, whereas P25-R2 exhibited a distinct TCR repertoire, which might be the result of different B2M mutation patterns and a subclonal neoantigen profile, implying the interaction between tumor evolution and TCR composition.

Multiple immune checkpoints, including PD-1, PD-L1, and CTLA4, were serially overexpressed from sparsely infiltrated tumors to highly infiltrated tumors (Supplementary Fig. S9F). Sparsely infiltrated tumors exhibited comparable proportions of PD-1⁺CD8⁺ T cells with corresponding tumor margins and nontumor livers, whereas highly and heterogeneously infiltrated tumors showed significantly increased proportions over paired tumor margins and nontumor livers, indicating that CD8⁺ T cells were gradually programmed into an exhaustion phenotype (Fig. 5B; Supplementary Fig. S9G). Also, a higher proportion of PD-L1⁺ stromal cells was observed in highly infiltrated tumors (Fig. 5B; Supplementary Fig. S9H), and CD8⁺ T-cell exhaustion scores sequentially increased from sparsely to heterogeneously to highly infiltrated tumors (Supplementary Fig. S9I). In addition, highly and heterogeneously infiltrated tumors displayed more CD4⁺FOXP3⁺ regulatory T-cell (Treg) and CD206⁺ M2-like macrophage infiltration than sparsely infiltrated tumors (Fig. 5C; Supplementary Fig. S9J and S9K). All three groups exhibited a significantly higher proportion of Tregs and a lower ratio of CD8⁺ T cells/Tregs in tumor regions than in corresponding tumor margins and nontumor livers, highlighting the widespread immunosuppressive role of Treg accumulation in iCCA (Fig. 5C; Supplementary Fig. S9L).

We then integrated immune infiltration status and immune evasion mechanisms to predict patient prognosis, considering that immune diversity (i.e., sparse, hetero, and high immune groups) itself had minimal prognostic significance. Patients were classified as the low evasion group (highly or heterogeneously infiltrated with no evidence of HLA LOH or extreme CD8⁺ T-cell exhaustion) and the high evasion group (sparsely infiltrated or exhibiting evidence of HLA LOH or extreme CD8⁺ T-cell exhaustion; Fig. 5D). In our cohort, the high evasion group had significantly shorter recurrence-free survival than the low evasion group (HR = 2.24, P = 0.040; Fig. 5E). We also divided the CPTAC iCCA cohort into low and high evasion groups, and patients with high evasion ability had shorter overall survival (HR = 1.90, P = 0.011; Fig. 5F). Notably, immune evasion ability remained as an independent prognosticator in multivariable analysis in both cohorts (Supplementary Table S5). Again, a similar result was found in The Cancer Genome Atlas (TCGA) iCCA cohort (n = 32, HR = 3.62, P = 0.046), corroborating the robustness of the prognostic classification system (Fig. 5G).

Intrahepatic metastasis (IM) is a common metastatic route of iCCA. We plotted phylogenetic trees and CNV profiles of four tumors with IMs and inferred the most likely tumor region origin for each metastasis (Fig. 6A and B). For example, P11-IM2 and P11-R2 had the closest clonal evolutionary trajectories and CNV profiles, indicating that P11-IM2 was disseminated from tumor clones in P11-R2 (Fig. 6A). There was no special enrichment of driver mutations or CNVs in metastases, and metastases could derive from tumor clones of different regions, suggesting that the whole tumor, rather than special clones, had the potential to disseminate.

Then, we performed unsupervised clustering based on the top 3,000 variant genes in each of the four patients (Supplementary Fig. S10A). There was no obvious separation between primary tumors and IMs, and the metastases even clustered together with the original tumor regions in 75.0% (6/8) of pairs. Unsupervised clustering using immune signatures of Danaher and colleagues (27) for each patient (Fig. 6C) showed that the IMs also clustered with original tumor regions in 75.0% (6/8) of pairs. Further immunostaining revealed that IMs exhibited similar infiltration of T cells and macrophages with their original tumor regions in P13 and P22 (Fig. 6D; Supplementary Fig. S10B). Overall, the metastatic clones with similar genomic and transcriptomic characteristics to primary clones within the same organ can shape a homologous immune microenvironment. Otherwise, distant metastases out of the original organ may generate distinct transcriptomic and immune profiles, despite stable genomic traits (45). Likewise, the TCR repertoire and high-frequency TCR clones were similar between paired primary tumors and IMs (Supplementary Fig. S10C and S10D), possibly due to the fact that they shared genomic features and consequently neoantigen profiles.

Neoantigens derived from driver mutations, such as TP53^R175H and KRAS^G12D, could activate a T-cell response, and neoantigen-specific T-cell clones could be used for adoptive cell therapy (46–48). Here, a total of 32 predicted neopeptides from driver mutations of TP53, KRAS, IDH1, IDH2, and BAP1 were synthesized to systematically estimate their potential as immunotherapeutic targets (Fig. 7A; Supplementary Fig. S1C; Supplementary Table S6). We used tetramer exchange experiments to evaluate the affinity of mutant peptides with HLA-A 02:01, HLA-A 24:02, and HLA-A 11:01 (the three most common HLA-A subtypes). The binding of peptides to HLA-A molecules was specific, where the same peptide could have completely different affinities with different HLA-A molecules (Fig. 7B; Supplementary Fig. S11A). The exchange efficiency was as high as 97.2% between peptide #3 and HLA-A 02:01, but as low as 7.8% between peptide #3 and HLA-A 11:01 (Fig. 7B), with 15.7% (13/83) of peptide–major histocompatibility complex (MHC) pairs defined as high affinity (exchange efficiency >75%; Fig. 7C).

In addition to the affinity with HLA molecules, the immunogenicity of peptides also determines their ability to stimulate immune cells. We isolated peripheral blood mono­nuclear cells (PBMC) from three healthy donors and performed an enzyme-linked immunospot (ELISpot) assay to quantify IFNγ secretion against different peptides (Fig. 7D). Strong reactivity against peptide #19 (LVVVGADGV, KRAS^G12D), peptide #29 (SQEQPRCHY, IDH2^R37C), and peptide #30 (RLFERDGLKV, BAP1^L183R) was detected in donor #1, whereas donors #2 and #3 showed only moderate reactivity against peptide #19 (Fig. 7D). Because dendritic cells are powerful antigen-presenting cells to prime T-cell responses (49), we further isolated peripheral monocytes from donor #4 and loaded them with antigen peptides (peptides #19, 20, 29, and 30 and the combined peptide pool) during the differentiation into dendritic cells. Paired T cells were then cocultured with mature dendritic cells. Subsequent flow cytometry analysis of T-cell activation markers 4-1BB and CD107a demonstrated the enhanced reactivity of both CD4⁺ and CD8⁺ T cells against peptide #19 (LVVVGADGV, KRAS^G12D), which was superior to other peptides and comparable with the combined peptide pool (Fig. 7E; Supplementary Fig. S11B).

To identify peptide #19–specific TCR repertoire, we isolated activated 4-1BB⁺ cells stimulated by peptide #19 or DMSO (control) from donors #1 to 3 and performed TCR sequencing. TRBV9-TRBJ2-1 and TRBV25-1-TRBJ2-1 rearrangements were significantly enriched in the top 100 CDR3s whose frequencies varied mostly between peptide #19 and DMSO groups (Fig. 7F). We then used the GLIPH algorithm (50) to cluster potential peptide #19–specific activated CDR3 library (Supplementary Fig. S11C), and most of the TCRβ sequences fell in one or few related TCR repertoires, among which the motifs RAGQ and KVD were identified (Fig. 7G). This preliminary screening indicated that LVVVGADGV (KRAS^G12D) had strong immunogenicity that can induce specific V–J rearrangements and TCRβ motifs, which could be further exploited for personalized immunotherapy.

Through the integration of genomic, transcriptomic, and immunophenotypic data, our study elucidates the dynamic interactions between iCCA evolution and the immune microenvironment. We observed intratumoral heterogeneous immune infiltration in more than half of iCCA, posing a great challenge for the immune classification of patients and the prediction of immunotherapeutic efficacy by single biopsy. Tumors were classified into three immune groups with distinct immune microenvironment characteristics and neoantigen depletion mechanisms (Supplementary Fig. S11D), providing a novel insight to guide future immuno-oncologic therapeutic strategies.

Despite extensive genomic ITH in iCCA, some recurrent gene alterations, such as TP53, KRAS, FGFR2, BAP1, PIK3CA, ATM, and IDH1/2, were primarily clonal, providing considerable reliability to assess these mutations from a single biopsy and highlighting their potential as therapeutic targets. Among them, we found that FGFR2-altered iCCA displayed higher FGFR2 expression and less immune infiltration, together with low TMB and TNB. Recently, Palakurthi and colleagues (51) combined an FGFR inhibitor (erdafitinib) and anti–PD-1 therapy to treat an FGFR2-mutated lung cancer mouse model and induced significant tumor regression, with increased cytotoxic T-cell infiltration and decreased Tregs, macrophages, and PD-L1 expression. The finding raises the possibility that FGFR2 inhibitors may drive the immunologic alterations and sensitize FGFR2-altered iCCA to anti–PD-1 therapy. Surprisingly, we revealed that IDH1/2 mutations may not uniformly lead to sparse immune infiltration, but rather may be reversed by hepatitis virus infection and comutation profiles. This is analogous to Epstein–Barr virus–associated iCCA, which had abundant lymphocyte infiltration and was sensitive to ICIs compared with non–virus-associated iCCA (52).

A tumor must have developed immune evasion mechanisms to survive and grow in vivo. Neoantigen depletion could reduce the immunogenicity of tumor cells and serve as a major escape mechanism, which was observed during tumor recurrence (53), metastasis (54), and immunotherapy (55). Consistent with the finding in lung cancer (40), active clonal neoantigen loss mainly occurred in low-infiltration tumor regions in iCCA. We further found that DNA immunoediting, which reflects the reprogramming of somatic mutations to produce fewer neoantigens (28), was complementary to clonal neoantigen loss in low-infiltration tumor regions. Both of them may represent tumor-intrinsic ways at the genomic level to immune escape. Previous studies in colorectal cancer and hepatocellular carcinoma have shown that DNA immunoediting was prevalent in tumors with abundant T-cell infiltration but always lacking in those with HLA LOH (44, 54). These phenomena collectively suggest that different immune escape mechanisms may be dominant in distinct contexts, and pervasive antigen presentation disruption may explain the relatively low level of DNA immunoediting in highly infiltrated iCCA in our study.

At the mRNA level, expression silence of clonal neoantigens was prevalent, especially in heterogeneously infiltrated tumors. Of note, highly infiltrated tumors displayed a comparable expression ratio of clonal neoantigens with sparsely infiltrated tumors. One hypothesis is that frequent antigen presentation defects may relieve the expression silence, because highly infiltrated tumors with HLA LOH exhibited a significant increase of expressed clonal neoantigens compared with their nonneoantigenic counterparts. We assume that highly infiltrated iCCA resemble microsatellite instability–high (MSI-high) colorectal cancers in the following aspects: (i) abundant immune infiltration, (ii) high neoantigen burden, (iii) frequent disruption to antigen presentation, and (iv) overexpression of immune checkpoints (56, 57). These similarities suggest that highly infiltrated iCCA could be suitable candidates for immune-checkpoint blockade. Considering that DNA hypermethylation is one of the important causes of neoantigen expression silence in iCCA (17, 58), heterogeneously infiltrated tumors may be optimal for the combination of epigenetic regulators and anti–PD-1 (Supplementary Fig. S11D). These conclusions are based on an in-depth analysis of the multiomics data, and further functional validation would solidify the findings.

Extrinsic immune escape mechanisms also played a critical role in tumor immune tolerance (59). T-cell exclusion was extensive in iCCA, engendering the formation of immune-privileged microenvironments, especially in sparsely infiltrated tumors. Imbalanced CD8⁺ T cell/Treg was an essential characteristic of iCCA regardless of the abundance of immune infiltration as previously reported (60). Unsurprisingly, high immune infiltration was also accompanied by immune-checkpoint overexpression and CD8⁺ T-cell exhaustion, explaining why tumors can still exist in “hot” microenvironments (61, 62). Meanwhile, abundant PD-L1⁺ myeloid cells and stromal cells such as cancer-associated fibroblasts increased in highly infiltrated tumors, which may further weaken antitumor–immune response (63–65). These immunosuppressive cells help promote the cultivation of a high-infiltrated but useless or even tumor-promoting microenvironment. As a result, immune infiltration alone may not be sufficient to distinguish the outcomes of patients with iCCA, and thus we developed a simple but robust prognostic classification system based on immune evasion ability.

Although driver mutations produced fewer and less immunogenic neoantigens, they underwent weak immunoediting compared with passenger mutations, which indicates that driver mutation expression is required for tumor cell fitness but also exposes vulnerability to targeted or immune therapy. Therefore, driver mutation–derived neoantigens were better targets than passenger mutation–derived neoantigens because the latter can be easily edited through neoantigen depletion. In this regard, we screened a panel of recurrent mutation–derived neoantigens and demonstrated strong immunogenicity and high affinity to HLA-A 02:01 of LVVVGADGV (KRAS^G12D). Recently, Kadam and colleagues (66, 67) combined anti–PD-1 and vaccination of LVVVGADGV peptide to treat a KRAS^G12D lung cancer mouse model and induced 80% tumor eradication, further supporting it as an ideal target for personalized immunotherapy (68). Specific V–J rearrangements and TCRβ motifs were also identified by us, which may help design engineered TCR T cells (69) for patients with KRAS^G12D in the future.

In conclusion, our study comprehensively illustrated tumor–immune interactions by integrating multiomics in iCCA, demonstrating the spatiotemporally heterogeneous immunogenomic features across patients and within tumors. Future studies of paired pre- and posttreatment samples using high-resolution spatial omics and single-cell profiling may provide a greater understanding of these processes. Based on current and ongoing knowledge, new combinational immunotherapies may ultimately renovate the treatment paradigm and improve clinical outcomes of this fatal malignancy.

We prospectively collected 207 tumor samples from 45 patients with iCCA who underwent curative resection from January 2019 to December 2019 at the Zhongshan Hospital of Fudan University. A cut fully bisecting the resected tumor was made, and representative spatially separated regions were collected at least 5 mm away from each other and tumor margin. Tumor margins (defined as the region within 500 μm on each side of the border; ref. 70) and nontumor liver tissues (>2 cm away from tumor margin) were sampled from formalin-fixed, paraffin-embedded tissues to construct a TMA. This study was approved by the Research Ethics Committee of Zhongshan Hospital (B2017-060R) with written informed consent from each patient and was conducted in accordance with ethical guidelines (Declaration of Helsinki).

WES and RNA-seq were performed by the DNBSEQ Platform (MGI Tech Co., Ltd.). Somatic mutations of tumor samples were detected using the blood samples or nontumor liver tissues as controls, and raw RNA data were processed as previously described (44). Detailed tissue preparation, data processing, and somatic mutation detection are described in Supplementary Materials and Methods.

A total of nine published in silico immune tools (27, 28, 71–77) were collected and tested in our data (Supplementary Fig. S4A and S4B; Supplementary Table S2) as previously described (40). The proportion of genes in each method that met the criteria was calculated and further compared with multiplex immunofluorescence. Due to the lack of a CD4⁺ T-cell estimate in Danaher and colleagues (27), we used the estimate of CD4⁺ T cells from Shen and colleagues (75). Pathologic TILs were estimated using internationally established guidelines (www.tilsincancer.org; ref. 78).

Multiplexed immunofluorescence was performed as we previously described (79). A total of four staining regimens were performed, and detailed information on antibodies is listed in Supplementary Table S2. Multispectral images of TMAs were scanned using Vectra 3.0 (PerkinElmer), and spectral unmixing was conducted using the inForm Advanced Image Analysis Software (v2.3 PerkinElmer) with built spectral libraries. Cutoff value for positivity was determined according to the staining patterns and intensities, and cell frequency was counted using R script.

Six-week-old, female FVB/N mice were ordered from the Shanghai Branch of Beijing Vital River Laboratory Animal Technologies Co. Ltd. The experiments were performed following the institutional guidelines strictly and were approved by the Institutional Animal Care and Use Committee (IACUC) of the Shanghai Branch of Beijing Vital River Laboratory Animal Technologies Co. Ltd. (2017-0014). pT3-EF1α-HA-myr-Akt (mouse, Addgene, 31789, RRID: Addgene_31789) and pT3-EF1a-YAP^S127A (human, Addgene, 86497, RRID: Addgene_86497) were obtained from Addgene (https://www.addgene.org/). The FGFR2::BICC1 fusion gene and KRAS^G12D gene were synthesized and cloned into the pT3 vector, and the CRISPR/Cas9 single-guide RNA (sgRNA)-p19 was designed, synthesized, and cloned to the pX330 vector (Shanghai Xitubio Biotechnology Co., Ltd.). For the construction of AKT/YAP, AKT/YAP/FGFR2, KRAS/p19, and KRAS/p19/FGFR2 murine iCCA models, plasmids were injected into mice by the hydrodynamic tail-vein injection at a 10:1 ratio of transposon to SB-luc transposase-encoding plasmid (20 μg of pT3-EF1α-HA-myr-Akt plasmid, 30 μg of pT3-EF1a-YAP^S127A plasmid, 10 μg of pT3-EF1α-HA-FGFR2 fusion plasmid, 25 μg of pT3-EF1a-KRAS^G12D plasmid) as well as 10 μg of the PX330-CRISPR/Cas9 sgRNA-p19 plasmid with the solution into the tail vein with a total volume corresponding to 10% of the body weight in 6 to 8 seconds. Vectors for hydrodynamic delivery were produced using the Qiagen Plasmid Plus Mega Kit. Equivalent DNA concentration between different batches of DNA was confirmed to ensure reproducibility among experiments. Mice were sacrificed 6 weeks after injection.

IHC and quantification of immune cells were performed as described previously (80). Primary antibodies used in the IHC staining of mouse liver tumors are listed in Supplementary Table S2. Five microscopic fields (magnification, × 200) of regions of interest were selected and captured. Positively stained cells were counted using Image-Pro Plus (81), and the density of positive cells was calculated by averaging.

Neoantigens were predicted by NetMHC (v4.0; ref. 82) and NetMHCpan (v4.1; ref. 83) with the following criteria: (i) the length of the peptide was 9 to 10 mer and (ii) minimum IC₅₀ <500 nmol/L in either software. Copy-number loss of clonal neoantigen was detected if the corresponding copy number of the SNV allele was 0. If the read count of the mutation site identified by RNA-seq was more than 5 with at least 3 mutated reads, the neoantigen was defined as expressed neoantigen and otherwise unexpressed as previously described (40). Neoantigens bound to the lost HLA allele were defined as loss of presentation.

Peptide–MHC tetramer assays were performed to estimate the affinity between peptide and HLA sites as we previously described (44). Neoantigen peptides were synthesized according to the standard procedure (ABclonal Inc.) with purity >95%. An HLA-A*11:01 tetramer kit (TB-7304-K1, MBL International), an HLA-A*24:02 tet­ramer kit (TB-7302-K1, MBL International), and an HLA-A*02:01 tetramer kit (TB-7300-K1, MBL International) were used according to the manufacturer's instructions.

PBMCs of donors were obtained using Lymphoprep (07851, Stemcell) and then cultured with ImmunoCult-XF T Cell Expansion Medium (10981, Stemcell) and 10 ng/mL human recombinant IL2 (78036, Stemcell). ImmunoCult Human CD3/CD28 T Cell Activator (25 μL/mL; 10971, Stemcell) was added to the cell suspension and incubated for 3 days. Then the volume of the cell suspension was expanded every 2 days according to the manufacturer's instructions.

ELISpot assays were performed according to the manufacturer's instructions. In brief, we incubated 5 × 10⁴ cells with 10 μg/mL of the synthesized peptide in each well at 37°C in ELISpot plates (2110003, DAKEWE) for 20 hours. The plates were washed and incubated with antihuman IFNγ and then streptavidin–HRP. Finally, the streptavidin–HRP was stained using the AEC solution, and ELISpot plates were scanned and counted using an ImmunoSpot plate reader and associated software (Cellular Technologies, Ltd.).

The differentiation of monocytes into dendritic cells was conducted using the ImmunoCult Dendritic Cell Culture Kit (10985, Stemcell) according to the manufacturer's instructions. We first obtained PBMCs from a healthy donor and then incubated the cells with a serum-free medium at 37°C for 2 hours. Suspension cells were collected for T-cell expansion, and ImmunoCult dendritic cell differentiation medium with differentiation supplement was added to induce the differentiation of adherent monocytes. The medium was replaced every 2 days. Synthesized peptides were added at 50 μg/mL and allowed immature dendritic cells to uptake and process on day 5. Maturation supplement was added on day 6 to induce the maturation of dendritic cells. Then mature dendritic cells were obtained on day 8 and cocultured with expanded T cells for 24 hours.

For flow cytometry, antibodies used were anti-CD3 (563800, SK7, BD Biosciences, RRID: AB_2744384), anti-CD8 (557945, RPA-T8, BD Biosciences, RRID: AB_396953), anti-CD107a (328608, H4A3, BioLegend, RRID: AB_1186040), and anti–4-1BB (309804, 4B4-1, BioLegend, RRID: AB_314783). Data were collected using the flow cytometers BD LSRFortessa and Celesta (BD Biosciences), and analyzed by FlowJo (TreeStar).

We incubated 10⁷ cells with 10 μg/mL of peptide #19 or DMSO in 2 mL culture media at 37°C for 24 hours. 4-1BB⁺ cells were obtained using the CD137 MicroBead Kit (130-093-476, Miltenyi Biotec) according to the manufacturer's instructions.

Genome DNA from peripheral blood, tumor-adjacent tissues, and tumor tissues was collected for TCR-seq. We used the Multiplex PCR (MPCR) primers as described previously (84, 85), which includes 30 forward V primers and 13 reverse J primers, to amplify the rearranged CDR3 regions of TCRs and the sequencing was performed using a BGISEQ-500 platform (MGI Tech Co., Ltd.) following the manufacturer's protocol. Detailed data processing and TCR clonotype detection are described in Supplementary Materials and Methods.

A simple prognostic classification system was developed based on immune evasion ability for iCCA. Immune infiltration status was determined based on the unsupervised clustering of immune signatures by Danaher and colleagues (27). HLA LOH was called using LOHHLA (86) and CD8⁺ T-cell exhaustion scores were calculated with single-sample gene set enrichment analysis using the signature by Kang and colleagues (62), in which the top 10% of tumor samples were defined as extreme CD8 exhaustion.

Statistical analyses were conducted using SPSS 25.0 software (IBM). The Chi-square test, Fisher exact test, Mann–Whitney U test, Kruskal–Wallis test, and Student t test were used appropriately. Kaplan–Meier plots (log-rank test) and Cox proportional hazards regression analysis were used to describe survival curves and identify independent prognostic factors, respectively. A two-tailed P < 0.05 was considered significant. The results of group comparisons presented as box plots followed these rules: The maximum and minimum are indicated by the extreme points, the median is indicated by the thick horizontal line, and the first and third quartiles are indicated by box edges; whiskers indicate 1.5 times the interquartile range above the first quartile and below the third quartile. The schematic diagrams were drawn using BioRender.

Six public data sets were analyzed. We collected multiomics data of the CPTAC iCCA cohort (ref. 19; https://www.biosino.org/node/project/detail/OEP001105), somatic mutation data of Zou and colleagues (21) from GeneBank (SRP045202), and the data sets from TCGA. The somatic mutation data of Xiang and colleagues (14), Dong and colleagues (8), and Chen and colleagues (17) were also used to verify our results.

Sequence data have been deposited at the National Omics Data Encyclopedia, which is hosted by the Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences (accession no. OEP002560), and at CNSA (CNGB Sequence Archive; https://db.cngb.org/cnsa/) of the China National GeneBank DataBase (CNGBdb), with the accession code CNP0002045.

No disclosures were reported.

Y. Lin: Software, formal analysis, validation, methodology, writing–original draft, writing–review and editing. L. Peng: Software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. L. Dong: Software, formal analysis, validation, investigation, methodology, writing–original draft, writing–review and editing. D. Liu: Software, formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. J. Ma: Investigation, methodology, writing–review and editing. J. Lin: Validation, investigation, methodology. X. Chen: Software, formal analysis, methodology. P. Lin: Software, formal analysis, methodology. G. Song: Validation, investigation, writing–review and editing. M. Zhang: Validation, investigation, methodology, writing–review and editing. Y. Liu: Formal analysis, validation, investigation, writing–review and editing. J. Rao: Software, formal analysis, methodology. C. Wei: Software, formal analysis, writing–review and editing. Y. Lu: Software, formal analysis, writing–review and editing. S. Zhang: Software, funding acquisition, writing–review and editing. G. Ding: Investigation, writing–review and editing. Z. Peng: Software, formal analysis, writing–review and editing. H. Lu: Software, formal analysis. X. Wang: Investigation, writing–review and editing. J. Zhou: Investigation, writing–review and editing. J. Fan: Conceptualization, resources, supervision, funding acquisition, methodology, writing–original draft, writing–review and editing. K. Wu: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, writing–review and editing. Q. Gao: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

This study was supported by the National Natural Science Foundation of China (nos. 82130077, 81961128025, and 82002514) to Q. Gao and L. Dong; the Research Projects from the Science and Technology Commission of Shanghai Municipality (grants 21JC1410100, 21JC1401200, 20JC1418900, 19XD1420700, and 20YF1407400) to J. Fan, Q. Gao, S. Zhang, and L. Dong; Guangdong Provincial Key Laboratory of Human Disease Genomics (2020B1212070028) to K. Wu; Science, Technology and Innovation Commission of Shenzhen Municipality under grant JCYJ20170817145208120 to L. Peng; and the China National GeneBank. We thank Mr. Gaoxiang Ge (Chinese Academy of Sciences, Shanghai, China), Mr. Fuqiang Li, Ms. Qing Zhou, and Mr. Chen Ye (BGI-Shenzhen, Shenzhen, China) for their assistance and suggestions on this project.

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.