1. Introduction

Esophageal cancer (EC) is one of the most common malignancies in the digestive system, characterized by high incidence and poor prognosis. The malignant lesion is formed by dysplasia of esophageal epithelial or squamous epithelium, and its incidence and mortality rank the eighth and the sixth respectively among all malignant tumors.[1] The incidence of EC exhibits significant regional disparities, primarily concentrated in Central Asia, East Asia, South Asia, Southeast Asia, East Africa, and South Africa.[2,3] Early symptoms of EC are often inconspicuous, leading to late-stage diagnoses, making detection, treatment, and prevention challenging. Although some progress has been made in drug therapies, the overall 5-years survival rate for EC is <30%, and for advanced cases, it drops to a mere 5%.[4] Surgery remains the primary treatment option; however, many cases go undetected until the optimal surgical window has passed. Hence, finding early diagnostic indicators and risk assessment markers is crucial. Developing alternative therapies is also crucial for improving EC treatment.

Immunotherapy, which activates the patient’s own immune system to combat cancer cells,[5–7] has shown significant success in various cancer types.[8] Immune checkpoint inhibitors (ICIs), as representatives of immunotherapy, have been successfully used in the clinical treatment of esophageal squamous cell carcinoma and achieved satisfactory results.[9] Therefore, researching the mechanisms of immune cells in EC holds significant importance in the immunotherapy of EC.[10] Targeting patient-specific immune cell phenotypes can enhance the immune system’s ability to attack cancer cells, thereby improving treatment outcomes and reducing adverse events associated with immunotherapy. Understanding biological markers or genetic factors that predispose individuals to EC can aid in early diagnosis and risk assessment. These markers can also be utilized in the development of specific ICIs. Targeting the immune cell phenotype of patients enhances the ability of the immune system to attack cancer cells, thereby improving treatment outcomes.

Mendelian randomization (MR) employs single nucleotide polymorphisms (SNPs) as instrumental variables (IVs) to infer potential causal relationships between exposures and outcomes.[11,12] It adheres to the Mendelian gamete random allocation principle as well as the principle of free combination. This approach provides reasonable causal relationships. Therefore, this study employs a 2-sample MR approach, utilizing large-scale Genome-wide Association Studies (GWAS) databases to explore the causal relationship between immune cell phenotypes and EC risk, providing the foundation for potential immune therapy strategies and immune checkpoint inhibitor development. To our knowledge, there haven’t been previous studies using MR to assess the causal impact of immune cell traits on the risk of EC. This study offers a new direction for the treatment and disease management of EC.

2. Materials and methods 2.1. Study design

We carried out a 2-sample MR analysis to explore the causal relationship between 731 immunophenotype and EC. MR leverages genetic variants as proxies for risk factors. Therefore, to ensure the validity of causal inference, the IVs relies on 3 fundamental assumptions: Genetic variants have a direct connection with the exposure; Genetic variants are not associated with potential confounding factors between exposure and outcome; and Genetic variants do not influence the outcome through pathways unrelated to the exposure.[13] Furthermore, we conducted a reverse MR analysis to investigate the potential for reverse causality. All MR analyses were conducted using publicly available summary statistics, thus eliminating the need for additional ethical approval or informed consent.

2.2. Data source

The GWAS summary statistics for EC were derived from the GWAS Catalog, specifically identified by the accession number GCST90018621.[14] This study conducted a GWAS involving 160,589 participants of East Asian descent, primarily of Japanese origin, comprising 1388 cases and 159,201 controls. Furthermore, the GWAS summary statistics for a diverse range of immune traits were retrieved from the GWAS Catalog, with accession numbers spanning from GCST90001391 to GCST90002121.[15] The original GWAS for immune traits utilized data from 3757 European individuals, ensuring the absence of any overlap in the study cohorts. These datasets encompass an extensive collection of 731 immunophenotypes, spanning various aspects, including: Median fluorescence intensities (MFI) reflecting surface antigen levels (n = 389). Relative cell counts (n = 192). Absolute cell (AC) counts (n = 118). Morphological parameters (n = 32).

