SOCS1 regulates senescence and ferroptosis by modulating the expression of p53 target genes

SOCS1 inhibition affects a selective group of p53 target genes

To investigate how SOCS1 modulates the p53 pathway we used a model of oncogene-induced senescence where p53 activation is dependent on SOCS1 [14,15]. In this model, a constitutively active allele of STAT5A (cS5A) is introduced into human fibroblasts IMR90-E7. These cells express the papillomavirus protein E7 to inactivate the retinoblastoma tumor suppressor, allowing us to focus on the p53 pathway contribution to cellular senescence. Expression of cS5A in these cells induced a cell cycle arrest and reduced the expression of proliferation markers (MCM6 and phospho-histone H3) that was rescued by expression of a small hairpin RNA against SOCS1 (Figure 1A-C). Consistent with the reported role of SOCS1 in p53 phosphorylation at serine 15, cS5A failed to induce this modification in cells expressing the shRNA against SOCS1 (Figure 1D). We then compared the transcriptome of IMR90-E7 cells expressing cS5A and co-expressing an shRNA against SOCS1 or a non-targeted shRNA using RNA purified from cells 7 days after introduction of cS5A. The microarray data was deposited in NCBI’s Gene Expression Omnibus (GEO) GSE98216. Genes whose expression changed +/- 1.5 fold were used for pathway analysis with the online bioinformatics platform DAVID. The most significant pathways regulated by SOCS1 in cells expressing cS5A are indicated in Figure 1E and S1. They include genes coding for secreted proteins (Figure S1A-C) and a few genes in the p53-signaling pathway (Figure 1E and Table 1). The bypass of senescence induced by cS5A and the inhibition of p53 activity was also confirmed by a second shRNA against SOCS1 (shS1b), (Figure S2). Therefore, although these results confirm the requirement for SOCS1 in p53 phosphorylation and the senescence response to aberrant STAT5A activation, the pattern of p53 target genes identified does not include the classic p53 target genes associated to the DNA damage response.

