Punicalagin Targets Atherosclerosis: Gene Expression Profiling of THP-1 Macrophages Treated with Punicalagin and Molecular Docking

Atherosclerosis is an important cause of cardiovascular disorders worldwide. Natural botanical drugs have attracted attention due to their antioxidant, anti-inflammatory, and antiatherogenic properties in the treatment of atherosclerosis. Punicalagin is the major bioactive component of pomegranate peel, and has been shown to have antioxidant, anti-inflammatory, antiviral, anti proliferation, and anticancer properties. To explore its antiatherogenic properties at a molecular level, we investigated the genome-wide expression changes that occur in differentiated THP1 cells following treatment with a non-toxic dose of punicalagin. We also conducted a molecular docking simulation study to identify the molecular targets of punicalagin.


Introduction
Cardiovascular disease (CVD) is the leading cause of death worldwide [1,2]. One of the most important causes of CVD is atherosclerosis. Atherosclerosis is a progressive, chronic inflammatory condition characterized by the initiation and perpetuation of atherosclerotic lesions, which may erode or rupture leading to clinical events such as angina, myocardial infarction, or cerebrovascular attack [3]. It is initiated through the activation of the arterial endothelium by several risk factors leading to the recruitment of immune cells-particularly T-lymphocytes and monocytes. The latter differentiate into macrophages-a process that is accompanied by increased expression of scavenger receptors-and then transform into lipid-loaded foam cells [4]. Immune competent cells producing proinflammatory cytokines are abundant in atherosclerotic lesions and are proposed to play an important role in the progression of the disease [4,5].
Polyphenols are gaining increasing acceptance as therapeutic agents for use in diverse diseases, including CVD [6,7]. Polyphenols naturally exist in plants and plant products, including fruits, vegetables, nuts, herbs, cocoa, and tea. Historically, biologic actions of polyphenols have been attributed to antioxidant activities, but recent evidence suggests sub-cultured when they reached approximately 80% confluence (8 × 10 5 cells/mL). To promote differentiation to macrophage, confluent THP-1 cells (5-8 × 10 5 cells/mL) were treated with 160 nM PMA for 24 h [22]. An in vivo mice study had used 140 µg/100 µL (1.4 µg/µL) of punicalagin as subcutaneous dose [12]. We chose the 10 µM concentration of punicalagin based on earlier published [23] in vitro work and from our preliminary dose range determination studies, as well as by extrapolations of the in vivo study. Therefore, THP-1 cells were treated with 10 µM punicalagin for 24 h and control cells were treated with vehicle (Dimethyl sulfoxide), and their levels were maintained at non-toxic levels (0.025-0.05%) to the cells.

Total Cellular RNA Preparation from the Control and Treated THP1 Cells
Total RNA was isolated from~2 × 10 6 treated (punicalagin) and control (DMSO) THP-1 cells after harvesting them by scraping. Total RNA was extracted using a Qiagen RNeasy Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. RNA concentrations were determined using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA).

Preparation of the Sense Strand DNA for Microarray Gene Expression Analysis
Transcriptional expression profiling was performed with Affymetrix GeneChips (Gene 1.0ST, Affymetrix, Santa Clara, CA, USA) according to the conventional Affymetrix eukaryotic RNA labelling protocol (Affymetrix). Briefly, 250 ng of total RNA isolated from the control (DMSO) and treated (punicalagin) THP1 cells was first converted to single-stranded sense strand DNA (cDNA) in two cycles using the whole transcript (WT) cDNA synthesis amplification kit and sample clean-up module. In the first cycle, the total RNA was converted to double-stranded cDNA using random hexamers tagged with a T7 promoter sequence. Each strand of the double-stranded cDNA was then used as a template to synthesize antisense RNA (cRNA). In the second cycle, the cRNA was reverse-transcribed into sense strand DNA in the presence of random hexamers (3 g/mL) and dNTPs mix (10 mM).

