Identification of molecular subtypes and immune infiltration in endometriosis: a novel bioinformatics analysis and In vitro validation

Introduction Endometriosis is a worldwide gynacological diseases, affecting in 6–10% of women of reproductive age. The aim of this study was to investigate the gene network and potential signatures of immune infiltration in endometriosis. Methods The expression profiles of GSE51981, GSE6364, and GSE7305 were obtained from the Gene Expression Omnibus (GEO) database. Core modules and central genes related to immune characteristics were identified using a weighted gene coexpression network analysis. Bioinformatics analysis was performed to identify central genes in immune infiltration. Protein-protein interaction (PPI) network was used to identify the hub genes. We then constructed subtypes of endometriosis samples and calculated their correlation with hub genes. qRTPCR and Western blotting were used to verify our findings. Results We identified 10 candidate hub genes (GZMB, PRF1, KIR2DL1, KIR2DL3, KIR3DL1, KIR2DL4, FGB, IGFBP1, RBP4, and PROK1) that were significantly correlated with immune infiltration. Our study established a detailed immune network and systematically elucidated the molecular mechanism underlying endometriosis from the aspect of immune infiltration. Discussion Our study provides comprehensive insights into the immunology involved in endometriosis and might contribute to the development of immunotherapy for endometriosis. Furthermore, our study sheds light on the underlying molecular mechanism of endometriosis and might help improve the diagnosis and treatment of this condition.


Introduction
Endometriosis is a chronic gynecological disorder characterized by the presence and infiltration of ectopic endometrial tissue (1,2).It affects nearly 5% to 10% of women of reproductive age and could cause pelvic pain (50% to 80%) and infertility (up to 50%) (3)(4)(5)(6).Other common symptoms include dysmenorrhea, dyspareunia, dysuria, dyschezia, and fatigue (2,4).Although various biomarkers such as IL-2, anti-PEP, CA125, and miR-150-5p have been studied as clinical diagnostic tools, there is generally no single indicator that can directly diagnose endometriosis (7).Surgery remains the gold standard for diagnosis (8), but it has complications, and early lesions might be missed during operations (9).Therefore, a thorough understanding of the molecular mechanisms underlying endometriosis and exploration of a non-invasive diagnostic indicator with high specificity is essential.
Sampson's retrograde stigma theory is the most widely accepted theory on the pathogenesis of endometriosis (6).However, endometriosis is now recognized as a systemic disease that can spread through lymphatic and hematogenous metastases (10), not just limited to pelvic lesions (6).In addition, recent research suggests that endometriosis is closely related to inflammation and autoimmunity (11,12).
Nowadays, the bioinformatics approach is widely used to reveal the molecular mechanisms associated with endometriosis.The underlying pathogenic factors are variable and could be related to the lesion microenvironment, individual differences, and environmental factors (4).As mentioned earlier, endometriosis is closely related to immunity, and some studies have revealed a part of the mechanisms.For instance, Wu et al. (13) reported that an increase in CD8 + T cells and a decrease in CD163 + macrophages might create a pro-inflammatory endometrial immune environment, leading to endometriosis.Zhong et al. (14) found that M2 macrophages significantly increase during the progression of endometriosis.In a previous study, it was discovered that PDLIM3, a specific biomarker in endometriosis, was correlated with multiple immune cells, such as M2 macrophages, activated NK cells, and CD8 + T cells (15).However, many unknowns still exist regarding the immune mechanisms of endometriosis.
In our study, we downloaded three endometriosis datasets (GSE51981, GSE6364, and GSE7305) from the GEO database to evaluate immune cell infiltration and identify immune correlation with endometriosis for further analysis.We identified endometriosis subtypes, explored differentially expressed genes (DEGs), and constructed a co-expression network.Central genes were identified, and subsequently, hub genes were comprehensively analyzed, with a particular focus on their immunological characteristics.Our study established a detailed immune network and aimed to systematically clarify the underlying molecular mechanisms of endometriosis, particularly from the perspective of immune infiltration.The results of our study could help facilitate the development of immune therapeutic approaches for endometriosis.

