Bidirectional Mendelian Randomization and Single‑Cell RNA Sequencing Reveal an NK Cell–Mediated Causal Link Between Endometriosis and Endometrial Cancer

Endometriosis is a benign estrogen-dependent inflammatory disorder that is characterized by the presence of ectopic endometrial implants. These implants most commonly occur within the pelvis but have also been identified in the upper abdomen, lungs, and central nervous system. 1 , 2 Endometriosis typically manifests as chronic pelvic pain (both cyclical and non-cyclical), dysmenorrhea, infertility, dyspareunia and dysuria. 3 The global Burden of Disease analyses indicate that although incidence and prevalence appear to have declined, overall burdens, such as mortality and disability-adjusted life years, have increased. Technological advances have reduced morbidity and disability among younger women; nevertheless, the mortality risk in older populations continues to rise. 4 The health-care costs, loss of work productivity, and reduced quality of life associated with endometriosis render it a gynecological public health priority. 5 , 6 An estrogen-driven endocrine milieu, chronic inflammation, immune dysregulation, and abnormalities in the epithelial–mesenchymal transition (EMT) and signaling pathways, such as Wnt and transforming growth factor (TGF-β), have been proposed to form a bidirectional network linking endometriosis and endometrial cancer, thereby influencing their pathogenesis and therapeutic responses. 7 Notably, immune imbalance, particularly perturbations in NK cell subsets, may constitute a mechanistic bridge between these two conditions. 8 , 9 As key effectors of innate immunity, NK cells in patients with endometriosis have been reported to exhibit increased cytotoxicity, and overexpression of inflammatory mediators such as SAMD9 and RGL2 has been associated with exacerbation of pain. 10 In contrast, in endometrial cancer, the dysregulated expression of inflammatory cytokines (eg, IL-1β and IL-6) and chemokines (eg, CXCL12 and CXCL27) appears to alter the tumor microenvironment in ways that may impair NK cell function and recruitment to tumor sites. 11 Collectively, alterations in NK cell subsets and functions constitute a critical immunological dimension linking research on endometriosis and endometrial cancer. 12 , 13 In recent years, single-cell RNA sequencing (scRNA-seq) has been widely employed to determine cellular lineage composition, cell states, and intercellular interaction landscapes within tissues at single-cell resolution, thereby enabling identification of disease-associated cell populations and alterations in transcriptional programs. 14–16 In parallel, genome-wide association studies (GWAS) have made significant progress at the population level in identifying susceptibility loci and regulatory variants associated with disease risk. 17 , 18 The integration of scRNA-seq with genomic approaches such as GWAS enables combining evidence across cell-type localization, molecular function, and genetic susceptibility. It may therefore help elucidate shared pathogenic mechanisms underlying endometriosis and endometrial cancer, and identify potential pathogenic pathways or therapeutic targets. In this study, scRNA-seq and MR approaches were employed to systematically investigate cellular and genetic commonalities between endometriosis and endometrial cancer (EC). By integrating single-cell sequencing data derived from patients, the proportions of each NK cell subset in both conditions were quantified, cell–cell communication analyses of NK cells were performed, and differential gene expression analysis of CD56dim NK cells was conducted. Pseudotemporal trajectory analysis was applied to delineate dynamic changes in key genes during cellular development. eQTL data were integrated with GWAS summary statistics for endometriosis and EC 19 and bidirectional MR analyses were performed to identify putative causal target genes affecting both diseases. Colocalization analysis was subsequently performed to confirm the shared localization of specific single-nucleotide polymorphisms (SNPs) with genes associated with endometriosis and EC, and candidate targets were validated in peripheral blood samples. Further analyses implicated BRI3 as a key gene in the pathogenesis of both the conditions.