Sense Strand DNA Labelling and Hybridization Gene 1.0 S.T Arrays
The sense strand DNA was cleaned up using the sample clean-up kit and then cleaved into small fragments using a mixture of UDP and apurinic/apyrimidinic endonuclease 1. The fragmented sense strand DNA was then end-labelled through a terminal transferase reaction that incorporates biotinylated di-deoxynucleotides using the WT terminal labelling kit. The fragmented biotinylated sense strand DNA was then hybridized to the Affymetrix Human Gene 1.0S.T array at 45 • C for 16 h in a hybridization Oven 640. After hybridization, the arrays were stained and then washed in the Affymetrix Fluidics Station 450 under standard conditions. The stained arrays were then scanned at 532 nm using an Affymetrix GeneChip Scanner 3000, and CEL files for each array were generated using the Affymetrix Gene-Chip ® Operating Software (GCOS). GCOS defines the probe cells and computes an intensity for each cell; complete probe array images were saved with a data image file extension (.dat, cel).

Microarray Data Normalization and Analysis
Affymetrix CEL files were used for raw data for image analysis and probe quantification using Partek Genome Suit 7.0. Normalization was done with the Robust Multi-chip Average (RMA) program that processes a group of CEL files simultaneously. The default options of background correction, quantile normalization, and log transformation were used. Calculated raw probe intensity data was used to derive fold change and p-value. Principal component analysis (PCA) was performed on all probes to visualize high-dimensional data and assess quality control, as well as overall variance in gene expression between the two treatments. Analysis of variance (ANOVA) was applied on the complete data set and the list of differentially expressed genes (DEGs) was then generated using a false-discovery rate (FDR) of 0.05 with 2-fold change cut-off. Unsupervised two-dimensional average linkage hierarchical clustering was performed using Spearman's correlation as a similarity matrix.

Molecular Pathway Analysis
Ingenuity Pathways Analysis (IPA) software version 338830M (Ingenuity Systems, Redwood City, CA, USA) was used to find significant canonical pathways, biological networks, biological functions, and phenotypes/disease associated with present study. DEGs along with their corresponding Affymetrix probe set ID/gene symbol/Entrez gene ID as clone identifiers, p-values, and fold change values were uploaded into IPA for functional analysis. The percentage and number of uploaded genes/molecules matching to genes of a canonical pathway were measured as Z-score, ratio, or Fisher's exact test for significance. The Molecule Activity Predictor was employed to predict the activation or suppression effects of a gene/molecule on other molecules of pathway.

Real-Time PCR Validation
Validation of the microarray data using quantitative real-time PCR (qRT-PCR) was carried out in StepOnePlusTM RT-PCR System (ThermoFisher, Waltham, MA, USA). Total RNA isolated and quality checked for the microarray experiments was used for cDNA synthesis for the qRT-PCR experiments. Manufacturer's protocol was followed for cDNA synthesis using ImProm-II™ Reverse Transcription System kit (Promega, Madison, WI, USA). The primer pairs listed in Table 1 were used for amplification of the appropriate target and endogenous control genes. All experiments were done in triplicates biological samples, data analyses were performed by Microsoft Excel 2010, and data were plotted by GraphPad Prism 8 (GraphPad Software, Inc., San Diego, CA, USA).

Molecular Docking
The X-ray crystal structure of all receptors selected for the study were downloaded from Protein Data Bank. Homology modeled structure of MSR1 macrophage scavenger receptor 1-3D structure was modeled based on homology using SWISS-MODEL. The three-dimension structures of CD36, TLR4, MSR1, LRP1, NR1H3, PPARγ, and TRAF1 were defined as target macromolecules in Autodock Vina Wizard of (PyRx, RRID: SCR_018548) software. The three-dimensional structure of compounds (palmitic acid, GW3965, pinocembrin, lobeglitazone and punicalagin) were optimized in Discovery Studio Visualizer software (version 2.5) and defined as the ligands for docking using PyRx. The active site analysis was conducted by calculating (x, y, and z coordinates) of palmitic acid, GW3965, pinocembrin, and lobeglitazone in the PyRx software for CD36, NR1H3, TLR4, and PPARγ, respectively. Entire receptors surfaces were used for identifying the active site in the case of MSR1, LRP1, and TRAF1. Blind docking was performed for the receptors MSR1, LRP1, and TRAF1 as no ligands were known. Coordinates of (centre_x, y, z) and dimensions in A (x, y, z) were maximized to cover the entire receptor structure to identify the active site/docking sites. Energy minimization was performed using the "uff" forcefield and conjugate gradients optimization algorithm for standard and test compounds before docking.

