Bioinformatical analysis identifies PDLIM3 as a potential biomarker associated with immune infiltration in patients with endometriosis

Gene expression profile

We searched the associated gene expression profiles of endometriosis in NCBI-GEO (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2005), using the keywords “endometriosis” and “Homo sapiens”. Three gene expression datasets containing microarray data from normal endometrial tissues and ectopic endometrial tissues, GSE11691 (Hull et al., 2008), GSE23339 (Hawkins et al., 2011) and GSE5108 (Eyster et al., 2007), were selected for further analysis. The GSE11691 microarray data consisted of nine normal human endometrium and nine ectopic human endometrium (two from the peritoneal wall, seven from the broad ligament) and was generated using the GPL96 platform Affymetrix Human Genome U133A Array. The GSE23339 microarray data consisted of nine normal human endometrium and 10 ectopic human endometrium (all from the endometriotic cysts of the ovaries) and was generated using the GPL6102 Illumina human-6 v2.0 expression beadchip. The GSE5108 microarray data consisted of 11 normal human endometrium and 11 ectopic human endometrium (five from the peritoneal wall, six from the endometriotic cysts of the ovaries) and was generated using the GPL2895 GE Healthcare/Amersham Biosciences CodeLink Human Whole Genome Bioarray. In addition, the GSE120103 (Bhat et al., 2019) dataset, containing 18 endometriosis (all from the endometriotic cysts of the ovaries) and 18 control samples, was used as the validation cohort using the GPL6480 platform Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F. GSE12768 (Ritchie et al., 2015) dataset, containing two endometriosis (both from the endometriotic cysts of the ovaries) and two control samples, was used as another validation cohort using the GPL7304 Institut Cochin HG18 60mer expression array 47K.

Identification of differentially expressed genes (DEGs)