Single-cell sequencing data were subjected to initial quality filtering to remove low-quality cells and genes ( Figure 1A ). Principal component analysis (PCA) was performed to extract major sources of variation and reduce data dimensionality. The Harmony algorithm was applied for batch-effect correction to minimize interbatch differences in the experimental data ( Figure 1B ). After batch correction, samples from different batches were better integrated into the data space and batch-associated variation was substantially reduced ( Figure 1C ). The extracted principal components accounted for a significant proportion of the observed variance ( Figure 1D and E ). Cells from different patient samples were clustered in the reduced UMAP space, and eight major lineages were identified based on the established marker gene sets: T cells, epithelial cells, NK cells, monocytes, macrophages, endothelial cells, smooth muscle cells, and B cells ( Figure 1F ). Figure 1 Preprocessing and initial quality control of single‑cell sequencing data. ( A ) Boxplots summarizing per‑cell quality‑control metrics: number of detected features (nFeature_RNA), total RNA counts (nCount_RNA), percentage of mitochondrial reads (percent.mt) and percentage of hemoglobin reads (percent.HB). ( B ) Principal component analysis (PCA) of the single‑cell RNA‑seq dataset. ( C ) Scatter plot of cells after batch‑effect correction using the Harmony algorithm. ( D ) Quantile–quantile (QQ) plot of p‑values computed for each principal component, used to assess deviation from the null distribution. ( E ) Elbow plot derived from PCA, used to determine the number of principal components to retain. ( F ) UMAP projection showing cell clusters annotated according to canonical marker genes. Image A shows violin plots for quality-control metrics: nFeatureRNA (0-10000), nCountRNA (0-150000) and percent.mt (0-75) across Control and Treatment groups, highlighting RNA and mitochondrial content differences. Image B is a PCA scatter plot with PC1 (0-100) and PC2 (-20 to 40), showing group separation. Image C displays a scatter plot using Harmony1 (0-80) and Harmony2 (-40 to 40), indicating batch-effect correction. Image D is a quantile plot with Original (0.000-0.100) and Theoretical Quantiles (0.0-0.3), showing deviation from the null distribution. Image E presents an elbow plot with Principal Component (0-30) and Standard deviation (2-10), suggesting retention of around 10 components. Image F includes four UMAP plots with umap1 (-10 to 10) and umap2 (-10 to 15), illustrating cell clusters like T and B cells, differentiated by color. The figure demonstrates preprocessing, PCA, batch correction and clustering of single-cell RNA data. Figure with violin, scatter, quantile, elbow plots and UMAP of single-cell RNA data. Preprocessing and initial quality control of single‑cell sequencing data. ( A ) Boxplots summarizing per‑cell quality‑control metrics: number of detected features (nFeature_RNA), total RNA counts (nCount_RNA), percentage of mitochondrial reads (percent.mt) and percentage of hemoglobin reads (percent.HB). ( B ) Principal component analysis (PCA) of the single‑cell RNA‑seq dataset. ( C ) Scatter plot of cells after batch‑effect correction using the Harmony algorithm. ( D ) Quantile–quantile (QQ) plot of p‑values computed for each principal component, used to assess deviation from the null distribution. ( E ) Elbow plot derived from PCA, used to determine the number of principal components to retain. ( F ) UMAP projection showing cell clusters annotated according to canonical marker genes. Following comparison of cell-cluster composition between endometriosis and EC patients and their respective normal controls, NK cells were found to represent the largest immune-cell fraction; therefore, NK cells were selected for focused analysis and were extracted for separate reprocessing. The same dimensionality reduction, clustering, and batch-correction procedures were applied to the NK cell subset, yielding good sample integration after correction ( Figure 2A and B ), and the PCA results were consistent with the global analysis ( Figure 2C and D ). Twelve NK cell subclusters were delineated ( Figure 2E ) and annotated according to canonical marker genes into CD56dim, CD56bright NK cells and NKT cells ( Figure 2G ). Comparative proportion analysis indicated that CD56dim NK cells were predominant in the endometrial compartments of both disease cohorts. NKT cell proportions decreased in both endometriosis and EC relative to controls. In contrast, the proportion of CD56bright NK cells diverged between conditions, increasing in endometriosis but decreasing in EC ( Figure 2F and H ). Figure 2 Preprocessing and initial handling of NK‑cell subset sequencing data. ( A ) PCA plot illustrating the distribution of samples by data source. ( B ) Scatter plot of NK‑cell data after batch‑effect correction using the Harmony algorithm. ( C ) QQ plot assessing the statistical significance of principal components. ( D ) Scree plot derived from principal component analysis. ( E ) UMAP projection displaying 12 identified NK‑cell subclusters annotated by canonical marker genes. ( F ) UMAP visualization showing the distributions of CD56bright NK cells, CD56dim NK cells and NKT cell subtypes within the clustering. ( G ) Expression profiles of marker genes across NK‑cell subpopulations. ( H ) Comparison of NK‑cell subpopulation proportions among endometriosis, endometrial cancer and healthy control groups. Panel A shows a PCA plot with PC 1 and PC 2 axes, illustrating sample distribution by data source. Panel B presents a scatter plot of NK-cell data post batch-effect correction using Harmony, with similar axes. Panel C displays a QQ plot assessing principal component significance, with expected and observed values. Panel D is a scree plot showing standard deviation across principal components, flattening after the initial few. Panel E includes UMAP projections for EC Control, EC Treatment, Endometriosis Control and Endometriosis Treatment, with UMAP 1 and UMAP 2 axes, highlighting 12 NK-cell subclusters. Panel F shows UMAP visualizations of CD56bright NK cells, CD56dim NK cells and NKT cells, indicating clustering differences by condition. Panel G is a dot plot with cell types on the x-axis and marker genes on the y-axis, showing expression levels and percentages. Panel H features a stacked bar chart comparing cell-type proportions across conditions, with CD56dim NK cells dominating. The panels collectively demonstrate dimensionality reduction, clustering and cell-type distribution differences among conditions, with CD56dim NK cells prevalent in disease cohorts and varying NKT cell proportions. NK-cell data: PCA, QQ, scree, UMAP, dot plot, stacked bar chart for clustering and cell-type proportions. Preprocessing and initial handling of NK‑cell subset sequencing data. ( A ) PCA plot illustrating the distribution of samples by data source. ( B ) Scatter plot of NK‑cell data after batch‑effect correction using the Harmony algorithm. ( C ) QQ plot assessing the statistical significance of principal components. ( D ) Scree plot derived from principal component analysis. ( E ) UMAP projection displaying 12 identified NK‑cell subclusters annotated by canonical marker genes. ( F ) UMAP visualization showing the distributions of CD56bright NK cells, CD56dim NK cells and NKT cell subtypes within the clustering. ( G ) Expression profiles of marker genes across NK‑cell subpopulations. ( H ) Comparison of NK‑cell subpopulation proportions among endometriosis, endometrial cancer and healthy control groups. To further investigate the determinants of the altered abundance and functional state of CD56dim NK cells, the intercellular communication of CD56dim NK cell subpopulations in endometriosis and endometrial cancer (EC) samples was evaluated using the CellChat package. Visualization of inferred interaction networks showed that in endometriosis, CD56dim NK cells exhibited pronounced interactions with monocytes and NKT cells, whereas in EC, these interactions were more extensive and tightly connected. Compared with endometriosis, associations between CD56dim NK cells and epithelial and endothelial cells were stronger in EC, suggesting enhanced cell–cell communication by CD56dim NK cells within the tumor microenvironment ( Figure 3A and B ). Ligand–receptor pathway analysis revealed distinct signaling enrichment between the conditions. In endometriosis, CD56dim NK cells were predominantly associated with the COL1A1–CD44 and COL1A2–CD44 signaling axes, and engaged monocytes, macrophages, NKT cells, and T cells through these pathways ( Figure 3C ). In contrast, in endometrial cancer, CD56dim NK cells were mainly enriched for MIF–CD74–CXCR4 and MIF–CD74–CD44 signaling, interacting primarily with monocytes and B cells while forming broader but comparatively weaker contacts with other immune populations via COL6A2–CD44 ( Figure 3D ). Figure 3 CD56dim NK‑cell communication analysis. ( A ) Cell–cell interaction network of CD56dim NK cells in lesional endometrial samples from endometriosis, visualized using Cell Chat. ( B ) Cell–cell interaction network of CD56dim NK cells in endometrial cancer lesion samples, visualized using Cell Chat. ( C ) Dot plot showing enrichment of ligand–receptor signaling pathways mediating interactions between CD56dim NK cells and other cell types in endometriosis single‑cell samples. ( D ) Dot plot of ligand–receptor pathway enrichment for CD56dim NK‑cell interactions in endometrial cancer samples. The image A showing a cell–cell interaction network diagram for CD56dim NK cells. No axes or units shown. Nodes are labeled Endothelial cells, Epithelial cells, Macrophages, Monocytes, NKT cells, Smooth muscle cells, T cells, B cells, CD56bright NK cells and CD56dim NK cells. Links originate from CD56dim NK cells to other nodes and each link has a numeric label: to CD56bright NK cells 43; to Endothelial cells 8; to Epithelial cells 7; to Macrophages 7; to Monocytes 14; to NKT cells 19; to Smooth muscle cells 10; to T cells 10; to B cells 7; and a self-loop at CD56dim NK cells labeled 6. The largest label is 43 to CD56bright NK cells. The image B showing a cell–cell interaction network diagram for CD56dim NK cells. No axes or units shown. Nodes are the same labeled cell types. Link labels from CD56dim NK cells are: to Epithelial cells 50; to Endothelial cells 43; to Macrophages 10; to Monocytes 41; to NKT cells 35; to Smooth muscle cells 28; to T cells 18; to B cells 14; to CD56bright NK cells 23; and a self-loop labeled 21. The largest label is 50 to Epithelial cells. The image C showing a dot plot of ligand–receptor pathway enrichment for CD56dim NK-cell interactions. X-axis label: Ligand minus Receptor; no units shown. Y-axis label: Cell minus Cell Communication; no units shown. A legend labeled Communication Probability ranges from 0.4, 0.5, 0.6, to 0.7. Dot color follows this probability scale; dot size varies and represents a second magnitude encoding, with larger dots indicating larger plotted magnitude. X-axis categories visible: APP minus CD74, COL1A1 minus CD44, COL1A2 minus CD44, COL6A2 minus CD44, MIF minus CD74 minus CD44 and MIF minus CD74 minus CXCR4. Multiple large, high-probability dots appear under COL1A1 minus CD44 and COL1A2 minus CD44 across several Cell minus Cell Communication rows; a single prominent dot appears under MIF minus CD74 minus CXCR4. The image D showing a dot plot of ligand–receptor pathway enrichment for CD56dim NK-cell interactions. X-axis label: Ligand minus Receptor; no units shown. Y-axis label: Cell minus Cell Communication; no units shown. The Communication Probability legend again ranges from 0.4 to 0.7 and dot size varies as a second magnitude encoding. X-axis categories visible include COL6A2 minus CD44, COL1A1 minus CD44, COL1A2 minus CD44, MIF minus CD74 minus CD44 and MIF minus CD74 minus CXCR4. Many small, low-probability dots align under COL6A2 minus CD44 across numerous Cell minus Cell Communication rows. Two large, high-probability dots appear under MIF minus CD74 minus CXCR4 and MIF minus CD74 minus CD44. Four plots of CD56dim NK-cell interaction networks and ligand–receptor enrichment dot plots. CD56dim NK‑cell communication analysis. ( A ) Cell–cell interaction network of CD56dim NK cells in lesional endometrial samples from endometriosis, visualized using Cell Chat. ( B ) Cell–cell interaction network of CD56dim NK cells in endometrial cancer lesion samples, visualized using Cell Chat. ( C ) Dot plot showing enrichment of ligand–receptor signaling pathways mediating interactions between CD56dim NK cells and other cell types in endometriosis single‑cell samples. ( D ) Dot plot of ligand–receptor pathway enrichment for CD56dim NK‑cell interactions in endometrial cancer samples. Eighteen key marker genes were identified by differential expression analysis of the patient-derived CD56dim NK cell subpopulations. These genes were treated as exposures and their corresponding eQTLs were employed as instrumental variables to model gene expression as quantitative traits. EC GWAS summary statistics from the ebi-a-GCST90018838 dataset were used as the outcome data. Five MR methods were applied: inverse-variance weighted (IVW), MR–Egger regression, simple mode, weighted median, and weighted mode. Furthermore, Cochrane’s Q test was performed to assess horizontal heterogeneity, while the MR-Egger intercept and MR-PRESSO global test were utilized to evaluate potential directional pleiotropy ( Supplementary Table S1 ). Scatter plots were generated to display the effect sizes and significance levels of candidate marker genes ( Figure 4A ). Figure 4 Mendelian randomization and colocalization analyses. ( A ) Volcano plot showing the Mendelian randomization (MR) results. ( B ) Forest plot summarizing MR estimates for all genes with IVW P < 0.05. (P-values less than 0.05 in IVW and WM are shown in bold). ( C ) Results of reverse Mendelian randomization analyses. ( D ) Forest plot showing SNP‑level effect sizes on BRI3 expression. ( E ) Funnel plot of SNPs associated with BRI3. ( F ) Scatter plot displaying the linear regression trend between genetically predicted BRI3 expression and endometrial cancer risk. ( G ) Colocalization analysis between BRI3 eQTLs and endometrial cancer GWAS signals. The infographic consists of multiple panels related to Mendelian randomization and colocalization analyses. A) Volcano plot showing Mendelian randomization results with log2(OR) and -log10(P-value) axes, highlighting genes like CAPG and BRI3. B) Forest plot summarizing MR estimates for genes with IVW P < 0.05, listing exposures, outcomes, SNP counts, methods, P-values and ORs with confidence intervals. C) Results of reverse Mendelian randomization analyses for endometrial cancer with exposures like BRI3, CAPG and ITGAX. D) Forest plot showing SNP-level effect sizes on BRI3 expression with MR-Egger and IVW methods. E) Funnel plot of SNPs associated with BRI3, displaying 1/SE against IVW. F) Scatter plot showing linear regression trend between SNP effect on BRI3 and SNP effect on endometrial cancer risk. G) Colocalization analysis between BRI3 eQTLs and endometrial cancer GWAS signals, with GWAS summary statistics and eQTL -log10(P-values) displayed. Infographic showing Mendelian randomization and colocalization analyses related to endometrial cancer risk. Mendelian randomization and colocalization analyses. ( A ) Volcano plot showing the Mendelian randomization (MR) results. ( B ) Forest plot summarizing MR estimates for all genes with IVW P < 0.05. (P-values less than 0.05 in IVW and WM are shown in bold). ( C ) Results of reverse Mendelian randomization analyses. ( D ) Forest plot showing SNP‑level effect sizes on BRI3 expression. ( E ) Funnel plot of SNPs associated with BRI3. ( F ) Scatter plot displaying the linear regression trend between genetically predicted BRI3 expression and endometrial cancer risk. ( G ) Colocalization analysis between BRI3 eQTLs and endometrial cancer GWAS signals. BRI3, CAPG, and ITGAX exhibited pronounced effects and statistical significance and were therefore selected for further investigation. The MR–Egger regression and Cochran’s Q tests did not detect evidence of horizontal pleiotropy or heterogeneity for BRI3, CAPG, or ITGAX (all P-values > 0.05). The results from the MR analyses are summarized in a forest plot to provide an overall visualization of the inferred causal relationships ( Figure 4B ). IVW estimates indicated that a higher genetically predicted expression of BRI3 was associated with an increased EC risk (BRI3: OR 1.1601, 95% CI 1.0564–1.2738; P = 0.002). In contrast, higher genetically predicted expression of CAPG and ITGAX was associated with reduced EC risk (CAPG: OR 0.8415, 95% CI 0.7525–0.9411; P = 0.002; ITGAX: OR 0.9308, 95% CI 0.8857–0.9783; P = 0.005). Reverse (bidirectional) MR analyses identified evidence of a reverse causal relationship only for ITGAX (IVW, P < 0.05), consistent with a potential disease-suppressive role. No reverse causality was observed for BRI3 or CAPG ( Figure 4C ). Given the magnitude of the estimated effect, BRI3 was prioritized for follow-up analyses. The MR results were consistent with the role of increased BRI3 expression in promoting EC risk ( Figure 4F ). Inspection of the single-SNP effect estimates ( Figure 4D and E ) showed that most SNP effects lay to the right of the zero axis, supporting the effect consistency and lending further support to the pathogenic contribution of BRI3. Funnel-plot inspection revealed a symmetric distribution of instrumental SNP effects without notable asymmetry, suggesting an absence of substantial heterogeneity ( Figure 4D and E ). Colocalization analysis identified rs4506131 as a significant eQTL for BRI3, which also showed a significant association in EC GWAS. Therefore, the spatial overlap of the eQTL and GWAS signals suggests that genetically determined variations in BRI3 expression are likely related to genetic susceptibility to EC ( Figure 4G ). Heatmap visualization showed that BRI3 expression was markedly increased in EC single-cell datasets from diseased samples, whereas CAPG and ITGAX exhibited the opposite pattern, with higher expression in control samples and lower expression in diseased samples ( Figure 5A ). Boxplot comparisons across cell subpopulations indicated that BRI3 expression in NK cells was generally elevated in the disease groups, CAPG and ITGAX expression levels in NK cells were higher in the control group ( Figure 5B–D, H ). Cell-type–specific expression mapping revealed that BRI3 expression was concentrated in CD56dim and CD56bright NK cells, whereas CAPG and ITGAX were predominantly expressed in NKT cells, indicating distinct cellular expression patterns of these candidate genes ( Figure 5E–G ). Interactions between BRI3-positive (BRI3+) and BRI3-negative (BRI3−) CD56dim NK cells were analyzed and visualized. Both BRI3+ and BRI3− CD56dim NK cells engaged with T cells, NKT cells, and monocytes. However, the interaction strength of BRI3− CD56dim NK cells was comparatively weaker. Notably, BRI3+ CD56dim NK cells also exhibited measurable interactions with endothelial cells ( Figure 6A and B ). Ligand–receptor and pathway-level analyses indicated that BRI3− CD56dim NK cells principally communicated with T cells, macrophages, and smooth muscle cells via the COL6A2–CD44 axis, and with B cells via the MIF–CD74–CXCR4 pathway ( Figure 6C ). Figure 5 Gene expression of BRI3, CAPG and ITGAX. (The asterisks represent the level of significance, where **p < 0.01, ***p < 0.001). ( A ) Heatmap displaying expression levels of BRI3, CAPG and ITGAX in healthy controls (blue) and patients with endometrial cancer (EC) (red). ( B ) Boxplot comparing BRI3 expression between healthy controls and the EC group. ( C ) Boxplot comparing CAPG expression between healthy controls and the EC group. ( D ) Boxplot comparing ITGAX expression between healthy controls and the EC group. ( E ) UMAP projection illustrating the distribution of BRI3 expression across NK‑cell subclusters. ( F ) UMAP projection illustrating the distribution of CAPG expression across NK‑cell subclusters. ( G ) UMAP projection illustrating the distribution of ITGAX expression across NK‑cell subclusters. ( H ) Dot plot summarizing mean expression and percentage of expressing cells for BRI3, CAPG and ITGAX across annotated cell types. The image A shows a heatmap of gene expression levels for BRI3, CAPG and ITGAX in control and treatment groups. BRI3 expression is higher in the treatment group, while CAPG and ITGAX are higher in the control group. The image B shows a boxplot comparing BRI3 expression across cell types, indicating higher levels in the treatment group. The image C shows a boxplot for CAPG expression, with higher levels in the control group. The image D shows a boxplot for ITGAX expression, also higher in the control group. The image E shows a UMAP projection of BRI3 expression across NK cell subclusters, highlighting CD56dim and CD56bright NK cells. The image F shows a UMAP projection for CAPG expression, focusing on NKT cells. The image G shows a UMAP projection for ITGAX expression, also focusing on NKT cells. The image H shows a dot plot summarizing mean expression and percentage of expressing cells for BRI3, CAPG and ITGAX across cell types in control and treatment groups. Analyze BRI3, CAPG, ITGAX in control/treatment via heatmap, boxplots, UMAP, dot plot. Figure 6 Potential physiological functions of BRI3 and other key genes. ( A ) Interaction network of BRI3‑positive (BRI3+) CD56dim NK cells with other cell types. ( B ) Interaction network of BRI3‑negative (BRI3−) CD56dim NK cells with other cell types. ( C ) Dot plot visualizing ligand–receptor pathways mediating interactions between BRI3+ and BRI3− CD56dim NK cells and other cell types. ( D ) Metabolic pathway enrichment analysis comparing BRI3+ versus BRI3− CD56dim NK cells. ( E ) Histogram showing the distribution of binarized gene expression values used for pseudotime analysis. ( F ) Expression trajectories of key genes along the pseudotime axis. ( G ) Scatter plot illustrating the relationship between overall gene expression levels and pseudotime. Image A shows BRI3-positive CD56dim NK cells interacting strongly with Endothelial cells (33) and Macrophages (24). Image B depicts BRI3-negative CD56dim NK cells with weaker interactions, strongest with Endothelial cells (19) and Macrophages (13). Image C presents a dot plot of communication probability across cell-cell and ligand-receptor pairs, ranging from 0.00 to 0.07. Image D illustrates pathways like Kaposi sarcoma-associated herpesvirus infection, with gene count sizes and significance from 0.0 to 0.8. Image E features a gene expression histogram, peaking near 0.0 with a threshold at 0.2. Image F is a scatter plot of gene labels over a pseudo-timeline, showing distribution across fit quality. Image G includes a scatter plot with marginal histograms, indicating a positive association between time and gene, with statistical details like r-student, p-value and confidence intervals. Seven plots: BRI3+/- CD56dim NK cell networks, ligand-receptor, pathways, gene expression, pseudotime, scatter. Gene expression of BRI3, CAPG and ITGAX. (The asterisks represent the level of significance, where **p < 0.01, ***p < 0.001). ( A ) Heatmap displaying expression levels of BRI3, CAPG and ITGAX in healthy controls (blue) and patients with endometrial cancer (EC) (red). ( B ) Boxplot comparing BRI3 expression between healthy controls and the EC group. ( C ) Boxplot comparing CAPG expression between healthy controls and the EC group. ( D ) Boxplot comparing ITGAX expression between healthy controls and the EC group. ( E ) UMAP projection illustrating the distribution of BRI3 expression across NK‑cell subclusters. ( F ) UMAP projection illustrating the distribution of CAPG expression across NK‑cell subclusters. ( G ) UMAP projection illustrating the distribution of ITGAX expression across NK‑cell subclusters. ( H ) Dot plot summarizing mean expression and percentage of expressing cells for BRI3, CAPG and ITGAX across annotated cell types. Potential physiological functions of BRI3 and other key genes. ( A ) Interaction network of BRI3‑positive (BRI3+) CD56dim NK cells with other cell types. ( B ) Interaction network of BRI3‑negative (BRI3−) CD56dim NK cells with other cell types. ( C ) Dot plot visualizing ligand–receptor pathways mediating interactions between BRI3+ and BRI3− CD56dim NK cells and other cell types. ( D ) Metabolic pathway enrichment analysis comparing BRI3+ versus BRI3− CD56dim NK cells. ( E ) Histogram showing the distribution of binarized gene expression values used for pseudotime analysis. ( F ) Expression trajectories of key genes along the pseudotime axis. ( G ) Scatter plot illustrating the relationship between overall gene expression levels and pseudotime. Comparative enrichment analysis of BRI3+ versus BRI3− CD56dim NK cells revealed distinct pathway-activity profiles. BRI3+ CD56dim NK cells were predominantly enriched in pathways related to cellular senescence and proteoglycans in cancer, and showed enrichment in MAPK and PI3K–Akt signaling. In contrast, BRI3− NK cells exhibited fewer and weaker enrichments, primarily in pathways, such as human cytomegalovirus infection and prion disease ( Figure 6D ). For dynamic analysis, gene expression data were binarized and pseudotime trajectories were inferred to investigate temporal expression patterns ( Figure 6E ). Pseudotemporal profiling showed that most immune-related genes (eg, GZMB, CCL5, and CD247) were highly expressed at early pseudotime stages, whereas genes implicated in ferroptosis and signal regulation (eg, PRDX1, YWHAE, and SLC7A11) were markedly upregulated in the intermediate to late stages. The candidate gene BRI3 reached a fitted peak at mid-pseudotime, suggesting its potential regulatory role in CD56dim NK cell activation or functional reprogramming ( Figure 6F ). Scatter plot analysis confirmed a significant positive correlation between gene expression levels and pseudotime progression (P < 0.01) ( Figure 6G ).