To further investigate whether SOCS1 modulates the p53 pathway, we performed Gene set enrichment analysis (GSEA) of the microarray data. We found that SOCS1-disabled cells have high levels of genes in the set TANG_SENESCENCE_TP53_TARGETS_DN that contain, genes upregulated by a dominant negative p53 in normal human fibroblasts (Figure 1F). This gene set mostly includes cell cycle regulated genes, which are known to be regulated by RB. However, the RB pathway was disabled in the IMR90-E7 cell line used in our experiments. These results thus imply a SOCS1- p53-dependent pathway that regulates cell cycle genes independently of the RB tumor suppressor. SOCS1-disabled cells also have a downregulation of the gene set MODULE_474 (http://robotics.stanford.edu/~erans/cancer/) containing several genes upregulated in senescent cells such as IGFBP2, 3, 5, 6 and 7, NOV and SOCS3 (Figure 1G). Further proof that SOCS1 controls the expression of p53 target genes was obtained using the platform DiRE that analyses the promoters of gene sets for signatures of transcription factors. P53 was one of the top transcription factors associated to transcriptome changes induced by an shRNA against SOCS1 (Figure 1H). In addition, this algorithm identified 12 new candidates for p53-target genes (Table 2). Many of these genes are poorly expressed in human cancers suggesting that they could be novel tumor suppressors (Table 3). The microarray data was validated using qPCR for either known p53-target genes or the new candidates suggested by the DiRE platform (Figure 1I). We also used data in TCGA from hepatocellular carcinoma, a tumor type where promoter DNA methylation often silences SOCS1 expression [16]. We found a significant correlation between the expression of SOCS1 and the SOCS1-dependent p53 target genes defined in human fibroblasts (Figure 2). For several genes (SLC7A11, SAT1, SERPINE1, IGFBP7, GADD45 and ATF3) these correlations decreased in samples from p53 mutant tumors (Figure S3). Taken together, our results show that SOCS1 controls a unique set of p53 target genes.

To investigate the biological importance of SOCS1-dependent p53 target genes we analysed their overlap with gene sets mediating specific p53-dependent responses. First, genes regulated by SOCS1 matched gene sets that regulate the response to chemotherapy (Figure S4A), which is largely influenced by the p53 pathway [17]. Second, SOCS1-disabled cells expressed high levels of genes upregulated by a dominant negative BRCA1 allele (Figure S4B), which also disrupts the p53 pathway [18]. Third, SOCS1-disabled cells have a decrease in the expression of genes that blocked angiogenesis in endothelial cells (Figure S4C). These genes include IGFBP3 [19] and COL4A2 which encodes the potent antiangiogenic factor canstatin [20]. It is well known that p53 inhibits angiogenesis [21,22], and our gene expression data suggest that SOCS1 modulates this p53 function as well. Finally, SOCS1-regulated genes also overlapped with a set of genes induced by oxidized phospholipids (Figure S4D), which has been recently linked to iron-dependent cell death or ferroptosis [23]. P53 sensitizes cells to ferroptosis by repressing the cystine transporter SLC7A11 [4] and inducing the polyamine metabolic enzyme SAT1 [24]. The regulation of those two genes in cS5A expressing cells was SOCS1-dependent (Table 1) suggesting a role for SOCS1 in ferroptosis.

The inhibition of SOCS1 expression in IMR90-E7 bypasses cS5A-induced senescence, raising the possibility that the defects in p53 target gene expression we described above are the consequence of senescence inhibition and are not directly linked to SOCS1. In IMR90 cells, where the retinoblastoma pathway remains intact, inactivation of SOCS1 does not bypass cS5A-induced senescence [14]. This is due to the known fact that the RB pathway is sufficient to regulate senescence in the absence of p53 [25,26]. We thus took advantage of this fact to investigate whether expression of p53 target genes still required SOCS1 in these cells. For most of the genes measured, disabling SOCS1 also inhibited the expression of p53 targets as well (Figure 3A-B) although cells remained senescent (Figure 3C-E). We conclude that the defects in p53 target gene expression we have seen after disabling SOCS1 are not the result of cell cycle re-entry after bypass of senescence and suggest a direct effect of SOCS1 on p53.

SOCS1 is sufficient to activate p53 and regulate the expression of its target genes

Previous research established a function for SOCS1 as an adaptor protein facilitating the phosphorylation of p53 at serine 15 by the DNA damage response activated kinases ATM and ATR [14]. In culture, cells are exposed to high concentrations of oxygen and growth factors, increasing the probability of DNA damage by either replication stress or reactive oxygen species. We reasoned that enforcing SOCS1 expression might cooperate with these factors and engage the p53 pathway. Indeed, expressing SOCS1, in U2OS cells (Figure 4A-B) or IMR90 cells (Figure 4C-E) increased the expression of most of the genes that were positively regulated by p53 in cS5A-induced senescence and decreased the expression of SLC7A11, a gene that is repressed by p53. As reported before, SOCS1 was sufficient to trigger senescence in close to 50% of the cells overexpressing the protein (Figure 4F). The increase in expression of p53 target genes by SOCS1 was, for the most part, dependent on p53 since knockdown of p53 (Figure 4G-H) abolished their stimulation and prevented senescence (Figure 4I). However, the decrease in SLC7A11 was not blocked after knockdown of p53 suggesting that SOCS1 regulates the expression of this gene by additional mechanisms. Of note, not all p53 target genes that required SOCS1 in cS5A-induced senescence were induced by SOCS1 expression alone. The following genes (IGFBP3, GFBP5, DUSP1, ALDH1L12, IGFBP7, DDIT3, PPP1R3C, NOLC1, NTN4 and SVEP1) were not significantly induced by SOCS1 in either U2OS or IMR90 cells. In addition, RGCC, ESM1, ATF3 and PROCR were induced in U2OS but not in IMR90 cells. These results indicate that SOCS1 partially reproduces oncogenic STAT5A signaling to p53 but additional STAT5A functions also contribute to p53 activation.

SOCS1 sensitizes cells to ferroptosis

The mechanisms by which p53 prevents tumor progression include apoptosis, senescence and ferroptosis [27,28]. The latter is an iron-dependent cell death mechanism that involves reactive oxygen species and lipid oxidation and that does not depend on caspases or Bcl2 family members Bak and Bax [28]. The SOCS1-dependent p53 target genes SLC7A11 [4] and SAT1 [24] play key roles in ferroptosis suggesting that SOCS1 could be a novel regulator of this process. As anticipated, expression of SOCS1 in U2OS or IMR90 cells sensitized them to the ferroptosis inducer tert-butyl-hydroperoxide (TBH). This effect of SOCS1 correlated with its ability to reduce the expression of the cystine transporter SLC7A11 (Figure 5A-F) and was efficiently blocked by treatment with deferoxamine, an iron chelator known to inhibit ferroptosis (Figure 5G-H and Figure S4). Ferroptosis depends on the actions of oxidized lipids, which are detoxified by the glutathione (GSH)-dependent enzyme GPX4 [29]. Cystine is used in the cells to generate cysteine, a metabolite required for the synthesis of glutathione (GSH). As expected, expression of SOCS1 reduced the levels of GSH (Figure 5I) explaining in part its ability to sensitize cells to ferroptosis.

The SOCS1-p53 axis is active in oncogene-induced senescence (OIS) but additional factors control ferroptosis in these cells

OIS is triggered by activated oncogenes and involves replication stress, mitochondrial dysfunction, reactive oxygen species and the DNA damage response [30–32]. We anticipated that SOCS1 is generally upregulated during OIS and not only in the specific case of cS5A-induced senescence. This is because senescent cells secrete a variety of inflammatory mediators that activate JAK-STAT signaling [33], the major regulator of SOCS1 expression [14,34]. We thus induced senescence in IMR90 normal human fibroblasts with constitutively active RAS (RASV12), cS5A or the tumor suppressor PML, that acts downstream of RAS and cS5A to induce senescence (Figure 6A). As published before for cells expressing cS5A [14], we found SOCS1 localizing with phospho-ATM [35] in DNA damage foci in senescent cells, which are known to contain phosphorylated p53 at serine 15 [36] (Figure 6B). Also, induction of senescence by RASV12 or PML led to induction of SOCS1-dependent p53 target genes and the downregulation of the p53-repressible gene SLC7A11 (Figure 6C). Hence, despite the differences in the mechanisms triggering senescence by RasV12, STAT5A, SOCS1 or PML, they regulate a common set of p53 target genes.

Next, we performed ferroptosis assays in cells expressing RASV12, STAT5A or PMLIV using the ferroptosis inducer TBH and the ferroptosis inhibitor deferoxamine. For this assay we plated the same amount of cells for each condition and found that neither RasV12 nor STAT5A sensitized cells to ferroptosis, despite both having an increase in endogenous SOCS1 expression and a reduction in SLC7A11. Intriguingly, PML expression turned cells highly resistant to ferroptosis providing a mechanistic insight into a ferroptosis resistance pathway in senescent cells (Figure 6D-E). Glutathione levels correlated with ferroptosis resistance in senescent cells (Figure 6F) suggesting compensatory mechanisms of GSH biosynthesis as a possible mechanism for resistance to ferroptosis. Targeting anti-ferroptosis pathways activated in senescent cells might selectively kill them in the same way as targeting anti-apoptotic proteins.

SOCS1 activates p53 via both phosphorylation and stabilization

The pattern of p53 target genes regulated by SOCS1, including ferroptosis regulators, is unique and anticipates unique mechanisms of p53 activation by SOCS1. So far, SOCS1 has been linked to the regulation of serine 15 phosphorylation of p53 [14,15], a modification that is also efficiently induced during an acute DNA damage response, which does not necessarily evolve into a permanent senescent cell cycle arrest [37,38]. Therefore, the ability of SOCS1 to induce senescence must involve additional mechanisms capable of changing the quality of the p53 response.

To investigate how SOCS1 modulates p53 activity, we first induced endogenous p53 using the DNA damaging drug doxorubicin and measured both p53 phosphorylation and total p53 levels. SOCS1 expression increased p53 levels during the course of a doxorubicin treatment without changing the levels of serine 15-phosphorylation of p53, which are highly induced by the drug (Figure 7A and B). In other words, in the context of oncogenic stress that involves relatively modest rates of DNA damage, SOCS1 stimulates serine 15 phosphorylation of p53 by promoting p53-ATM interactions as described before [14]. In cells treated with doxorubicin, high levels of DNA damage maximally stimulate serine 15 phosphorylation of p53 in a SOCS1-independent manner but SOCS1 can still stabilize p53 in this condition.

In order to find additional mechanisms of p53 activation by SOCS1, we looked at the SOCS1 interactome [15]. Interestingly, KAP1, a repressor of p53 [39], immunoprecipitated with SOCS1 [15]. To confirm this finding we expressed wild type SOCS1 in U2OS cells and immunoprecipitated endogenous KAP1 using a specific antibody. We found SOCS1 in KAP1 immunoprecipitates but not in immunoprecipitates obtained with a control antibody (Figure 7C). The reciprocal co-IP analysis confirmed that endogenous KAP1 interacted with SOCS1 (Figure 7D). Next, we used a GST-pull-down assays with two fragments of KAP1 and an in vitro translated and ³⁵S labelled SOCS1. We found that the Bromo/PHD domain of KAP1 directly interacted with SOCS1 (Figure 7E-G). We thus propose that SOCS1 activates the senescence functions of p53 via two mechanisms (Figure 7H): 1) facilitating serine 15 phosphorylation of p53 and 2) p53 stabilization by interfering with KAP1. Together, these mechanisms can contribute to convert an acute p53 response into a chronic response characterized by the expression of a unique pattern of p53 target genes.