Normalization and Visualization
Microarray data were initially visualized using a principal components analysis scatter (PCA) plot ( Figure 1). The data from the control group (n = 3) are shown in red, whereas those in the treated group (n = 3) are indicated in blue. Each ball represents gene expression data generated from a sample that was applied on a gene chip. The results of PCA of transcriptomic data as per their overall expression pattern showed that the samples from the same type of treatment cluster together. From the plot, sample outliers were not detected both in the control and in treated group. The control group is also distinguishably separated from those of the treated samples. Figure 2 shows hierarchical clustering of differentially expressed genes in THP1 macrophages in response to treatment with punicalagin. Each lane represents an array. Upregulated genes in treated samples are shown in red, and down-regulated genes are shown in purple. The pattern of gene expression was similar in the biological replicates. Comparison of the genome-wide expression of treated and control groups by ANOVA revealed 373 differentially expressed genes, including 347 up-regulated and 26 down-regulated with a cut-off of p value < 0.05 and fold change more than ±2 (Supplementary Table S1). The top 20 up-and down-regulated genes are represented in Tables 2 and 3, respectively. FABP4, CD36, MSR1, LPL, and TDO2 were upregulated more than ten-fold in the treated population. Genes that were significantly downregulated include metallothionein genes MT2A and MT1P1.

Identification of Canonical Pathways and Associated Disease Functions
Ingenuity pathway analysis (IPA) on the 373 upregulated genes was carried o elucidate the affected biological processes. This analysis revealed canonical path such as Cholesterol Biosynthesis I, II, III, super pathway of cholesterol biosynthesis LXR/RXR activation, (Table 4, Figure 3) (Supplementary Table S2). Network intera between the 15 most differently expressed genes in (LXR/RXR) the activation pathw represented in (Figure 4), showing 20 related genes and 301 total links related t LXR/RXR pathway. IPA also revealed cell movement, migration of cells, leukocyte m tion, non-hematological solid tumor, non-hematologic malignant neoplasm, develop

Identification of Canonical Pathways and Associated Disease Functions
Ingenuity pathway analysis (IPA) on the 373 upregulated genes was carried out to elucidate the affected biological processes. This analysis revealed canonical pathways such as Cholesterol Biosynthesis I, II, III, super pathway of cholesterol biosynthesis, and LXR/RXR activation, (Table 4, Figure 3) (Supplementary Table S2). Network interaction between the 15 most differently expressed genes in (LXR/RXR) the activation pathway is represented in (Figure 4), showing 20 related genes and 301 total links related to the LXR/RXR pathway. IPA also revealed cell movement, migration of cells, leukocyte migration, non-hematological solid tumor, non-hematologic malignant neoplasm, development

Identification of Canonical Pathways and Associated Disease Functions
Ingenuity pathway analysis (IPA) on the 373 upregulated genes was carried out to elucidate the affected biological processes. This analysis revealed canonical pathways such as Cholesterol Biosynthesis I, II, III, super pathway of cholesterol biosynthesis, and LXR/RXR activation, (Table 4, Figure 3) (Supplementary Table S2). Network interaction between the 15 most differently expressed genes in (LXR/RXR) the activation pathway is represented in (Figure 4), showing 20 related genes and 301 total links related to the LXR/RXR pathway. IPA also revealed cell movement, migration of cells, leukocyte migration, non-hematological solid tumor, non-hematologic malignant neoplasm, development of vasculature, cancer, solid tumor, malignant solid tumor, synthesis of lipid, etc., as most associated diseases or functions (Supplementary Tables S3 and S4).  Tables S3 and S4).