Data for Mendelian randomization and colocalization analyses were obtained from the IEU OpenGWAS project ( https://gwas.mrcieu.ac.uk/ ). These comprised eQTL data corresponding to 18 genes that were significantly differentially expressed between CD56dim NK cells and other endometrial cell subpopulations, together with EC GWAS summary statistics from the ebi-a-GCST90018838 dataset. For colocalization analyses, eQTL data for BRI3, CAPG, and ITGAX and GWAS summary statistics for LA were retrieved from the IEU OpenGWAS project. scRNA-seq and microarray datasets were sourced from the Gene Expression Omnibus (GEO) database. Single-cell RNA sequencing datasets were retrieved from the GEO database. Specifically, two independent cohorts were included: GSE213216 , related to endometriosis, and GSE173682 , related to endometrial cancer. To ensure consistency and statistical rigor, three diseased samples and one healthy control sample were curated from each dataset for downstream analysis. Comparative transcriptomic analyses were subsequently performed between the disease groups and their respective healthy controls to identify differentially expressed features. Endometrial cell expression profiles were derived from scRNA-seq data. Dimensionality reduction and downstream processing were performed in R, using the Seurat package. Initial quality-control filtering retained cells with between 200 and 4000 detected genes and a mitochondrial gene proportion below 10%; cells and genes failing these criteria were removed as low quality. The filtered high-dimensional dataset was normalized using standard procedures implemented in Seurat. Principal component analysis (PCA) was applied to extract major sources of variation and convert the data into a lower-dimensional representation. The results were projected onto a two-dimensional space for visualization. Batch effects were corrected using the limma and sva R packages. Based on the elbow plot, the first 30 principal components were selected at the inflection point for the subsequent analysis. Graph-based clustering was performed using the Louvain algorithm and the cluster structure was visualized in two dimensions using the UMAP algorithm. Among all endometrial cell populations, the most abundant immune cell population, identified as NK cells by canonical marker expression, was extracted for focused analysis. The dimensionality reduction and clustering procedures confirmed the use of 30 principal components for this subset. Differential gene expression analysis among NK cell subtypes was conducted using the Seurat package, and differentially expressed genes were identified based on statistical significance (P-value) and log2 fold change. Each cell cluster was annotated based on differential expression results and previously reported marker genes. UMAP visualization was used to display the cluster organization and corroborate the biological relevance of the annotations. The intercellular communication between CD56dim NK cells was analyzed using the CellChat R package. Ligand–receptor interactions were inferred from expression data to identify signaling trends and key pathways involved in signal transmission between CD56dim NK cells and other cell types. Putative ligand–receptor pairs were predicted, and their interaction strength and directionality between different cell populations were quantified. The results were visualized using pseudotemporal trajectory plots and communication network maps to reveal the complex signaling relationships among cell types. Differential expression analysis of CD56 dim NK cell subpopulations was performed using the Seurat package. Genes with an absolute log2 fold change (|log2FC|) > 2 and an adjusted p-value (p_adj) < 0.05 were considered significantly differentially expressed. A total of 18 genes meeting these criteria were defined as key marker genes and selected as candidate exposures. These genes were treated as exposures, and eQTLs were used as instrumental variables to model gene expression as a quantitative trait and examine the relationship between genetic variation and gene expression. Summary statistics from the ebi-a-GCST90018838 dataset were used as outcome data for endometrial cancer (EC). Mendelian randomization (MR) analyses were conducted using the TwoSampleMR R package, applying inverse-variance weighted (IVW), MR-Egger regression, and weighted median methods. Heterogeneity and horizontal pleiotropy were assessed to evaluate the robustness of the findings. Comparative analyses between patients with EC and healthy controls were conducted to assess the association between CD56^dim NK cell–specific marker genes and EC. Following the initial MR result visualization, BRI3 was selected from the 18 candidate genes for further focused MR analyses based on its statistical significance and effect size. Forest plots, scatter plots, and funnel plots were generated to visualize and evaluate the results. Reverse MR analyses were also performed to assess the possibility of reverse causation and test the consistency of the causal estimates. For genetic-association validation and colocalization analysis, the BIMR and coloc R packages were employed to identify significant loci and test whether eQTL signals and GWAS signals shared a common causal variant. Colocalization analyses were specifically applied to the BRI3 eQTLs and EC GWAS data. Regional association plots were used to visualize the colocalization results, revealing the relative association strength of the loci under different conditions and the probability of shared genetic causal variants. This study exclusively utilized publicly available datasets from the Gene Expression Omnibus and IEU OpenGWAS. Although no new human subjects were directly involved, the use of these data was conducted under the oversight of the Institutional Review Board of Guangdong Medical University (Approval No. PJB2025-362-01). All data were handled in accordance with institutional and national ethical guidelines.

By integrating bidirectional Mendelian Randomization and single-cell RNA sequencing data, this study systematically investigated the potential causal relationships and shared pathogenic mechanisms linking endometriosis and endometrial cancer. Focusing on the CD56dim NK cell subset, we identified BRI3, CAPG, and ITGAX as key marker genes, with genetically predicted BRI3 expression positively associated with EC risk, while CAPG and ITGAX showed inverse associations. Functional analyses suggest that BRI3 may regulate immune responses and mediate the functional reprogramming of CD56dim NK cells, contributing to immune alterations that facilitate disease progression. These findings provide a mechanistic basis at the cellular and genetic levels for the link between endometriosis and EC and offer a theoretical foundation for the development of targeted diagnostic and therapeutic strategies.

Endometriosis is a common and often recalcitrant gynecological disorder that imposes chronic, intractable pain and markedly reduces quality of life, adversely affecting daily functioning and mental health, and, in some cases, fertility. EC, a malignant neoplasm arising from the endometrium, may progress rapidly if not detected early, increasing therapeutic complexity and mortality risk and substantially impairing patient well being. Therefore, elucidating the biological mechanisms underlying endometriosis and EC is of clinical importance. Although the pathogenic mechanisms linking endometriosis and EC remain incompletely defined, both conditions share notable similarities, including strong endometrial involvement, influence of endocrine factors, and obesity. Previous studies have indicated that immune microenvironmental dysregulation is present in the endometrial tissues of both endometriosis and EC, 20 , 21 and that such alterations contribute to disease initiation and progression. Clinical data further suggest an elevated risk of EC development in patients with endometriosis. 22 , 23 Genomic investigations have identified cancer-associated alterations in some endometriotic lesions, supporting the hypothesis that endometriosis may be a precursor for endometrioid and clear-cell ovarian carcinomas. 24–26 The common etiological contributors to endometriosis and EC include estrogenic stimulation and chronic inflammation. 27 Dorien et al reported the aberrant expression of NK‑cell receptors and altered cytokine production in the pelvic environment of patients with endometriosis, reinforcing the role of NK‑cell dysfunction in disease pathogenesis. 28 Using flow cytometry, Eidukaite et al demonstrated reduced levels of cytotoxic CD56+ NK cells in peritoneal fluid from endometriosis samples, increased Fas (CD95) expression in early stage (I/II) endometriosis, and elevated early activation marker CD69, indicating substantial shifts in NK‑cell subset composition relative to normal endometrium; upregulation of pro-apoptotic genes likely contributes to these changes. 29 NK cells also play a critical role in the development of EC. A case report described the successful application of cord‑blood–derived NK (CB‑NK) cell immunotherapy in an EC patient, which enhanced NK‑cell activity and suppressed tumor growth and metastasis. 30 Degos et al found that intratumoral CD103+ NK cells in EC expressed higher levels of co-inhibitory molecules, resulting in reduced effector function and impaired recruitment to tumor sites. 11 In advanced EC, myeloid-derived suppressor cell (MDSC) infiltration hinders CD8+ T cell accumulation, and MDSC‑associated overexpression of arginase and inducible nitric oxide synthase (iNOS) reduces NK‑cell antibody-dependent cellular cytotoxicity (ADCC). 31 Tumor-derived chemokine alterations in the tumor microenvironment (TME) further promote tumor growth and dissemination; for example, increased utilization of CXCL12 by CXCR4‑expressing EC cells may deplete CXCL12 available for recruitment of antitumor cytotoxic NK cells. 11 Consequently, although NK cells are present within the EC TME, their antitumor responses are less effective than those of NK cells in the adjacent healthy tissues. Compared with peritumoral healthy tissue, EC lesions contain fewer NK cells with lower cytotoxic potential, and a subset of tumor-associated NK cells exhibits a phenotype resembling decidual NK cells (expressing CD56, CD9, CD49a, and CXCR3), thereby adopting immunosuppressive properties that favor tumor growth rather than cell death. 32 Approximately 40% of NK cells within the EC TME express CD103, a marker of tissue-resident NK cells with attenuated cytotoxicity. Diminished chemotactic signals in the TME reduce NK‑cell recruitment and further promote an immunosuppressive milieu. 11 Accumulating evidence indicates that alterations in NK‑cell number, phenotype, and function constitute important immunological links between endometriosis and EC. Elucidation of these changes may improve our understanding of the overlapping pathobiology of the two conditions and have potential implications for the diagnosis, risk stratification, and development of targeted therapies. We investigated the potential interactions between endometriosis and EC with the aim of developing novel theoretical frameworks and therapeutic strategies for both diseases. Using single-cell transcriptomic data, we focused on NK cells, particularly CD56dim NK cells and identified 19 endometriosis-associated marker genes. MR analyses were performed using eQTLs as instrumental variables to assess the relationships between candidate genes and the risk of endometrial cancer. The MR results indicated a positive association between genetically predicted BRI3 expression and EC risk, whereas the genetically predicted expression of CAPG and ITGAX was inversely associated with EC risk. BRI3 encodes BRI3 protein and is located on chromosome 7p14.3. Initially identified in the brain tissue, BRI3 is expressed in multiple organs, including the brain, liver, lungs, kidneys, and heart. 33 BRI3 has been primarily studied in neurobiology and implicated in the neurodegenerative features of Alzheimer’s. 34 Recently, several groups have begun to explore BRI3 in the field of oncology. In hepatocellular carcinoma models, the activation of Wnt/β-catenin signaling has been associated with the upregulation of BRI3, suggesting a functional role in liver cancer. 35 Jun Chen et al combined bioinformatic analyses with experimental validation to show that BRI3 regulates lipid autophagy in glioblastoma, thereby maintaining lipid homeostasis and metabolic adaptability and promoting tumor cell survival under metabolic stress; this work supports BRI3 as a potential metabolic therapeutic target. 36 Collectively, these findings indicate that BRI3 may serve as a biomarker for certain malignancies and could inform prognosis and therapeutic responses, although research on BRI3 in cancer remains at an early stage. In our in situ single-cell analyses of endometrial tissue from endometriosis and EC samples, batch-corrected data revealed concordant clustering trends across immune cell subpopulations, with NK cells comprising the dominant component. Prior studies have reported the enrichment of CD56dim NK cells in the peripheral blood 37 and the predominance of CD56bright NK cells within the endometrium. 38 Paradoxically, our data showed a clear predominance of CD56dim NK cells within the endometrial NK‑cell compartment in both endometriosis and EC. Guided by this observation, we characterized the roles of NK cells under these two conditions. In endometriotic lesions, CD56dim NK cells primarily communicate with other immune cells via extracellular matrix–receptor axes (notably, COL1A1/COL1A2–CD44, supplemented by MIF–CD74–CXCR4 and MIF–CD74–CD44), forming weak‑to‑moderate networks centered on matrix remodeling and inflammatory regulation with monocytes, macrophages, NKT cells, and T cells. In EC, CD56dim NK‑cell communication was markedly enhanced: interactions with monocytes/NKT cells were stronger, and significant crosstalk with epithelial and endothelial cells was observed, implicating CD56dim NK cells in angiogenesis, epithelial–immune interplay, and immune-escape reprogramming within the tumor microenvironment. When EC single-cell data were stratified by BRI3 expression, BRI3+ CD56dim NK cells exhibited stronger interactions with T, NKT, and monocyte populations, and enhanced communication with endothelial cells. BRI3+ CD56dim NK cells were enriched in pathways related to cellular senescence, proteoglycans in cancer, MAPK signaling, and PI3K–Akt signaling, consistent with an increased capacity for signal integration and tissue remodeling in the tumor milieu. Pseudotime analysis indicated an early phase dominated by cytotoxicity and recruitment programs that progressed toward mid‑to‑late phases characterized by stress responses and metabolic regulation; BRI3 expression peaked at intermediate pseudotime, suggesting that BRI3 may function as a regulatory node mediating the transition of CD56dim NK cells from inflammatory responses toward tumor‑adaptive functional reprogramming. Mendelian randomization leverages the random allocation of genotypes to infer the causal effects of exposures on outcomes, and has the advantage of mitigating confounding and reverse causation because genetic variants are fixed at conception. 39 Limitations of MR include generally small effect sizes of genetic variants on most risk factors, which can reduce statistical power and increase the risk of false-negative findings; accordingly, MR results should be interpreted alongside complementary experimental and epidemiological evidence to increase robustness. 40 In this study, colocalization analysis combining BRI3 eQTLs and EC GWAS summary statistics revealed that the single-nucleotide polymorphism rs4506131 shows strong signals in both BRI3 eQTL and EC GWAS datasets, suggesting that this locus may influence EC susceptibility through modulation of BRI3 expression. 41–44

First, the study relies on publicly available scRNA-seq and GWAS datasets, which, as secondary data sources, may be subject to selection bias, limited sample diversity, and inconsistencies in data quality, and the lack of prospective cohort studies limits the ability to capture temporal relationships and disease progression. Second, although MR provides statistical inference of causality, the observed associations, including immune alterations in CD56^dim NK cells, cannot definitively distinguish between causal mediation and mere correlation, and further experimental validation is required to clarify the directionality and functional impact of these immune changes. Finally, unmeasured genetic or environmental factors may influence the results, so these findings should be interpreted as exploratory and hypothesis-generating, providing a basis for future longitudinal and functional studies.