Cell lines, reagents, growth analysis and senescence

U2OS were purchased from the American Type Culture Collection, (Manassas, VA) and normal human diploid fibroblasts IMR90 were purchased from the Coriell Institute (Camden, New Jersey, USA). IMR90 were cultured in DMEM supplemented with 10% fetal bovine serum (FBS; Wisent, Montréal, QC, Canada) and 1% penicillin G/streptomycin sulphate. U2OS were supplemented with 5% FBS (Wisent) and 5% Newborn Calf serum (Wisent), with 1% penicillin G/streptomycin sulphate and with 2mM L-glutamine. Tert-butyl Hydroperoxide was purchased from Sigma (cat #458139) and used at 350 µM in U2OS cells and 88 µM in IMR90 cells. Doxorubicin was purchased from Sigma. Growth curves and the senescence associated β-galactosidase were performed as previously described [49].

Plasmid constructions and viral gene transfer

Ca-STAT5A was previously described in [50]. pLPC-SOCS1 was previously described in [14]. pLPC PML IV and pBabe RASV12 were previously described in [49]. GST-KAP1 constructs were a kind gift of Dr. Xavier Mascle. Lentiviral shRNAs against SOCS1 were purchased from Sigma-Aldrich in the pLKO vector. ShSOCS1a (shS1) (TRCN0000356245, Sigma) has the following sequence: CCGGCTGGTTGTTGTAGCAGCTTAACTCGAGTTAAGCTGCTACAACAACCAGTTTTTG and shSOCS1b (TRCN0000057065, Sigma) has the following sequence: CCGGCTTCCGCACATTCCGTTCGCACTCGAGTGCGAACGGAATGTGCGGAAGTTTTTG. shNTC has the following sequence: CCGGCAACAAGATGAAGAGCACCAACTCGAG TTGGTGCTCTTCATCTTGTTGTTTTT. pLXSN E7 was a kind gift from Dr. D. Galloway. Retroviral gene transduction was performed as previously described [49]. Lentiviral gene transduction was performed by co-transfecting 6 µg of the lentiviral vector with the packaging vectors in 293T cells as follows: 3 µg of the VSV-G envelope protein expression plasmid pMD2, 1.5 µg of the regulator of virion expression (REV) expression plasmid pRSV and 1.5 µg of gag/pol elements expression plasmid pMDLg/pRRE. 24 hours post transfection, media were changed and the supernatants (viral soup) were collected 48 hours post transfection, filtered through 0.45 µm filters (Sarstedt) and added to target cells together with DMEM and 4 µg/mL polybrene. Viruses were removed after 8 hours and replaced by fresh medium. Selection was started 24 hours post-infection. For triple infections in IMR90 cells, G418 (300 µg/mL, 6 days) and puromycin (1.5 µg/mL, 3 days) were used first. After puromycin selection was over, hygromycin (80 µg/mL) was added for 4 days. Cells were selected for a total of 7 days before RNA collection. For U2OS and IMR90 single infections, puromycin was used at 3 µg/mL for 3 days.

Ferroptosis cell death assays

U2OS and IMR90 cells were seeded 24 hours prior to treatments in 10 cm plates (Corning) at 60% confluence. Cells were then treated with 350 µM (U2OS) or 88 µM (IMR90) TBH (Sigma #458139) for 16 hours. Cell death was assessed by imaging the cells under a white-light microscope and by counting live and dead cells using Trypan-Blue (BioRad) and a cell counter. Supernatants were collected and added to trypsinized cells before counting. All assays were performed at least three times. For IC50 assays, 20 000 IMR90 cells of each condition were plated per well of a 48 well plate. Cells were plated 24 hours prior to TBH treatment. Cells were then treated with either: 0, 10, 20, 40, 60, 80, 90, 100, 120, 160, 180 or 200 µM of TBH for 16 hours. Cells were then fixed with 1% gluteraldehyde in PBS and stained with Crystal Violet (Sigma #C0775). The dye was then resuspended in 10% acetic acid and dosed with a spectrophotometer. GraphPad Prism was used to generate IC50 curves and determine IC50 values.

Western blotting

Protein analysis was performed by lysing cells in Cell Lysis Buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 2.5 mM sodium pyrophosphate, 1 mM β-glycerolphosphate and a cocktail of protease inhibitors (Roche). Quantification of protein content was performed with the Bradford method. Extracts were prepared with Laemmli buffer (1X final). For endogenous SOCS1 detection, 150 µg of total extract was loaded on a 12% SDS-PAGE gel and transferred on an Immobilon P membrane (Millipore). For SOCS1 overexpression and other protein detection, 50 µg of total cells extract was loaded on a 12% SDS polyacrylamide gel. Western blots were performed as described previously [14]. The following primary antibodies were used: anti-SOCS1, clone 4H1 MBL (cat #K0175-6) used at 1:1000 dilution overnight at 4 °C. Anti-phospho-p53 (Serine15), NEB (cat #9284) used at 1:1000 dilution overnight at 4 °C. Total p53 (DO1 clone), Santa Cruz biotechnology (cat# sc-126) used at 1:1000 overnight at 4 °C. Anti-MCM6, Bethyl (#A300-194A) used at 1:1000 overnight at 4°C. Anti-phosphorylated-Rb (Serine 795), Cell Signaling (#9301) used at 1:1000 overnight at 4°C. Anti-phosphorylated Histone H3 (Serine 10), Millipore (#6570) used at 1:1000 overnight at 4°C. Alpha Tubulin antibody clone B-5-1-2, Sigma (cat # T6074) used at 1:20,000 1 hour at Room Temperature and served as loading control. Signals were revealed by using secondary antibodies coupled to peroxidase (BioRad laboratories) and ECL (GE Healthcare, cat# RPN2106) or Clarity ECL (BioRad #cat 1705061).

Co-immunoprecipitation

U2OS cells were collected in Cell Lysis Buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 2.5 mM sodium pyrophosphate, 1 mM β-glycerolphosphate and a cocktail of protease inhibitors (Roche)) and protein concentration measured using the Bradford method. 2 mg of cell extract from each condition was used for IP and 50 ug was loaded as input (Whole Cell Lysate). Immunoprecipitation of either KAP1 (Bethyl #A300-274A) or SOCS1 (SantaCruz #sc-7005R) were performed at an antibody dilution of 1:200 overnight at 4°C. Protein A Sepharose 4B fast flow beads (Sigma #P9424) were used to immunoprecipitate antibody-protein complexes for 1 hour at 4°C. Beads were then washed four times with Cell Lysis Buffer and complexes were eluted by adding Laemmli 2X buffer directly to the beads. Samples were heated for 5 minutes at 100°C prior to loading on an SDS-PAGE gel for western blotting.

GST-pull-down

BL21 E. coli strain harboring each of the KAP1-GST-fusion vectors were grown at 37°C to an OD of 0.6 in 400 ml YTA 2X medium (16g/l tryptone, 10g/l yeast extract, 5g/l NaCl) and induced with 0.4 mM IPTG for 4 hr at 30°C. Cell pellets were resuspended in STE buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA, 100 mM NaCl) supplemented with DTT (5mM) and a Protease Inhibitor Cocktail (Roche). Cells were lysed by adding 1mg/mL lysozyme and incubating 45 minutes on ice prior to sonicating with a microtip five times for 10 s. Extracts were centrifuged (13 000 rpm 10 minutes) and supernatants were incubated for 2 hr at 4°C with glutathione Sepharose-4B (Amersham) and washed four times with NETN buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA, 100 mM NaCl, 0.5% NP-40). Beads were resuspended in PBS after the last wash.

Briefly, ³⁵[S] Flag-SOCS1 was produced by in vitro transcription/translation (Promega #L1170) and incubated with GST, GST-KAP1, GST KAP1 PHD/BROMO or GST-KAP1 ΔN purified proteins bound to glutathione Sepharose-4B for 2 hr at 4°C. Equal amounts of GST fusion proteins were used as judged from SDS gel electrophoresis and Coomassie blue staining. Beads were washed 5 times with NETN buffer. Precipitates were eluted in 20 µl of SDS Sample Buffer and 10 µl of Bromophenol Blue and boiled for 5 min. Eluates (30 µl) were separated on SDS-PAGE, dried and exposed for autoradiography.

GSH quantification

Cells were collected in Cell Lysis Buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 2.5 mM sodium pyrophosphate, 1 mM β-glycerolphosphate and a cocktail of protease inhibitors (Roche)). GSH was then quantified with the QuantiChrom Glutathione (GSH) Assay Kit (BioAssay Systems, DIGT-250) following manufacturer’s instructions. Quantification was normalized to total protein concentration of each sample.

Real-Time Quantitative PCR (qPCR)

Total RNA was isolated using Trizol (Invitrogen). Reverse transcription was performed using 5X All-In-One RT MasterMix (Abmgood) on 2 µg of total RNA in 20 µL final volume according to the kit’s instructions. Reverse transcription products were diluted 10-fold in RNAse free water before proceeding to qPCR. QPCR was performed using primers and probe sets from Roche Universal Probe Library (https://lifescience.roche.com/en_ca/brands/universal-probe-library.html?_ga=1.38905443.192324701.1470126343#overview). 96-well plate formats with SYBR-green technology was used as described previously [26]. Relative target-gene quantification was obtained by using the ΔΔCT method in a light Cycler 480 (Roche). The mRNA expressions were measured relative to the mRNAs of two housekeeping genes: HMBS and TBP. All forward and reverse primers are listed in Table 4.

Microarray analysis

RNA was collected from IMR90 cells expressing oncoprotein E7, seven days after co-infection with pWZL ca-STAT5A (5A) and pLKO shNTC or pWZL 5A and pLKO shSOCS1 a. Total RNA samples were isolated with the RNeasy mini kit (Qiagen) and sent to the Genome Quebec facility at McGill University for cRNA amplification and subsequent hybridization on GeneChIP Human Gene 2.0 ST Array Affymetrix DNA Chip. Data were analyzed using Affymetrix Expression Console Software and Transcriptome Analysis Console (www.affymetrix.com). Each condition was analysed in biological triplicates and the cut-off applied for analysis was an ANOVA p-value <0.05 for genes that have a fold change ≥±1.5. Data are available at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98216 Further analysis were conducted with DiRE (https://dire.dcode.org/), Gene-Set Enrichment Analysis (GSEA) (http://software.broadinstitute.org/gsea/index.jsp) and DAVID database (https://david.ncifcrf.gov/). Genes identified by DiRE were further analysed with Oncomine (https://www.oncomine.org/resource/login.html).

Gene expression analysis in hepatocellular carcinoma (HCC) specimens

The correlation between SOCS1 gene expression and that of p53 target genes related to ferroptosis was analysed using The Cancer Genome Atlas (TCGA) [51] provisional dataset containing 373 hepatocellular carcinoma (HCC) specimens, which was accessed via the cBioportal (http://www.cbioportal.org) [52]. The transcript levels were expressed as RNASeq Ve RESM (RNA-Seq by Expectation Maximization) [53]. The downloaded data were plotted using the GraphPad Prism software to determine the Spearman correlation (ρ) and statistical significance (p).