PharmMapper (http://www.lilab-ecust.cn/pharmmapper/) was used to identify molecular targets based on pharmacophore mapping approach (Liu et al., 2010). The pharmacophore database contains 7,000 models by extracting targets from TargetBank, DrugBank, BindingDB and PDTD. 2D structures of EGCG and ProEGCG molecules in SDF format were submitted to the server and the best mapping poses of the molecules against all the pharmacophore models in the database were identified. Default settings were used for other parameters; the proteins species was limited to Homo sapiens only. The fit score of EGCG and ProEGCG to each pharmacophore were calculated and normalised to 0–1 range. The larger the fit score, the better the protein-molecule complex. Z-score was generated by combining the fit score to the database score matrix. It considered both pharmacophore models score, as well as the statistical factors such as the significance of the binding between protein targets to the compound. It compared the fit score of a specific pharmacophore to the scores level among all scores in the pharmacophore library. The larger the positive z-score, the more significance of the target to the compound and vice versa.

Top 300 potential protein targets of EGCG and ProEGCG, obtained from PharmMapper, were uploaded separately on the DAVID Bioinformatics Resources 6.8 server (https://david.ncifcrf.gov/summary.jsp) (Dennis et al., 2003) for GO functional and KEGG pathway enrichment analysis. The identifier was set to gene symbol, and Homo sapiens was selected to limit annotations in the gene list and background list. A significant value of p < 0.05 was set as the cutoff criterion. Common and different GO terms in biological process, cellular component, molecular function and top 25 enriched KEGG pathway were identified for further analysis.

Proteins involved in the enriched KEGG pathways of ProEGCG were selected. RNA and protein expressions in endometrial stromal cells were analyzed for their ranges of RNA and protein expressions. Data were obtained from Human Protein Atlas Database (https://www.proteinatlas.org) (Uhlen et al., 2015; Thul et al., 2017). RNA and Protein expression levels were achieved by RNA-sequencing data and immunohistochemical staining patterns from internal, external sources and available protein and gene characterization public data. Gene expressions in endometriosis were obtained from Mammalian Gene expression database for uterus tissue (http://resource.ibab.ac.in/MGEx-Udb/) (Bajpai et al., 2012).

Molecular docking was employed to predict the non-covalent binding of the target proteins to EGCG and ProEGCG. Autodock Tools, version 4.2 and Autodock Vina version 1.1.2 were used (Morris et al., 2009). The pdb structure of target proteins was downloaded from protein data bank (https://www.rcsb.org). The pdb structure of EGCG and ProEGCG were converted and exported from OpenBabel. Ligand binding sites of each protein were identified from metaPocket2.0 (Huang, 2009) (https://projects.biotec.tu-dresden.de/metapocket/index.php). Each ligand was docked individually to each protein. The grid box was set at a size of 20 × 20 × 20 points (x,y,z). It covered the binding site residues and allowed the ligand to move at a flexible condition. Grid spacing was set to 1 Å. Default settings were used for other parameters. The predicted ligand binding affinities were shown as a negative Gibbs free energy (kJ/mol). The more negative the binding affinity, the stronger the bonding of complex. Only proteins that show gene expression in endometrium, protein expression in the endometrial stromal cells, are expressed in endometriosis patients, and form protein-ligand complex with the most negative Gibbs free energy were selected as final potential protein targets for further chemical binding and functional analysis.

Ligplot (https://www.ebi.ac.uk/thornton-srv/software/LIGPLOT/) (Wallace et al., 1995) was used for virtual screening after the molecular docking which gave a 2D interaction scheme. The software was used to observe the binding poses and to identify residues that involved hydrogen bonding and hydrophobic interaction in specific bonding distance as interaction radii in the complex. These were the non-covalent intermolecular interactions that influenced the binding affinity and the biological activity of drug.

String version 11 (https://string-db.org) (Szklarczyk et al., 2019) was used for protein network analysis. NMNAT1 and NMNAT3 were submitted to the multiple protein search. Medium confidence with a threshold of 0.4 was used to define the protein-protein network. We also added a total of 10 proteins to show a network around the input by default enrichment. All active interaction sources were included based on the available evidence from literature, open databases and reported experiments. Enrichment analysis with the whole genome as a background dataset was used to identify the enriched terms and pathways.

EGCG (Sigma Aldrich Cat. E4143, purity ≥ 95%) was purchased from Sigma. Prodrug of EGCG (ProEGCG) was a kind gift from Polytechnic University of Hong Kong and prepared according to literature procedure (Lam et al., 2004). ProEGCG (purity ≥ 99.8%, determined by LCMS/MS) is a peracetate derivative of EGCG, and DMSO was used as a solvent to dissolve EGCG and ProEGCG. The solution was diluted to the desired concentration with culture medium when used. Final concentration of DMSO was less than 0.1%.

This study was conducted as approved by the Hospital Authority of Hong Kong and The Chinese University of Hong Kong, under ethics number CRE-9264. Primary human endometrial stromal cells from a healthy woman (n = 1) without endometriosis (hEsc) were obtained from endometrial biopsies for cell culture experiments. The patient gave informed consent for donating eutopic endometrial tissue for endometriosis research purposes. She was 42 years old, diagnosed with subserosal fibroids but without endometriosis, as confirmed by laparoscopy. The stromal cells were isolated by using collagenase Crude: type 1 (Sigma, cat. C2674) and purity was assessed by morphological determination. The Cell culture medium was DMEM/F-12 (Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12) (Thermo Fisher, cat. 21041025), with 10% FBS (Gibco, cat. 12676-029), 1% Amphotericin-B (Sigma, cat. A2942), 1% Pen-Strep (Gibco, cat. 1514066) and 1% L-Glutamines (Gibco, cat. 25030). Cells were passaged every 2 days at 80% of the confluency. RIPA buffer (Thermo Scientific, cat. 89900) with protease inhibitor (Thermo Scientific, cat. 78439) was used to lyse the cells for in-vitro assays. Protein concentrations were determined using DC™ Protein Assay (BioRad, cat. 5000112) with BSA as standard.

hEsc was treated with EGCG and ProEGCG in a range of 0–300 μM for 48 h. Untreated cells were used as a negative control, while reagents with no cells was used as blank control. MTT (3-(4,5-Dimethylthiazol-2-yl)-2,5-Diphenyltetrazolium Bromide; ThermoFisher cat. M6494) was added to each well and incubated at 37°C for 3 h. DMSO was used as a solubilizing agent to dissolve the formazan crystals, and quantified with absorbance measurement at 570 nm. GraphPad Prism version 8.4.2. (https://www.graphpad.com/updates/). generated the MTT cell proliferation plot and calculated value of the IC50 concentration.

CESTA validated the protein-ligand binding without requiring small molecule modification, and solely depended on drug bindings, instead of drug actions (Jafari et al., 2014). Protein lysates of the cultured hEsc (2 μg/μl) were incubated with EGCG or ProEGCG at their IC50 concentrations or DMSO as vehicle control for 30 min at 25°C on a microplate incubator shaker. Complexes were heated for 3 min in seven different temperatures, ranging from 40°C to 70°C. Heat-induced unbound proteins to unfold while remaining soluble proteins were quantified. Soluble fraction of protein was isolated from the lysates by centrifugation for 15 min at 20,000 x g. Samples were then treated with SDS loading buffer for 10 min at 75°C, followed by western blotting analysis.

Denatured proteins were separated by 12.5% SDS/polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to a PVDF membrane. The membrane was blocked with 5% skim milk in TBST for 30 min at room temperature. It was then probed with the corresponding primary antibodies, 1:200 Alpha tubulin (Santa cruz, cat. Sc-8035), 1:200 Beta actin (Santa cruz, cat. Sc-47778), 1:100 NMNAT3 (Santa cruz, cat. Sc-390433) and 1:100 NMNAT1 (Santa cruz, cat. Sc-271557) with gentle agitation overnight at 4°C. α-Tubulin was used as a positive control and ß-actin was used as a negative control protein (Chakrabarty et al., 2015). On the second day, the membrane was incubated with secondary antibodies, 1:500 Goat anti Mouse IgG1 HRP (Abcam, cat. A10551) with gentle agitation for 1 h at room temperature. The blots were then developed with ECL™ Prime Western Blotting Detection Reagent (GE Healthcare Amersham™, cat. 45-002-401). Western blot analysis were performed in triplicate. Protein bands were quantified by Fiji ImageJ software version 2.1.0/1.53c (https://imagej.net/software/fiji/). GraphPad Prism version 8.4.2. generated melting curves. The data were normalized in a range between 0 and 100, and was plotted as a sigmoidal curve of temperature vs percentage of unfolded protein.

All animal works were approved by the Animal Experimentation Ethics Committee of the Chinese University of Hong Kong and performed under ethics 07/010/MIS in accordance with institutional guidelines. 8-week-old C57BL/6 female mice provided by the Laboratory Animal Service Center were housed in a pathogen-free animal house with chow and tap water. Ovariectomy was performed 7 days before transplantation. 100 μg/kg estradiol-17β (E2; Sigma, Cat. E2758) were intramuscularly injected at the time of ovariectomy and on every 5 days to synchronize the E2 level. Mice in the donor group were sacrificed to obtain the endometriotic tissues of the uterine fragments, which was then cut into endometrial fragments of 2 mm in diameter using a dermal biopsy punch (Miltex). Endometriosis model was established on day 0 by subcutaneous transplantation. A 3 mm skin incision on the abdominal wall of each recipient mice was made, thus subcutaneous pockets were created. Three endometrial tissues were placed into the pocket and then the skin incision was closed with surgical thread (Wang et al., 2013; Xu et al., 2013). Anesthesia with 100 mg/kg ketamine, 10 mg/kg xylazine in 0.9% saline water was applied on all surgery. The mice were then randomly assigned to one of the following treatments: Vehicle as the negative control, EGCG and ProEGCG, injected intraperitoneally at 50 mg/kg every day for 21°days until termination for lesions collection and after mice were sacrificed. To assess the effects of EGCG and ProEGCG (n = 10 per treatment group) in inhibiting lesions growth, length, and width of the area of the endometriotic lesions were measured using a digital caliper. Lesions area (mm²) was determined using the formula, length (mm) × width (mm) × π/4.

Endometrial implants were collected in RNAlater solution at 4°C and incubated for 24 h, followed by removing the supernatant, and the tissues were stored at −80°C for long term storage. These were used for whole genome expression microarray analysis of mouse (n = 3 per treatment group). Total RNAs were extracted by the RNeasy mini kit (Qiagen) and were subjected to gene expression profiling using the Agilent microarray platform. For all RNAs, 260/280 ratio was ∼2, while RIN was ∼9. The experiment was performed according to previous procedures (Wang et al., 2013; Xu et al., 2013). Data were acquired by Feature Extraction 10.5 and were analyzed with GeneSpring GX 11.0 software (Agilent). T-test analysis was performed and p values of <0.05 were considered as significant. Multivariate variable pattern recognition of the target and related downstream proteins of NMNAT1 and NMNAT3 was performed. We transposed the dataset to have single animals as variable and gene expression as statistical units, plotted by Devtools and Ggbiplot from Bioconductor in R studio, with the animals on loading space (Gabriel, 1971; RStudioTeam, 2020). Differentially expressed genes were clustered using Cluster 3.0 (http://bonsai.hgc.jp/%7Emdehoon/software/cluster/) (de Hoon et al., 2004), and heatmap was generated and presented by Java Treeview (http://jtreeview.sourceforge.net) (Saldanha, 2004).

To validate the microarray results and quantify gene expression levels of NMNAT1 and NMNAT3 in lesions of endometriosis mice model after EGCG or ProEGCG treatments, qRT-PCR was performed (n = 4–5 per treatment group). TB Green Premix Ex Taq (Tli RNaseH Plus) (Takara; cat RR420A) was used to amplify the gene expression and quantified by LightCycler 480 quantitative PCR system (Roche, Switzerland). Mouse GAPDH genes were used as house-keeping genes to normalize the gene expressions of studied proteins. PCR primers for NMNAT1, NMNAT3 and GAPDH were listed in Supplementary Table S1.

Lesions tissues were fixed in 10% buffered formalin for 24 h and were embedded in paraffin blocks, cut into slices of 4 μm thickness. Sections were deparaffinized, incubated with 3% H₂O₂ in the dark to block the endogenous peroxidase activity. The sections were then incubated in boiled 10 mM sodium citrate buffer, pH 6 for antigen retrieval, and was boiled in a conventional microwave for 2 min for two times. It was then cooled and followed by blocking with 5% goat serum in 1% BSA/PBS and goat F(ab) anti-mouse IgG H&L (ab6668, 1 mg/ml, Abcam) to reduce non-specific signal and non-specific binding of the mouse serum to mouse tissue. The tissues were then incubated with primary antibodies 1:100 dilution of NMNAT3 (Santa cruz, cat. sc390433) and 1:100 dilution of NMNAT1 (Santa cruz, cat. sc271557) overnight, followed by secondary antibodies Goat Anti-Mouse IgG H&L (Abcam, cat. ab97023) at a dilution of 1:200 for 1 hour and developed in 3,3-diaminobenzidine (DAB) for 3–5 min, all in the appropriate dilution. After mounting, Leica DM6000B microscope was used for image capturing. For semi-quantitative analysis, NMNATs protein expression were quantified by H-score. The scores were generated blindly and independently by two authors (S. W. H and R. Z. Z). The H-score was calculated by summing up the percentage of different intensity level-stained cells: strongly stained cells (3+), moderately stained cells (2+), weakly stained cells (1+) and unstained cells (0) of a fixed field, on a continuous scale of 0–300, H-score = [1 × (% 1+) + 2 × (% 2+) + 3 × (% 3+)] (Fedchenko and Reifenrath, 2014; Parris et al., 2014). Each field was selected randomly and was repeated to capture five fields. The final cell count was reported as average of the five fields, and H-score was reported as whisker plot with median shown.

In the animal experiments, between groups comparisons were made using one-way ANOVA analysis. Tukey post-hoc test was used for multiple comparisons. Statistical analyses were performed by Graphpad Prism. p-value of ≤0.05 was considered statistically significant.