2.3. Instrumental variable selection

According to recent research, we performed association studies for genetic variations within each immune phenotype, with a significance level set at 1 × 10−5.[15,16] To mitigate bias stemming from linkage disequilibrium, we enforced the condition that the R2 value of SNPs closely associated with the exposure must be <0.001, within a 10,000 kb distance threshold.[17] We defined F > 10 as the absence of weak instrumental bias[18,19] in order to ensure a strong association of selected SNPs with EC risk, and avoid any potential weaknesses in the IVs. The formula for calculating F is: F=(N−K−1K)×(R21−R2), where N is the sample size in the exposure database, K is the sample size, R2 is the proportion of variance explained by SNPs in the exposure database.[20] The formula for calculating R2 is: R2=2×EAF×(1−EAF)×β2[(2×EAF×(1−EAF)×β2)+(2×EAF×(1−EAF)×K×SE×β2)], where EAF is the effect allele frequency, β is the allele effect value, K is the sample size, and SE is the Standard Error. The exposure and outcome datasets were merged, and incompatible alleles and palindromic SNPs were eliminated.[21] Phenotypes related to the remaining SNPs were identified through the Human Genotype-Phenotype Association Database (http://www.phenoscanner.medschl.cam.ac.uk/).[22] The remaining SNPs served as the ultimate IVs for exposure.[23]

2.4. Statistical analysis

The univariable MR analysis was performed using the R packages “TwoSampleMR.” Our analysis primarily utilized the inverse variance weighting (IVW) method.[24] The IVW method is an optimal estimation approach that assumes all genetic variations function as effective IVs, providing robust causal detection capabilities. However, the IVW method specifically requires that genetic variations solely influence the target outcome through exposure in the study. Despite our efforts to exclude known confounding SNPs, it is important to note that there may still exist numerous unknown confounding factors that could introduce genetic pleiotropy and introduce bias in effect estimates. Therefore, we employed 4 additional methods to assess the reliability and consistency of our results. These methods include MR-Egger regression,[25] the weighted median estimator,[26] mode-based simple estimation,[27] and mode-based weighted estimation. For each immune cell signature, we conducted MR analysis sequentially, and if the 5 different MR models produced consistent effect estimates, we considered the causal relationship between this immune phenotype and EC to be both stable and dependable.

The Cochran Q test and its corresponding P-values were employed to assess heterogeneity among the individual genetic variation estimates as well as among the selected IVs. If there was no significant heterogeneity (P > .05), we applied a fixed-effects model; otherwise, a random-effects model was employed.[20] We used the MR-Egger intercept to assess the effect of horizontal pleiotropy.[28] Sensitivity analysis was conducted using the “leave-one-out” method to assess the impact of individual SNPs on the causal relationship. The Mendelian randomization pleiotropy RESidual sum and outlier (MR-PRESSO) method were utilized as indicators to evaluate and correct the level of pleiotropy,[29] the MR-PRESSO was conducted using the R package “MRPRESSO” (Fig. 1).

Figure 1.:

The steps of Mendtlian randomization (MR) analysis. LD = linkage disequilibrium, IVW = inverse varianee weighting.

3. Results 3.1. Screening results of immunophenotypes associated with esophageal cancer risk

Using the IVW method in MR analysis, we conducted an analysis and selected 31 immunophenotypes associated with EC out of a total of 731 immune cell traits. Assessed 31 immunophenotypes and found weak IVs with F < 10, present in 3 immunophenotypes, whose ID is respectively ebi-a-GCST90001605 (rs76019287, F = 7.3461), ebi-a-GCST90001421(rs41494, F = 8.8348), and ebi-a-GCST90002071 (rs73018616, F = 4.2861) (Table S1, Supplemental Digital Content, https://links.lww.com/MD/M343). After removing weak IVs, performing MR analysis by the IVW method, we obtained 29 immunophenotypes associated with EC. To eliminate the influence of confounding factors on causality, we conducted a search on the Phenoscanner website for SNPs related to these 29 immunophenotypes and found no SNPs associated with confounding factors. Horizontal pleiotropy was assessed through MR-Egger regression, and the MR-Egger intercept for all 29 immunophenotypes had P-Values > .05, indicating no evidence of horizontal pleiotropy. The Cochran’s Q test assesses the heterogeneity of SNPs, and both the Q_P-values for IVW and MR-Egger are >.05, indicating the absence of heterogeneity. The results of the MR-PRESSO method all show P-values > .05, and no outliers were detected.

Next, we employed the MR-Egger regression method, WME, Mode-Based Simple Estimator, and Mode-Based Weighted Estimator to assess the reliability and consistency of the results. The analysis revealed inconsistencies for CD20− AC, CD25++ CD8br %CD8br, CD25hi CD45RA+ CD4 not Treg %T cell, CD66b on Gr MDSC, and DP (CD4+CD8+) %T cell compared to the results obtained using the IVW analysis. For CD20− AC, the IVW analysis showed P = .0377, OR = 1.0656, 95% CI = 1.004 to 1.1314, while the Simple Mode analysis showed P = .4670, OR = 0.8715, 95% CI = 0.6069 to 1.2514. For CD25++ CD8br %CD8br, the IVW analysis showed P = .0437, OR = 1.1867, 95% CI = 1.0049 to 1.4015, while in MR-Egger, P = .6422, OR = 0.8897, 95% CI = 0.5517 to 1.4350. For CD25hi CD45RA+ CD4 not Treg %T cell, the IVW analysis showed P = .0290, OR = 1.2049, 95% CI = 1.0193 to 1.4243, while in MR-Egger, P = .4720, OR = 0.8587, 95% CI = 0.5781 to 1.2754. CD66b on Gr MDSC had a P-value of .0154 in the IVW method, with OR = 0.9115 and 95% CI = 0.8456 to 0.9824. In the Simple Mode analysis, P = .1849, OR = 1.1528, 95% CI = 0.9447 to 1.4067. Finally, for DP (CD4+CD8+) %T cell, the IVW method resulted in P = .0167, OR = 0.6553, 95% CI = 0.4636 to 0.9265, while MR-Egger showed P = .8248, OR = 1.3178, 95% CI = 0.1337 to 12.9922 (Fig. 2, Table S2, Supplemental Digital Content, https://links.lww.com/MD/M344). Due to the inconsistent direction of the results of odds ratio (OR), which may be the result of potential outliers, the causal relationship between the above 5 immunophenotypes and EC needs to be verified by more studies in the future.

Figure 2.:

Forest plots showing the results of the 5 MR analysis methods showed an inconsistent causal relationship between 5 immune cell traits and EC. CI = confidence interval, EC = esophageal cancer, MR = Mendelian randomization.

3.2. Exploration of the causal effect of immunophenotypes on esophageal cancer

Using the IVW analysis, we found that 9 immunophenotypes types were positively associated with the risk of EC: CD20− %B cell, IgD− CD38dim %B cell, CD20− %lymphocyte, CD25 on IgD+ CD24+, CD25 on IgD− CD27−, CD27 on IgD+ CD24+ (B cell panel). CD3 on TD CD8br (Maturation stages of T cell panel). Mo MDSC AC (Myeloid cell panel). CD28+ CD45RA− CD8br AC (Treg panel). Specifically, The genetically predicted CD20− %B cell was found to have a positive association with EC, as indicated by the IVW method (OR = 1.1690, 95% CI = 1.0263–1.3315, P = .0187). Positive causal correlation between CD20% lymphocytes and EC. The OR of IVW was 1.0719 (95% CI = 1.0197–1.1268, P = .0064).

In the case of CD25 on IgD− CD27−, there was a positive association with EC (IVW OR = 1.1476, 95% CI = 1.0178–1.294, P = .0246). Similarly, CD25 on IgD+ CD24+ also is thought to increase the risk of EC (IVW OR = 1.1419, 95% CI = 1.0132–1.287, P = .0296). Moreover, CD27 on IgD+ CD24+ had a positive causality (IVW OR = 1.1006, 95% CI = 1.0034–1.2071, P = .0421). In addition, CD28+ CD45RA− CD8br AC showed a positive association with EC (IVW OR = 1.2406, 95% CI = 1.0555–1.4582, P = .0089). CD3 on TD CD8br was positively associated with EC, and the IVW analysis showed OR = 0.0471 (95% CI = 1.0016–1.2701, P = .0471). For IgD-CD38dim %B cells, a positive correlation with the EC. (IVW OR = 1.1914, 95% CI = 1.0043–1.4134, P = .0446). Finally, The analysis of Mo MDSC AC showed a IVW OR of 1.1001 (95% CI = 1.0022–1.2075, P = .0449), suggesting that Mo MDSC AC may contribute to the occurrence of EC. (refer to Fig. 3 and Table S3, Supplemental Digital Content, https://links.lww.com/MD/M345).

Figure 3.:

Forest plots showing positive causal relationships between 9 immune cell traits and EC. CI = confidence interval, EC = esophageal cancer.

All the results obtained from various methods were consistently in the same direction (refer to Fig. 3 and Table S3, Supplemental Digital Content, https://links.lww.com/MD/M345). Notably, there was no evidence of confounding heterogeneity, as confirmed by the Cochran’s Q test (P > .05, Table 1), and the leave-one-out test also supported this observation (Figure S1, Supplemental Digital Content, https://links.lww.com/MD/M348). Our findings further demonstrated that there was no significant evidence of horizontal pleiotropy impacting the causal relationship (MR-Egger intercept P > .05, Table 1). Additionally, no outliers were detected when using MR-PRESSO (P > .05, as shown in Table 1). The estimates of the causal effect were visually represented in scatter plots (Figure S2, Supplemental Digital Content, https://links.lww.com/MD/M349). Furthermore, funnel plots, which offer a visual assessment, showed no heterogeneity in the causal effect (Figure S3, Supplemental Digital Content, https://links.lww.com/MD/M350).

Table 1 - Sensitivity analysis of positive causal relationships between 9 immune cell traits and EC. Exposure Outcome Method Cochran’s Q Pleiotropy Q Q_df Q_pval Egger intercept SE P-value CD20- %B cell Esophageal cancer MR-Egger 12.84 11 0.30 0.025 0.034 .48 IVW 13.46 12 0.34 MR-PRESSO .41 CD20- % lymphocyte Esophageal cancer MR-Egger 7.26 11 0.78 0.025 0.018 .20 IVW 9.08 12 0.70 MR-PRESSO .57 CD25 on IgD- CD27- Esophageal cancer MR-Egger 14.86 10 0.14 -0.017 0.038 .66 IVW 15.17 11 0.18 MR-PRESSO .28 CD25 on IgD + CD24+ Esophageal cancer MR-Egger 8.07 11 0.71 -0.007 0.032 .83 IVW 8.12 12 0.78 MR-PRESSO .80 CD27 on IgD + CD24+ Esophageal cancer MR-Egger 10.14 12 0.60 -0.002 0.047 .97 IVW 10.14 13 0.68 MR-PRESSO .69 CD28 + CD45RA- CD8br AC Esophageal cancer MR-Egger 19.94 16 0.22 -0.020 0.032 .55 IVW 20.40 17 0.25 MR-PRESSO .27 CD3 on TD CD8br Esophageal cancer MR-Egger 10.72 10 0.38 -0.039 0.023 .12 IVW 13.82 11 0.24 MR-PRESSO .32 IgD- CD38dim %B cell Esophageal cancer MR-Egger 25.49 21 0.23 -0.009 0.022 .69 IVW 25.69 22 0.27 MR-PRESSO .26 Mo MDSC AC Esophageal cancer MR-Egger 4.32 12 0.98 -0.039 0.033 .26 IVW 5.72 13 0.96 MR-PRESSO .96

IVW = inverse variance weighted, MR-PRESSO = MR-pleiotropy residual sum and outlier.

A total of 15 immunophenotypes were identified as having a negative causal relationship with EC by the IVW method. This includes 4 were in the B cell panel, 1 in the cDC panel, 3 in the TBNK panel, 4 in the Maturation stages of T cell panel and 3 in the Treg panel. IgD+ CD38− %B cell was found to have a negative association with EC (IVW OR = 0.7807, 95% CI = 0.6161–0.9892, P = .0404). IgD− CD24− %lymphocyte also showed a negative causality through IVW (OR = 0.881, 95% CI = 0.7770–0.9990, P = .0482). CD19 on IgD− CD38dim exhibited a negative causality with EC, with an IVW (OR = 0.8181, 95% CI = 0.6725–0.9953, P = .0447). Similarly, CD20 on IgD+ CD24+ had a negative causality on EC (IVW OR = 0.8996, 95% CI = 0.8165–0.9911, P = .0323). CD62L- myeloid DC AC was negatively associated with EC (IVW OR = 0.8263, 95% CI = 0.6918–0.9869, P = .0352). CD4+ AC also demonstrated a negative association with EC (IVW OR = 0.8245, 95% CI = 0.6805–0.9989, P = .0487). Furthermore, Lymphocyte %leukocyte had a negative causality on EC (IVW OR = 0.9157, 95% CI = 0.8536–0.9823, P = .0139). CD3 on HLA-DR+ T cell was negatively causally linked to EC IVW OR = 0.9244 (95% CI = 0.8713–0.9807, P = .0092). Similarly, CD3 on CD45RA− CD4+ exhibited a negative causality with EC, having an IVW OR of 0.9012 (95% CI = 0.8202–0.9902, P = .0303). HVEM on naive CD4+ AC was negatively associated with EC (IVW OR = 0.8261, 95% CI = 0.7183–0.9502, P = .0075). And HVEM on CD45RA− CD4+ showed a negative causality with EC, having an IVW OR of 0.9244 (95% CI = 0.8618–0.9916, P = .0281). Moreover, CD4 on TD CD4+ also had a negative causality with EC, showing an IVW OR of 0.9111 (95% CI = 0.8393–0.9891, P = .0263). CD4 on CD4 Treg showed a negative causality on EC (IVW OR = 0.8634, 95% CI = 0.7802–0.9554, P = .0045). As well as the IVW method detected a negative correlation between CD4 on CD39+ resting Treg and EC (OR = 0.796, 95% CI = 0.6373–0.9942, P = .0443). Finally, for CD4 on activated & secreting Treg, there was a negative causality with EC (IVW OR = 0.9057, 95% CI = 0.8271–0.9917, P = .0324). (refer to Fig. 4 and Table S4, Supplemental Digital Content, https://links.lww.com/MD/M346).

Figure 4.:

Forest plots showing a negative causal relationship between 15 immune cell traits and EC. CI = confidence interval, EC = esophageal cancer.

In the same vein, all results obtained from various methods consistently pointed in the same direction (Fig. 4 and Supplementary Table 4, https://links.lww.com/MD/M346). There was no evidence of confounding heterogeneity (Cochran’s Q test P > .05, Table 2), and the leave-one-out test further supported this observation (refer to Figure S4, Supplemental Digital Content, https://links.lww.com/MD/M351). And the same, there was also no significant evidence of horizontal pleiotropy, as confirmed by the MR-Egger Intercept (P > .05, Table 2). Moreover, MR-PRESSO detected no outliers (P > .05, as shown in Table 2). A scatter plot depicting the causal effect see Figure S5, Supplemental Digital Content, https://links.lww.com/MD/M352. Furthermore, provide funnel plots that visually assess the heterogeneity of the causal effects refer to Figure S6, Supplemental Digital Content, https://links.lww.com/MD/M353.

Table 2 - Sensitivity analysis of negative causal relationship between 15 immune cell traits and EC. Exposure Outcome Method Cochran’s Q Pleiotropy Q Q_df Q_pval Egger intercept