Validation of Microarray Data
The expression of selected up-regulated genes NR1H3 (LXR), MSR1 (SRA), and CD36 relative to that of GADPH was also quantitated using real-time RT-PCR to validate findings from the microarray data. As shown in Figure 5, a robust increase in mRNA derived from these three genes was observed in the punicalagin treated population vs. control. CD36 showed approximately a 13-fold increase in transcript quantities, whereas MSR1 and NR1H3 mRNA increased by 140 folds and 5 folds, respectively, in the treated THP1 macrophage population as compared to the control cells, thus corroborating the increase in mRNA levels of these genes observed following microarray analysis.

Validation of Microarray Data
The expression of selected up-regulated genes NR1H3 (LXR), MSR1 (SRA), and CD36 relative to that of GADPH was also quantitated using real-time RT-PCR to validate findings from the microarray data. As shown in Figure 5, a robust increase in mRNA derived from these three genes was observed in the punicalagin treated population vs. control. CD36 showed approximately a 13-fold increase in transcript quantities, whereas MSR1 and NR1H3 mRNA increased by 140 folds and 5 folds, respectively, in the treated THP1 macrophage population as compared to the control cells, thus corroborating the increase in mRNA levels of these genes observed following microarray analysis. Figure 5. Validation of the microarray data using quantitative real-time PCR. The bar chart shows the gene expression patterns (presented as fold change) of selected significantly up-regulated genes NR1H3 (LXR), MSR1 (SR-A), and CD36 calculated using quantitative real-time PCR. All PCR data were normalized to the intensity of GAPDH as a housekeeping gene. Data represented as mean + SEM following three independent experiments. Statistical analysis was performed using a one-way ANOVA with a Sidak's test analysis where **: p < 0.01, ****: p < 0.0001.