The raw data of three datasets, GSE11691, GSE23339 and GSE5108, was merged and removed batch effects by “sva” package (http://www.bioconductor.org/) (Ritchie et al., 2015). Then, the “limma” package (http://www.bioconductor.org/) was used to attain DEGs between endometriosis samples and normal samples. The raw data of GSE12768 was analyzed separately to acquired the DEGs list, and compared with the list of three datasets in order to verify the reliability of our analysis. The significant DEGs were identified according to adjusted P value < 0.05 and ∣log fold change (FC)∣ > 2.

Functional and pathway enrichment analysis of DEGs

To reflect gene functions, gene ontology (GO) (http://geneontology.org/) (Thomas, 2017) biological processes were shown in three aspects: biological processes (BP), cellular component (CC) and molecular function (MF). And the analysis process was achieved by “clusterProfiler” package (Wu et al., 2021; Yu et al., 2012) in the present study, using the cut-off criterion P value < 0.05. ToppGene (https://toppgene.cchmc.org/enrichment.jsp) (Chen et al., 2009) is a free tool for uncovering the biological significance behind a cluster of genes. The REACTOME (https://reactome.org/) (Jassal et al., 2020) pathway enrichment analysis, was exactly performed using this tool. A P value < 0.05 was considered as the threshold points.

Moreover, to get more information of mechanism behind DEGs, gene set enrichment analysis (GSEA) was employed to explore the potential feature between the endometriosis and control groups. The “c2.cp.kegg.v7.4.symbols.gmt” from the Molecular Signatures Database (MSigDB) was downloaded and used as the reference gene set. An adjusted P < 0.05 and false discovery rate (FDR) < 0.05 were considered as significant.

Protein–protein interaction (PPI) networks construction

STRING (https://string-db.org/) (Szklarczyk et al., 2015) is a tool designed to achieve a comprehensive and objective global network and visualize PPI information. Proteins associated with DEGs were chosen based on information in STRING (confidence score > 0.4), and then PPI networks were constructed using Cytoscape software (http://cytoscape.org/) (Shannon et al., 2003).

Screening of specific gene biomarkers

The least absolute shrinkage and selection operator (LASSO), a regression analysis algorithm, was applied to evaluate the value of DEGs in this study. The “glmnet” package in R did the above analysis and returned candidate potential genes. And the expression levels of candidate biomarkers were further tested in the GSE120103 dataset.

Receiver operating characteristic (ROC) curve analysis

To test the predictive accuracy of the identified biomarkers, we generated ROC curves using “pROC” package (Robin et al., 2011) based on the candidate genes and compare the predictive value of candidate genes by the area under the ROC curve (AUC). And then, the predictive value of these genes was also validated in the GSE120103 dataset.

Relationship between identified genes and infiltrating immune cells in endometriosis

We used the CIBERSORT (https://cibersortx.stanford.edu/) (Newman et al., 2015; Robin et al., 2011) to calculate immune cell infiltrations of endometriosis and control tissues. Then, we applied “corrplot” and “vioplot” R package to visualize the results and explore correlation of infiltrating immune cells. In addition, the relation of the candidate gene markers with the levels of infiltrating immune cells was calculated and visualized by “ggpubr” package.

Quantitative real-time PCR (qRT-PCR)

Total RNA was isolated from ovarian endometriosis and eutopic endometrium (NC). These operations were approved by Ethics committee of Shanghai First Maternity and Infant Hospital (KS21198). We used Trizol (RNAiso Plus; Takara, Tokyo, Japan) to isolate RNA from these two kinds of tissues according to the manufacturer’s instructions. Then, one microgram of RNA was reversely transcribed into cDNA using a PrimeScript^™ RT reagent kit (Takara, Tokyo, Japan). Amplification was performed using gene-specific primers (Sangon, Shanghai, China) and TB Green^® Premix EX TaqTM II (Takara, Tokyo, Japan) on a qRT-PCR device (QuantStudio5, Thermo Fisher Scientific, Waltham, MA, USA). β-actin was used as an internal control. The relative expression of the genes was calculated using the 2^–ΔΔCT method. The primers used were as follows: human β-actin-Forward: 5′-GTGCTATGTTGCTCTAGACTTCG-3′; human β-actin-Reverse: 5′-ATGCCACAGGATTCCATACC-3′; human PDLIM3-Forward: 5′-GAGTCGGACGTGTACCGGA-3′; human PDLIM3-Reverse: 5′-TGCCACTCCCACATTTGTCAC-3′.

Protein extraction and western blot

Normal endometrium tissues and ovarian endometriosis tissues were lysed in RIPA lysis buffer (PC101; EpiZyme, Shanghai, China) with protease and phosphatase inhibitor cocktail. The protein concentration was detected by a BCA protein assay. The same amounts of total protein (32 µg) were loaded and separated onto 12.5% SDS-PAGE and was then transferred to PVDF membranes. After blocking with protein-free rapid blocking buffer (PS108P; EpiZyme, Shanghai, China) for 0.5 h at room temperature, the membranes were incubated with primary antibodies at 4 °C overnight. The primary antibodies rabbit anti-PDLIM3 antibody (1:2,000; Servicebio) was used and rabbit anti-β-Actin (1:5,000; ABclonal) served as the internal control.

The membranes were probed with secondary antibodies (1:3,000, Servicebio) for 1.5 h at room temperature and three washes in TBST were processed to remove excess secondary antibody. We then used an ECL western blotting kit (NCM Biotech, Suzhou, China) to visualize and image the targeted protein bands.

Statistical analysis

Most statistical calculations were conducted using R (version 4.1.1). The results of qRT-PCR and western blot were shown as the mean ± SD and GraphPad Prism 8.0 (GraphPad Software, San Diego, CA, USA) were used for statistical analysis. All the comparisons between two groups were analyzed using Student’s t-test for normally distributed variables. If the variables distribution was abnormal, the Mann–Whitney U-test would be used. The relationship between the expression of biomarkers and infiltrating immune cells was analyzed using Spearman’s rank correlation analysis. Two-sided with P < 0.05 (95% confidence interval) and false discovery rate (FDR) < 0.05 were considered as statistically significant.

Identification of DEGs in endometriosis

Data from a total of 30 endometriosis samples (14 peritoneal endometriosis samples, 16 ovarian endometriosis samples) and 29 control endometrium samples from three GEO datasets (GSE11691, GSE23339 and GSE5108) were analyzed. The DEGs were calculated using the “limma” package after removing the batch effects. We obtained a total of 42 DEGs, including upregulated 23 genes and downregulated 19 genes (Table 1). All DEGs were shown in the heatmap and volcano map (Fig. 1). In addition, the DEGs list from GSE12768 dataset was displayed in Table S1.

GO and pathway enrichment of DEGs in endometriosis

The above list of DEGs was used to evaluate if they were enriched into specific GO or reactome pathways. GO analysis identified that the DEGs were significantly enriched in BP, including humoral immune response, complement activation, alternative pathway, antimicrobial humoral response, and so on (Table 2). In terms of MF, DEGs were mainly enriched in Toll-like receptor binding, G protein-coupled receptor binding, long-chain fatty acid binding, and chemokine receptor binding (Table 2). Moreover, CC demonstrated that filamentous actin, lipid droplet, sperm flagellum and actin filament were the most significantly enriched GO term (Table 2). The results were displayed in Figs. 2A and 2B. In addition, the GSEA results illustrated that the enriched pathways mainly involved cytokine–cytokine receptor interaction, systemic lupus erythematosus, chemokine signaling pathway, complement and coagulation cascades, and leishmania infection (Figs. 2C and 2D).

REACTOME pathway enrichment analysis was also used to obtain the signaling pathways for different genes. In this study, these DEGs were mainly involved in the regulation of complement cascade, suprofen pathway, indomethacin pathway and mefenamic acid pathway (Table 3).

PPI networks construction and validation of predictive feature biomarkers

We uploaded a list of proteins associated with DEGs, removed the disconnected nodes and then acquired a PPI network using Cytoscape software (Fig. 3A). There were 26 nodes and 50 edges on this network.

Furthermore, the LASSO algorithm was used to narrowed down the DEGs, leading to the identification of seven variables as specific biomarkers for endometriosis (Fig. 3B). The seven genes were FMO1, PDLIM3, FAM129A, CLDN11, C3, TMPRSS4 and DEFB1. And then, we used ROC curve to estimate the predictive value of genes. As shown in Fig. 3C, the ability of PDLIM3 in discriminating endometriosis from the control samples demonstrated a favorable diagnostic value, with an AUC of 0.955 (95% CI [0.894–0.997]). Moreover, this result was confirmed in the GSE120103 dataset with an AUC of 0.836 (95% CI [0.691–0.948]). Other six genes showed no significant difference in the GSE120103 dataset (Fig. 3D).

At the same time, in order to generate more accurate and reliable results, the GSE120103 dataset was employed to verify the expression levels of the seven features. The expression levels of PDLIM3 in endometriosis tissue were notably higher than those in the normal tissue (Fig. 4A). There was no significant difference between the two groups in other six gene expression (Fig. 4A). Therefore, PDLIM3 became one of the potential specific biomarkers in endometriosis.

QRT-PCR and western blot

To validate the expression of PDLIM3, qRT-PCR and western blot were used to measure transcriptional and protein expression in both ovarian endometriosis and normal eutopic endometrium tissue from endometriosis-free women. We found PDLIM3 had higher expression level in endometriosis group than normal control (qRT-PCR: n = 6 vs. 6; Western blot: n = 7 vs. 4) (Figs. 4B and 4C). The patients of endometriosis ranged in age from 26 to 43 years (average age 32 years). The BMI was 21.76 ± 2.42 (mean ± SD). And, their rASRM endometriosis stage was III or IV. Normal group samples were from non-endometriosis patients, ranged in age from 31 to 41 years (average age 34 years). The BMI was 21.23 ± 1.21. There was no significant different in their general characteristics.

Correlation analysis between PDLIM3 and infiltrating immune cells

The composition of immune cells in endometriosis tissues and normal control tissues were shown in Fig. 5A. The proportions of naive B cells (P < 0.001), plasma cells (P = 0.012), CD8+ T cells (P = 0.013), follicular helper T cells (P < 0.001), regulatory T cells (Tregs) (P = 0.004), resting NK cells (P = 0.037), and activated NK cells (P < 0.001) in endometriosis tissues were significantly lower than those in control tissues. However, the proportions of CD4+ memory resting T cells (P < 0.001), gamma delta T cells (P = 0.031), M1 Macrophages (P = 0.020), M2 Macrophages (P < 0.001), and neutrophils (P = 0.002) in endometriosis tissues were significantly higher than those in normal tissues (Fig. 5B). However, in GSE120103, only the proportion of Tregs cell (P = 0.037) in endometriosis tissues were significantly lower than that in control tissues (Fig. S1).

The association of 22 types of immune cells was calculated and drawn in Fig. 5C. Then, we found PDLIM3 was positively correlated with M2 Macrophages (r = 0.70, P < 0.001), neutrophils (r = 0.45, P < 0.001), CD4+ memory resting T cells (r = 0.44, P < 0.001), gamma delta T cells (r = 0.34, P = 0.0088), resting mast cells (r = 0.31, P = 0.019) and M1 Macrophages (r = 0.28, P = 0.037). However, PDLIM3 was negatively correlated with follicular helper T cells (r = –0.66, P < 0.001), activated NK cells (r = –0.46, P < 0.001), CD8+ T cells (r = –0.41, P = 0.0015), Tregs cells (r = –0.39, P = 0.0028), naive B cells (r = –0.37, P = 0.0047), plasma cells (r = –0.31, P = 0.018) and resting NK cells (r = –0.29, P = 0.031) (Fig. 6). Besides, PDLIM3 was significantly negatively correlated with activated mast cells (r = –0.69, P = 0.038) in GSE120103 (Fig. S1).