Data collection and preprocessing
To identify DEGs, the expression profiles of GSE51981, GSE6364, and GSE7305 were obtained from the GEO database using the GPL570[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.We used the GEOquery ( 16) package of R software with the keyword "endometriosis".GSE5198 included 77 endometriosis samples, GSE6364 included 21 endometriosis samples, and GSE7305 included 10 endometriosis samples.We selected and combined the endometriosis samples from the three datasets, which underwent removal of batch effects, standardization and annotation of probes, and other data cleaning processes.We used the "removeBatchEffect" function to remove batch effects using the "limma" R package and employed "normalizeBetweenArrays" function to standardize the data.The overall analytical flow diagram is shown in Figure 1.

Evaluation and analysis of immune cell infiltration and correlation between immune cells
We used the CIBERSORT algorithm (17) to convert the normalized gene expression matrix into a matrix of 22 types of immune cells.We uploaded the gene expression matrix data to CIBERSORT and filtered out samples with P < 0.05 to obtain the immune cell infiltration matrix.A total of 108 samples were uploaded, and 32 samples were filtered out based on the criteria of P < 0.05.The Ggplot2 (18) package was used to draw bar charts showing the distribution of the 22 immune cell infiltrates in each sample.The Corrplot package was used to plot the correlation heatmap to visualize the correlations of these 22 immune cell infiltrates.The Ggplot2 package was used to visualize the expression of immune-related HLA family and KIR family genes.

Gene set variation analysis
GSVA ( 19) is a non-parametric, unsupervised analysis that calculates enrichment scores for specific sets of genes in each sample.GSVA contributes to the construction of pathway-central models of biology and meets the need for bioinformatic methods for RNA-seq data (20).In our study, we used the expression profile data of endometriosis patients based on the c2.cp.v7.4.symbols.gmtdataset to analyze the corresponding biological characteristics and observe changes in immune-related pathways to explore the underlying mechanism of pathogenesis in endometriosis.The c2.cp.v7.4.gm geneset is a collection of gene sets in the Molecular Signatures Database (MSigDB) version 7.4.It contains coexpression-based gene sets that capture potential functional relationships among genes across diverse tissues or cell types (https://docs.gsea-msigdb.org/#MSigDB/Release_Notes/MSigDB_7.4/#_top).

Construction and analysis of related molecular subtypes
Consensus clustering involves the identification of molecularly related subtypes using the Consensus Cluster Plus algorithm (21).The consensus matrix is determined by consensus clustering to classify the samples, with the k value of the cluster number set between 2 and 9.When the cumulative distribution function index reaches the approximate maximum value, the optimal K value is determined.The classification is then validated by principal component analysis (PCA) of mRNA expression profiles.
The construction of immune cell-related subtypes was based on the consensus clustering of immune cell infiltration data.The results were used to define Cluster 1 and Cluster 2. DEGs were identified by limma (22) with the criteria of P-value < 0.05 and |log2FC| > 0.3.We used the ggplot2 package to draw a volcano plot and a heatmap of the DEGs to visualize their expression.

Weighted gene co-expression network based on immune cell subsets
We used WGCNA based on immune cell subsets (23) to process data.We used the c2.cp.v7.4.symbols.gmtdataset, which includes immune-related pathways and genes selected from the dataset obtained through GSVA.First, the soft threshold value of network construction was selected, and the adjacency matrix was a continuous value between 0 and 1, ensuring that the constructed network followed a power-law distribution and was closer to the real biological network state.Second, a scale-free network was constructed using the function of block modules, and the co-expression modules were identified by block partitioning analysis, grouping genes with similar expression patterns.All modules were summarized by modular characteristic genes, which were the most important major components of each module and were defined as synthetic genes representing the expression profile of all genes in a given module.Flow diagram of methodologies applied to explore the biological characteristics of endometriosis.

Functional enrichment analysis
The ClusterProfiler package (24) was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs, and the common genes were analyzed using WGCNA.The critical value of false discovery rate < 0.05 was considered statistically significant.The GO and KEGG results were visualized using the pathview package (25), and the first two KEGG enrichment results were evaluated.
To investigate the differences in biological processes among different subgroups, we used gene set enrichment analysis (GSEA) based on the gene expression profile dataset of patients with endometriosis.The C2.cp.v7.4.symbols.gmtgene set was downloaded from MSIGDB for GSEA, and a P-value < 0.05 was considered statistically significant.

PPI network construction
We used the STRING database (26) (http://string-db.org,version 11.09) online tool to obtain PPI data.Candidate targets were input into STRING, and genes with scores greater than 0.3 were selected to construct a network model using Cytoscape (V3.7.2) (27).Cytohubba plug-in (28) was used to screen the top 10 hub genes based on the score.Pearson correlation analysis was conducted between the hub genes and immune cells, and both cor > 0.55 and > 0.7 were tried.An interaction score of 0.55 and P < 0.05 was finally set as the significance criteria for further analysis.The correlation between hub genes and immune checkpoint genes was assessed as previously using corrgram (29).

Quantitative real-time PCR
Total RNA was isolated from 6 ovarian endometriosis tissue samples and 9 eutopic endometrial tissue samples from endometriosis-free women (normal control, NC) using Trizol (RNAiso Plus; Takara, Japan) according to the manufacturer's instructions.One microgram of RNA was reverse-transcribed into cDNA using a PrimeScript ™ RT reagent kit (Takara) according to the manufacturer's protocol.Amplification was performed on a quantitative real-time (qRT)-PCR device (QuantStudio5; Thermo Fisher Scientific, Waltham, MA, USA) using TB Green ® Premix EX TaqTM II (Takara) and gene-specific primers (Sangon, Shanghai, China).b-Actin was used as the internal control.The relative expression levels of the genes were calculated using the 2 -DDCT method.

Protein extraction and western blotting
Ovarian endometriosis tissues and normal eutopic endometrium tissues were lysed in RIPA buffer containing protease and phosphatase inhibitor cocktail.After centrifugation at 12,000 rpm and 4°C for 15 min, the protein concentration in the supernatant was determined using a BCA protein assay.Equal amounts of total protein (30 µg) were loaded and separated on 12.5% SDS-PAGE and then transferred to PVDF membranes.After blocking with protein-free rapid blocking buffer for 0.5 h at room temperature, the membranes were incubated with primary antibodies overnight at 4°C.Rabbit anti-RBP4 antibody (1:5000; Abcam) and rabbit anti-b-actin (1:1000, cell signaling technology) were used as primary antibodies, with b-actin serving as the internal loading control.Then, the membranes were then incubated with secondary antibodies (1:3000) for 1 h at room temperature.Excess secondary antibody was removed by three washes in TBST.The targeted protein bands were visualized and imaged using an ECL Western blotting kit (NCM Biotech, China).

Statistical analysis
All data processing and analysis were conducted using R software (version 4.1.0).The statistical significance of the normally distributed continuous variables was estimated using the independent Student t test, while the differences between nonnormally distributed continuous variables were analyzed using the Mann-Whitney U test (Wilcoxon rank-sum test).The Chi-square test or Fisher's exact test was used to compare and analyze the statistical significance between categorical variables in two groups.All statistical P values were bilateral, and P < 0.05 was considered statistically significant.

Evaluation of immune cells infiltration and immune correlation analysis in endometriosis
To evaluate the extent of immune infiltration in patients with endometriosis, we used the CIBERSORT algorithm to analyze the proportions of the immune subset in endometriosis.We constructed immune cell maps (Figure 2A) and correlation maps of these 22 types of immune cells (Figure 2B) in endometriosis samples.The 22 types of immune cells are listed in Supplementary Table 1.We used GSVA to calculate the enrichment scores for all genes in each immune cell type.We observed changes in specific gene sets associated with 10 immune-related pathways in the c2.cp.v7.4.symbols.gmtdataset, including Toll-like receptor signaling pathway, Nod-like receptor signaling pathway, Rigi-like receptor signaling pathway, Jak-stat signaling pathway, NK cellmediated cytotoxicity, T cell receptor signaling pathway, B cell receptor signaling pathway, autoimmune thyroid disease, allograft rejection, and graft versus host disease the gene sets and (Figure 2C).The HLA and KIR family genes were immunerelated gene families, and the expression of each gene in the family showed different manifestations.HLA-A, HLA-B, and HLA-C showed relatively high expression, while the expression of the KIR family was relatively low (Figure 2D).

Identification of two endometriosis subtypes and exploration of DEGs based on immune characteristics
To determine the biological differences between different immune subtypes in endometriosis, we classified the endometriosis samples into cluster 1 and cluster 2 based on the consensus clustering of immune cell infiltration data.The two subtypes were clearly distinguished (Figure 3A).Differential gene expression was identified with limma (P-value < 0.05 and |log2FC| > 0.3).A total of 140 DEGs were obtained, including 29 low-expression genes and 111 high-expression genes (Supplementary Table 7).The volcano plot indicated a higher number of genes with upregulated expression compared to genes with downregulated expression (Figure 3B).In the heat map, the genes in cluster 1 exhibited an overall downregulated trend, while the genes in cluster 2 showed an overall upregulated trend (Figure 3C).GSEA and GSVA were applied for all the DEGs between cluster 1 and cluster 2, and the results were shown in Supplementary Tables 2, 3, respectively.The results of GSEA were visualized in Supplementary Figure 1A, indicating that immunoregulatory interactions and the complement system play crucial roles in the pathology of endometriosis.The results of GSVA (Supplementary Figure 1C) revealed that the B cell receptor signaling pathway and CD8 TCR downstream pathway are important pathways involved in endometriosis.These findings indicate that the immunology system plays an essential role in the pathologies of endometriosis.Frontiers in Immunology frontiersin.org

Co-expression network construction and central genes identification
To explore the central genes involved in immune infiltration in endometriosis, we constructed a gene co-expression network using the WGCNA algorithm.The gene expression profile of the endometriosis immune subtype was evaluated using cluster analysis (Figure 4A).To ensure that the network was scale-free, we selected a soft threshold of b = 2 (Figure 4B).We then transformed the representation matrix into an adjacency matrix and a topological matrix and used the average linkage hierarchical clustering method to cluster the genes.We set a minimum number of 50 genes in each gene network module, based on the standard of the hybrid dynamic clipping tree.We used the dynamic shearing method to determine the gene module and calculated the characteristic gene values for each module.
We then conducted Cluster analysis on the modules and merged modules that were close to each other using the following parameters: height = 0.25, deep split = 4, and minimum module = 50, resulting in a total of 11 modules (Figure 4C).The module with the highest correlation with immune characteristics was the green module (r = 0.6, P = 5e-12; Figure 4D).We identified 285 genes in the green module.We then compared these genes to the 140 DEGs based on immune characteristics and identified 52 common genes (Figure 4E).Frontiers in Immunology frontiersin.org

GO and KEGG enrichment analysis of cerntral genes in endometriosis
To further investigate the functions of the 52 common genes, we performed GO and KEGG analyses.We used a cutoff p value of 0.5 for these analyses.GO analysis showed that the common genes were mainly related to leukocyte-mediated cytotoxicity, negative regulation of growth, NK cell-mediated cytotoxicity, NK cellmediated immunity, and positive regulation of protein transport (Figures 5A, C; Supplementary Table 4).KEGG analysis showed that the pathway was mainly related to graft-versus-host disease, NK cell-mediated cytotoxicity, antigen processing and presentation, complement and coagulation cascades, allograft rejection, and NK cells in both resting and activated states (Figures 5B, D; Supplementary Table 5).Pathway enrichment analysis also showed that the common genes were significantly expressed in the complement and coagulation cascade (Figure 5E) and NK cellmediated cytotoxicity (Figure 5F) pathways.Frontiers in Immunology frontiersin.org

Pathway enrichment of central genes in endometriosis
To address potential biases in the enrichment analysis of intersection genes, we performed GSEA on all 140 DEGs related to immune subtypes in the endometriosis expression profile.We used the C2.cp.v7.4.Symbols.gmtreference gene set for this analysis.The results showed that DEGs were significantly enriched in the chemokine signaling pathway, cytokine cytokine receptor interaction, NK cell mediated cytotoxicity, cell cycle checkpoints, cell cycle mitotic, homology directed repair, and reactome mitotic prometaphase (Figures 6A-G).We also performed GSVA using the C2.cp.v7.4.symbols.gmtreference gene set to explore the potential mechanism underlying the pathogenesis of endometriosis.We observed significant differences in 10 gene sets based on immune subtypes, including innate immune system, biocarta csk pathway, graft versus host disease, WP oxidative damage, neutrophil degranulation, signaling by CSF3 G CSF, WP microglia pathogen phagocytosis pathway, WP chemokine signaling pathway, WP IL3 signaling pathway, and interferon-gamma signaling.Most of these DEG sets were immunerelated pathways (Figure 6H).Frontiers in Immunology frontiersin.org09

Construction of PPI network
To explore the protein functions of the 52 candidate genes (Supplementary Table 6), we constructed a PPI network using the STRING database and Cytoscape software.With a confidence score of 0.55, we identified 45 genes that showed close interactions with each other (Supplementary Figure 2A), while only 8 genes remained with a cut-off of > 0.7 (Supplementary Figure 3).We utilized an interaction score of 0.55 to identify the hub genes.Using the cytoHubba plugin, we calculated the value of the PPI network and selected the top 10 hub genes (Supplementary Figure 2B) and their extensions: GZMB, PRF1, KIR2DL1, KIR2DL3, KIR3DL1, KIR2DL4, FGB, IGFBP1, RBP4, and PROK1 (Supplementary Figure 2C).To further explore the upstream regulatory relationships, we predicted miRNA interactions with the hub gene by selecting the results of "Functional MTI" and "positive" evidence, based on the multiMiR and TarBase v.8 databases.We found that 8 mRNAs in the hub genes interacted with 11 miRNAs (Supplementary Figure 2D).

Identification of three endometriosis subtypes through unsupervised clustering based on hub genes
To investigate the biological differences among different subtypes based on the characteristics of endometriosis, we used the ConsensusClusterPlus software package to construct subtypes based on the expression profile of the 10 hub genes.When k = 3, the classification was reliable and stable (Figures 7A-C).The samples were divided into cluster 1, cluster 2, and cluster 3. PCA confirmed that cluster 1, cluster 2, and cluster 3 showed significant differences (Figure 7D).We also analyzed the correlation of immune cell fractions by disease subtypes and found that the three components showed different expressions in different immune cells.NK.cells.activated,T.cells.CD8, and NK.cells.restingshowed a large difference in expression (Supplementary Figure 4G).For the hub genes, we calculated the correlation between the hub gene expression and immune cell fraction (Pearson's coefficient).Most hub genes were significantly correlated with immune markers.We Frontiers in Immunology frontiersin.orgobserved a significant positive correlation between GZMB and NK.cells.activated(P = 1.23e-11,R = 0.59) (Supplementary Figure 4A), KIR2DL1 and NK.cells.resting(P = 2.51e-8, R = 0.50) (Supplementary Figure 4B), KIR2DL3 and NK.cells.resting(P = 1.55e-11,R = 0.59) (Supplementary Figure 4C), KIR2DL4 and NK.cells.resting(P = 4.26e-12, R = 0.60) (Supplementary Figure 4D), PRF1 and NK.cells.activated(P = 3.54e-10, R = 0.56) (Supplementary Figure 4E), and PRF1 and NK.cells.resting(P = 1.15e-11,R = 0.59) (Supplementary Figure 4F).Further analysis of the correlation between hub genes and immune checkpoint genes revealed that hub genes were strongly positively correlated with CD40, IDO1, LAG3, TNF, and TNFRSF18 immune checkpoint genes (Supplementary Figure 4H).These results showed that the expression of hub genes was significantly positively correlated with immune characteristics.

Expression validation of hub genes
To validate the transcriptional and protein expression of the hub genes, we used qRT-PCR and Western blot analyses to detect their expression in both ovarian endometriosis tissue and normal eutopic endometrial tissue from the same patients (Figure 8).We found that, except for FGB, KIR2DL1, and KIR2DL3, all other hub genes showed significantly different expression levels in the endometriosis and normal control groups (Figure 8A).Specifically, the expression of KIR2DL4 was remarkably downregulated, while the expressions of PROK1, IGFBP1, RBP4, G2MB, KIR3DL1, and PRF1 were significantly upregulated in ovarian endometriosis.Moreover, the protein expression of RBP4 was remarkably higher in ovarian endometriosis than in normal control (Figure 8B).

Discussion
Endometriosis is a common disease that affects up to 10% of women.It is characterized by pain and infertility (31).Up to 90% of cases of chronic pelvic pain in women of reproductive age are associated with endometriosis (32).Endometriosis is also a risk factor for ovarian cancer (33).However, endometriosis is not curable, and its underlying etiology remains unclear.In a previous study, we identified PDLIM3 as a specific biomarker for endometriosis that was associated with multiple immune cells (15).However, that analysis focused on the relationship between a single gene and the immune cell environment.
In the present study, we conducted a more comprehensive analysis of the immune cell environment and explored immune- Frontiers in Immunology frontiersin.orgassociated genes in endometriosis.We merged three mRNA microarray datasets (GSE51981, GSE6364, GSE7305) into one dataset with 108 samples of endometriosis and conducted several analyses.Our findings indicated that hub genes, enriched modules, and pathways have genetic effects on endometriosis.Moreover, the predicted genes in the merged dataset might interact with each other and coregulate endometriosis.HLA and KIR are immune-related gene families.KIRs are expressed by NK cells and can recognize class I HLA on target cells.HLA-C is a ligand of KIRs and is an essential regulator of NK cell activity (34,35).Modulation of the HLA-KIR axis provides new prospects for cancer treatment (36) and immunotherapy of other diseases such as COPD (37).A previous study revealed that the expression of HLA-C*03:03*01 was increased in endometriosis.Aberrant expression of KIR2DL1 impaired NK cell cytotoxicity in endometriosis (38).KIR2DS5-positive women with endometriosis had a lesser chance of peritoneal disease, and KIR2DS5 might be a protective factor for endometriosis (34).In this study, we constructed 22 immune profiles from endometriosis samples.GSVA was used to calculate enrichment scores for specific gene sets in each sample, and variations of 10 immune pathways, including Toll-like receptor signaling pathway, Nod-like receptor signaling pathway, and Rig I-like receptor signaling pathway, were observed.HLA-A, HLA-B, and HLA-C expression levels were relatively high, and KIR expression was relatively low, which is consistent with previous reports (34,35).
To investigate the biological functions of the differentially expressed genes, we analyzed 52 common genes that were common to both immunogrouping and WGCNA analyses and found several functions were associated with endometriosis, including leukocyte-mediated cytotoxicity, negative regulation of growth, NK cell-mediated cytotoxicity, NK cell-mediated immunity, and positive regulation of protein transport.Further analysis using GO and KEGG pathways showed that the complement and coagulation cascade pathways and NK cellmediated cytotoxicity pathways were significantly associated with endometriosis.Consistent with our findings, the abnormal immune function of NK cells is closely related to endometriosis (39).In particular, impaired NK cell cytotoxic activity has been linked to the development of endometriosis (40,41).To validate our results, we performed GSEA and GSVA analysis on differential gene sets related to immune subtypes.Our analysis showed common enrichment in NK cell-mediated cytotoxicity and graft versus host disease, and most differential gene sets were immune-related pathways.These results further support the important role of NK cells in the development of endometriosis.
The immune system plays a crucial role in the pathogenesis of endometriosis, with immune cells infiltrating the ectopic endometrial lesions.Immune infiltration points in endometriosis include the involvement of various immune cell types, such as macrophages, neutrophils, NK cells, and CD8 + T cells.Macrophages, found in increased numbers in the peritoneal fluid and at the site of the lesions, contribute to the development and maintenance of endometriotic lesions by producing growth factors, angiogenic factors, and pro-inflammatory cytokines (42)(43)(44).Previous study suggests that M2 macrophages may play a role in the development and recurrence of endometriosis (45).Additionally, microvesicles secreted by M1 macrophages in endometriosis can induce polarization of M2 macrophages towards an M1 phenotype, potentially inhibiting the development of endometriosis (46).These studies indicate the presence of both M1 and M2 macrophages in endometriotic lesions and suggest that repolarizing M2 macrophages towards an M1 phenotype could be a potential therapeutic strategy.However, further research and clinical studies are needed to fully understand the role of macrophage subtypes and their potential for the treatment of endometriosis.Similarly, neutrophils are present in the peritoneal fluid of women with endometriosis and promote angiogenesis and tissue remodeling in the formation of endometriotic lesions (47, 48).Moreover, NK cells in peripheral blood might help to destroy retrograde endometrial cells, thus preventing endometriosis development.The reduction in NK cell cytotoxic activity might result in endometriosis (49).The alterations in the immunological parameters of NK cells not only included lower NK cell cytotoxicity but also involved a shift in the balance between type 1 and type 2 NK cells, as well as changes in the percentages of inhibitory and activating NK cell receptors (39).Thus, abnormalities in the function of NK cells, such as reduced cytotoxicity, altered activity and phenotype, and imbalances in receptor expression, might contribute to the development and pathology of endometriosis.In our study, the NK cell-mediated cytotoxicity pathway is implicated in several analyses, aligning with previous findings.CD8 T cells also play a role in the development and progression of endometriosis.An increased number of CD8 T cells has been observed in the peritoneal fluid of women with endometriosis compared to those without the disease, suggesting a possible role of CD8 T cells in the local immune response (50).CD8 T cells in the peritoneal fluid of endometriosis patients exhibited a reduced cytotoxic capacity, potentially contributing to impaired immune surveillance and facilitating the establishment of endometrial implants in the peritoneal cavity (51).
To further elucidate the pathogenesis of endometriosis, we performed PPI analyses to identify the endometriosis-associated hub genes.We selected 10 DEGs as hub genes, namely, GZMB, PRF1, KIR2DL1, KIR2DL3, KIR3DL1, KIR2DL4, FGB, IGFBP1, RBP4, and PROK1.Among these, KIR2DL1, KIR2DL3, KIR3DL1, and KIR2DL4 belong to NK cell inhibitory receptors (39).Their functions are mentioned above.GZMB is a component of cytolytic granules within NK cells (52) and was reported to be correlated with cervical cancer (53).RBP4 belongs to the lipocalin family and is the major transport protein of the hydrophobic molecule retinol, which is known as vitamin A (54).A recent study showed that RBP4 is involved in the pathological process of endometriosis and could promote the activity, proliferation, and invasion of ectopic cells (55).We validated that RBP4 was overexpressed in endometriosis patients.In addition, there are no reports about RBP4 and other uterine pathologies, which indicates that this signature is unique to endometriosis.PROK1 is a secreted peptide that belongs to the prokineticin family and performs a wide range of functions in angiogenesis, modulation of inflammatory responses, and regulation of hematopoiesis (56).Tiberi et al. discovered that the expression of PROK1 in normal female tissues was much higher than that in women with endometriosis, suggesting that women with endometriosis might show abnormal vascular function (57).
In our study, all the hub genes were found to be highly correlated with endometriosis.For example, GZMB (Granzyme B) belongs to the granzyme family of serine proteases that are found in the cytotoxic granules of cytotoxic T lymphocytes and NK cells (58).Granzymes play a critical role in the immune system by inducing apoptosis in target cells, such as virally infected cells or tumor cells (59).GZMB is a therapeutic target for endometriosis, based on its role in immune response, tissue remodeling, and angiogenesis.Targeting GZMB might alleviate endometriosisassociated pain and inflammation while preserving fertility (60).Apart from their association with endometriosis, some of the hub genes have also been linked to other uterine pathologies.For instance, the KIR family (61), IGFBP1 (62), and PROK1 (63) have been correlated with recurrent implantation failure, while GZMB (64) and FGB (65) have been implicated in cervical carcinoma.However, there is no report about the relationship between PRF1 and RBP4 and other uterine pathologies, which suggests that these two might be unique signatures of endometriosis.
Our study has some limitations.Firstly, as the pathogenesis of endometriosis is multidimensional, discussing it only from the aspect of immune infiltration might not be comprehensive enough.Secondly, all study data were obtained from publicly available databases, and additional clinical characteristics of endometriosis patients should be included in subgroup analysis.Thirdly, although we validated the expression of hub genes through qRT-PCR and western blot, further multicenter and prospective studies are required to evaluate the possible applications of molecular signatures.Additionally, more in vivo and in vitro experiments are needed to elucidate the molecular mechanisms of hub genes for clinical applications.Furthermore, while we have analyzed gene expression data for HLA and KIR genes, we acknowledge that further investigation in specific cell types and disease contexts through sequencing or genotyping approaches is necessary to better understand the heterogeneity of these gene families.Lastly, scRNA-seq analysis may provide more valuable insights into immune cell heterogeneity compared to bulk RNAseq, and we plan to conduct scRNA-seq analysis in future studies.
This study presents several novel aspects compared to previous similar research.Firstly, we incorporate a novel approach by considering the immune phenotype of patients with endometriosis.By analyzing the immune cell composition and immune-related pathways in endometriosis, we provide valuable insights into immune dysregulation associated with the disease.Secondly, we utilize diverse datasets to comprehensively analyze the immune landscape, revealing new relationships between immune cell subgroups, gene expression patterns, and immune pathways in endometriosis.Advanced analysis methods, such as CIBERSORT, GSVA, and WGCNA, further deepened the exploration, identifying key genes associated with immune subtypes.Notably, the study uncovers two previously unknown genes, KIR2DL3 and FGB, in endometriosis, contributing to the understanding of potential genetic factors in the disease.
In conclusion, this study provides comprehensive and reliable evidence for the pathogenesis of endometriosis from the perspective of immune infiltration.The results showed that most of the DEGs are involved in immune-related pathways.Ten hub genes, which show significant correlation with immune markers, might be potential diagnostic and therapeutic targets for endometriosis.
Our study provided a reliable and deep analysis of the mechanism underlying endometriosis from the perspective of immune infiltration.A total of 10 hub genes were identified and were validated by qRT-PCR and western blotting.The majority of hub genes showed significantly different expression between endometriosis patients and normal controls.Of these genes, PROK1 (Prokineticin 1) has been suggested to be associated with endometriosis due to its role in promoting angiogenesis (71).Altered expression of IGFBP1 (Insulin-like growth factor-binding protein 1) has been found in the endometrium of women with endometriosis, suggesting a potential association (72).KIR2DL4 and KIR3DL1 are members of the killer cell immunoglobulin-like receptor family and might be associated with the pathogenesis of endometriosis through their involvement in the regulation of NK cell activity (35).No direct link has been established between RBP4, G2MB, or PRF1 and endometriosis, necessitating further investigation.institutional requirements.The participants provided their written informed consent to participate in this study.

2
FIGURE 2 Evaluation of immune cell infiltration and immune correlation analysis.(A) Barplot shows the proportion of 22 types of immune cells in endometriosis samples.The column of the graph is sample.(B) Correlation heat map of 22 immune cell infiltrates; Blue represents positive, red represents negative.The depth of the color indicates the strength of the correlation.(C) Changes in enrichment of 10 immune-related pathways.(D) Expression differences of HLA and KIR family genes.

3
FIGURE 3 Differential analysis based on immune cell subtypes.(A) PCA analysis based on immune cell subtypes, cluster1 is green, cluster 2 is yellow.(B) Volcano map of differential analysis based on immune cell subtype, blue dots represent low expression and red dots represent high expression.(C) Heat map based on differential analysis of immune cell subsets, cluster1 is green, cluster 2 is red; blue dots represent low expression and red dots represent high expression.PCA, principal component analysis.

4
FIGURE 4 Construction and module analysis of the weighted co-expression network.(A) Sample clustering tree based on Euclidean distance.(B) Network topology analysis under various soft threshold power settings.Left: The X-axis represents the soft threshold power, and the Y-axis represents the fitting index of scale-free topological model.Right: The X-axis represents the soft threshold power, and the Y-axis reflects the average connectivity (degree).(C) Clustering tree of genes with different similarities based on topological overlap, and assigned module colors.(D) Association with Module-trait.Each row corresponds to a module, and each column corresponds to a feature.Each cell contains the corresponding correlation and P values.This table is color-coded according to the correlation of the color legends.(E) The intersection Venn diagram of core modules and central genes of WGCNA.WGCNA, weighted gene co-expression network based on immune cell subsets.

5
FIGURE 5 Central genes in GO, KEGG, and PATHWAY enrichment analysis.(A) Dot plot of central genes GO analysis.(B) Dot plot of central genes KEGG analysis.(C) Loop diagram of central genes GO analysis.(D) Loop diagram of central genes KEGG analysis.(E) Enrichment diagram of central genes in the complement and coagulation cascade pathways.(F) Enrichment diagram of central genes in the natural killer cell-mediated cytotoxicity pathway.GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

7
FIGURE 7 Characteristics based on endometriosis-related molecular typing.(A) Consensus CDF.(B) Delta region.(C) Consistent matrix at k = 3.The rows and columns of the matrix represent samples.(D) The PCA diagram verifies the stability and reliability of classification.PCA, principal component analysis.

8
FIGURE 8Expression validation of hub genes.(A) qRT-PCR was performed to determine mRNA expression of hub genes in normal eutopic endometrial tissue and ovarian endometriosis tissue.(B) Western blotting was performed to determine the protein expression of RBP4 in normal eutopic endometrial tissue and ovarian endometriosis tissue.EMS, endometriosis.Data are expressed as the mean ± SEM. *P < 0.05, **P < 0.01; ns, non-significant.