Molecular Docking Analysis
To identify the molecular target of punicalagin, we conducted molecular docking simulation study. According to the microarray analysis, we focused on seven key receptors: CD36 (Cluster of Differentiation (Table 5). We found a significant binding affinity of punicalagin with CD36 (−9.3 Kcal/mol) by establishing hydrogen bond interactions. Binding of CD36 was computed with its well-described ligand palmitic acid, and observed lower binding affinity compared to punicalagin (−6.8 Kcal/mol). Interestingly, it was found that these two CD36-binding molecules exhibit distinct interaction profiles with different binding sites ( Figure 6A,B).
In the case of TLR4, pinocembrin-an antagonist for TLR4-showed binding affinity as (−8.2 Kcal/mole), whereas punicalagin showed a slightly higher value as (−9.0 Kcal/mole). Interaction profiling of pinocembrin and punicalagin showed a maximum of seven common hydrogen bonding interactions ( Figure 6C,D). This observation suggests that punicalagin is a high affinity inhibitor of TLR4 receptors. As MSR1 structure is not experimentally solved by X-ray crystallography and nuclear magnetic resonance technique, the computational technique of homology modelling was used for 3D structure determination of MSR1. Punicalagin binding affinity was observed to be (−7.4 Kcal/mol). It was observed that punicalagin directly interacted with hydrogen bonding to amino acid residues (Figure 6E).
The binding affinity of punicalagin to LRP1 was estimated to be (−6.4 Kcal/mole), involving hydrogen bonding interactions with amino acid residues in LRP1 ( Figure 6F). The binding affinity of punicalagin to the NR1H3 receptor was (−7.1 Kcalcal/mole). In comparison, the binding affinity of the NR1H3 agonist GW3965 was (−13.9 kcal/mole). Interestingly, there was not a single common interaction observed between these two molecules ( Figure 6G,H), suggesting that punicalagin could interact with NR1H3 in a different way than other NR1H3 agonists.
In the case of PPARγ, lobeglitazone-an agonist for PPARγ-demonstrated binding affinity of (−8.6 Kcal/mole), whereas punicalagin had a slightly higher value of (−8.7 Kcal/mole). Interaction profiling indicated a maximum of three common hydrogen bonding interactions between lobeglitazone and punicalagin ( Figure 6I,J). This indicated that punicalagin might be acting as an agonist in the same way as a lobeglitazone against Figure 5. Validation of the microarray data using quantitative real-time PCR. The bar chart shows the gene expression patterns (presented as fold change) of selected significantly up-regulated genes NR1H3 (LXR), MSR1 (SR-A), and CD36 calculated using quantitative real-time PCR. All PCR data were normalized to the intensity of GAPDH as a housekeeping gene. Data represented as mean + SEM following three independent experiments. Statistical analysis was performed using a one-way ANOVA with a Sidak's test analysis where **: p < 0.01, ****: p < 0.0001.

Molecular Docking Analysis
To identify the molecular target of punicalagin, we conducted molecular docking simulation study. According to the microarray analysis, we focused on seven key receptors:  (Table 5). We found a significant binding affinity of punicalagin with CD36 (−9.3 Kcal/mol) by establishing hydrogen bond interactions. Binding of CD36 was computed with its well-described ligand palmitic acid, and observed lower binding affinity compared to punicalagin (−6.8 Kcal/mol). Interestingly, it was found that these two CD36-binding molecules exhibit distinct interaction profiles with different binding sites ( Figure 6A,B).
In the case of TLR4, pinocembrin-an antagonist for TLR4-showed binding affinity as (−8.2 Kcal/mole), whereas punicalagin showed a slightly higher value as (−9.0 Kcal/mole). Interaction profiling of pinocembrin and punicalagin showed a maximum of seven common hydrogen bonding interactions ( Figure 6C,D). This observation suggests that punicalagin is a high affinity inhibitor of TLR4 receptors. As MSR1 structure is not experimentally solved by X-ray crystallography and nuclear magnetic resonance technique, the computational technique of homology modelling was used for 3D structure determination of MSR1. Punicalagin binding affinity was observed to be (−7.4 Kcal/mol). It was observed that punicalagin directly interacted with hydrogen bonding to amino acid residues ( Figure 6E).
Curr. Issues Mol. Biol. 2022, 2, FOR PEER REVIEW 10 PPARγ. On the other hand, punicalagin binding affinity against TRAF1 was observed to be (−9.0 Kcal/mol). Its interaction profiling showed hydrogen bonding interaction as in ( Figure 6K).  The binding affinity of punicalagin to LRP1 was estimated to be (−6.4 Kcal/mole), involving hydrogen bonding interactions with amino acid residues in LRP1 ( Figure 6F). The binding affinity of punicalagin to the NR1H3 receptor was (−7.1 Kcalcal/mole). In comparison, the binding affinity of the NR1H3 agonist GW3965 was (−13.9 kcal/mole). Interestingly, there was not a single common interaction observed between these two molecules ( Figure 6G,H), suggesting that punicalagin could interact with NR1H3 in a different way than other NR1H3 agonists.
In the case of PPARγ, lobeglitazone-an agonist for PPARγ-demonstrated binding affinity of (−8.6 Kcal/mole), whereas punicalagin had a slightly higher value of (−8.7 Kcal/mole). Interaction profiling indicated a maximum of three common hydrogen bonding interactions between lobeglitazone and punicalagin ( Figure 6I,J). This indicated that punicalagin might be acting as an agonist in the same way as a lobeglitazone against PPARγ. On the other hand, punicalagin binding affinity against TRAF1 was observed to be (−9.0 Kcal/mol). Its interaction profiling showed hydrogen bonding interaction as in ( Figure 6K).

Discussion
Macrophages are central to the initiation and progression of atherosclerosis and can be highly appropriate targets for antiatherogenic therapy. Punicalagin was demonstrated to exert promising anti-atherosclerotic effects on macrophages, but molecular targets remain to be identified. In this study, we performed the first large-scale gene expression analysis to gain insight into the mechanism of action of punicalagin. Key regulators of identified pathways were then subjected to in silico docking analysis to predict the molecular target of punicalagin. Our results showed that punicalagin anti-atherosclerotic properties rely on multiple targets in THP-1 macrophages. Indeed, we predicted that punicalagin regulates the activity of CD36, TLR4, MSR1, LRP1, NR1H3, PPARγ, and TRAF1. These predictions were strongly supported by transcriptomic analyses showing that punicalagin upregulates the expression of genes involved in cholesterol efflux and lipid metabolism, and reduces the expression of genes involved in inflammation, proliferation, cell migration, and adhesion.
The results of the transcriptomic analyses showed that among the top 20 significantly upregulated genes in response to punicalagin, several were associated with cholesterol or lipids metabolism (FABP4, CD36, LPL, FADS2). This was further supported by ingenuity pathway analysis that revealed that the four most modified canonical pathways in response to punicalagin treatment were associated to cholesterol synthesis (cholesterol biosynthesis pathway I, II, and III and super pathway of cholesterol biosynthesis) and that lipid synthesis is one of the top modified functions. These results are in line with previous studies that revealed a reduction of circulating cholesterol levels following punicalagin treatment, resulting from a reduction of cholesterol synthesis and accumulation in macrophage, associated with an increase of LDL influx to macrophages [12,15,23]. Molecular docking analysis revealed four high-probability targets of punicalagin involved in cholesterol and lipid metabolism: CD36, LRP1, NR1H3, PPARγ. The interaction of punicalagin with CD36 is of particular interest, as the binding is predicted to be higher than with palmitate-the traditional ligand of CD36 [27]. Future studies are warranted to support this finding.
Transcriptomic analyses also support previous observations showing a reduction of proliferation, migration, and adhesion in cancer cells [22,28,29] and macrophages [24] in response to punicalagin treatment. Indeed, among the top 20 significantly upregulated genes in response to punicalagin, several were associated with cell proliferation and migration (CEMIP, GPC4, FCER2, TDO2, MTUS1, PDGFRL, MMP1) or adhesion (MSR1). This was further supported by ingenuity pathway analysis that revealed that some modified canonical pathway in response to punicalagin treatment were associated to proliferation, migration, and adhesion (Granulocyte Adhesion and Diapedesis, Colorectal Cancer Metastasis signalling, agranulocyte Adhesion and Diapedesis Leukocyte Extravasation Signalling). Importantly, cell movement and migration of cells were among the top modified functions. Molecular docking analysis predicted a high likelihood for punicalagin binding to MSR1. This could represent a key target for punicalagin, as this receptor has been showed to regulate the PI3K/AKT/GSK3β/β-catenin pathway, which is involved in lipid metabolism, proliferation, adhesion, and inflammation [22,[30][31][32].
Finally, our transcriptomic analyses also confirm the anti-inflammatory role of punicalagin [18,20,33], as evidenced by the fact that among the top 20 significantly downregulated genes, 10 were associated with inflammation (TRAF1, FKBP5, MT2A, MT1P1, PLAGL1, NR4A2, BCL3, MT1X, PIM2, and HSPA4L). Molecular docking analyses predicted that TLR4 and TRAF1, two important regulators of inflammation [34,35], could be targeted by punicalagin. Further studies are needed to experimentally confirm these results. Interestingly, while the antiatherogenic effect of punicalagin was suggested to rely on its antioxidant capacity [14], transcriptomic analyses did not reveal important modifications of genes involved in the antioxidant response. This result supports previous molecular docking analyses that revealed no significant binding of punicalagin to key antioxidant enzymes [36], and reinforces the importance of punicalagin's effect on lipid metabolism and cell inflammation, proliferation, migration, and adhesion.

Conclusions
Overall, our data identified key pathways associated with the anti-atherosclerotic properties of punicalagin in THP-1 macrophages. Punicalagin's mechanism of action relies on the regulation of multiple pathways including lipid metabolism, cell inflammation, proliferation, and migration. Whilst future studies are required to support our findings, we identified with high likelihood key regulators of these process that could be targeted by punicalagin.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cimb44050145/s1, Table S1: Differentially expressed genes (DEGs), Included 347 up-regulated and 26 down-regulated with a cut-off of p value < 0.05 and fold change more than ±2. Table S2: Canonical pathways generated by ingenuity pathway analysis (IPA). Table S3: Disease functions generated by ingenuity pathway analysis (IPA). Table S4: Top10 predicted z scores of diseases or functions (inhibition or activation).