Pirin is a prognostic marker of human melanoma that dampens the proliferation of malignant cells by downregulating JARID1B/KDM5B expression

Originally considered to act as a transcriptional co-factor, Pirin has recently been reported to play a role in tumorigenesis and the malignant progression of many tumors. Here, we have analyzed the diagnostic and prognostic value of Pirin expression in the early stages of melanoma, and its role in the biology of melanocytic cells. Pirin expression was analyzed in a total of 314 melanoma biopsies, correlating this feature with the patient’s clinical course. Moreover, PIR downregulated primary melanocytes were analyzed by RNA sequencing, and the data obtained were validated in human melanoma cell lines overexpressing PIR by functional assays. The immunohistochemistry multivariate analysis revealed that early melanomas with stronger Pirin expression were more than twice as likely to develop metastases during the follow-up. Transcriptome analysis of PIR downregulated melanocytes showed a dampening of genes involved in the G1/S transition, cell proliferation, and cell migration. In addition, an in silico approach predicted that JARID1B as a potential transcriptional regulator that lies between PIR and its downstream modulated genes, which was corroborated by co-transfection experiments and functional analysis. Together, the data obtained indicated that Pirin could be a useful marker for the metastatic progression of melanoma and that it participates in the proliferation of melanoma cells by regulating the slow-cycling JARID1B gene.


Pirin expression regulates the proliferation rates of melanoma cell lines.
To further study the role of Pirin in the pathogenesis of melanoma, we evaluated the expression of Pirin in melanocytes and eleven melanoma cell lines by RT-qPCR and in Western Blots. The relative Pirin mRNA levels reflected the generalized weaker expression of the PIR gene in melanoma cell lines relative to the melanocytes from neonatal foreskin (p value = 7.04.e−08: Fig. 2a). Notably, heterogeneous mRNA levels were detected among the melanoma cell lines. Accordingly, when we analyzed the amounts of Pirin protein (Fig. 2b), less Pirin was evident in all the melanoma cell lines analyzed relative to the melanocytes.
These results were consistent with the generally stronger Pirin expression observed in benign nevi relative to the melanoma biopsies. The differences in Pirin expression observed between benign and malignant lesions, as well as in the cell lines, suggested a need to assess the biological processes regulated by Pirin in the context of melanoma. As such, MeWO and A2058 melanoma cells were transfected with a plasmid to overexpress PIR, and, essentially biological processes in cancer progression such as migration and cell proliferation were analyzed.
Overexpression of PIR was verified by Western Blots, achieving a 5 to 15-fold increase in Pirin protein relative to the parental A2058 and MeWO melanoma cell lines, respectively (Fig. 3a, Supplementary file Fig. S2). When the wound-healing and proliferation of these cells was assessed, no significant differences were detected www.nature.com/scientificreports/ in the migration capacity between the control cells and those that overexpress Pirin (Fig. 3b) which was in line with the lack of consistent modifications in migration-sustaining EMT markers (Supplementary file, Fig. S3). By contrast, increased expression of Pirin led to a decrease in the viability of both cell lines, which was evident after 48 h (Fig. 3c). The reduction of viability was due to a decrease on the proliferation rate as determined by CFSE assay (Fig. 3d).
Transcriptome modulation in healthy skin melanocytes after PIR down-regulation. Hav-ing assessed the capacity of Pirin to regulate basic tumor-related processes like proliferation, and in order to reveal the molecular mechanisms regulated by this protein, we performed a RNA-seq analysis. The strong Pirin expression observed in melanocytes and the lack of information regarding the role of Pirin in non-transformed cells led us to select this cell model to downregulate Pirin gene expression and to determine the transcriptomic alterations provoked. Pirin-specific shRNA expression induced a 90% reduction in PIR mRNA and a 30% loss of Pirin protein in HEMn-LP melanocytes (Supplementary file, Fig. S4a-b). Two independent replicates of control  Fig. S4c). This demonstrated the reliability of the data produced and illustrated   4a shows the genes with most variable expression after applying the logFC > (2) and logFC < (− 2) filter).
To validate the sequencing results, we focused on two genes linked to the malignant properties of melanoma cells, both of which underwent a significant increase in their mRNA levels in our dataset (Fig. 4b): Lumican (LUM), a gene related to growth and metastasis [20][21][22] ; and THY-1, a cell surface antigen (Thy-1/CD90) that acts as an adhesion molecule for the extravasation of endothelial tissues and as a tumor suppressor 23,24 . We analyzed the expression of these two genes by RT-qPCR in HEMn-LP PIR-silenced cells, and in MeWO and A2058 metastatic melanoma cells overexpressing Pirin (Fig. 4c). The reduction in PIR mRNA was correlated with higher levels of LUM and THY-1 mRNA. Conversely, overexpression of PIR mRNA was associated with weaker LUM and THY-1 gene expression in transfected melanoma cells relative to the control cells. Together, PIR modulation observed appears to modify the transcriptome of cells, provoking the up and down-regulation of different genes.
To identify the biological processes that might be modified in response to PIR silencing, a functional enrichment analysis was performed on RNA sequencing datasets using the BioJupies Interactive Notebook, obtaining the Biological Processes (BPs) of that Gene Ontology (GO) analysis (Fig. 4d). Among the upregulated genes the BPs were mainly associated with extracellular matrix (ECM) organization, the regulation of migration, the regulation of proliferation and the cell cycle, and cell responses to type II interferon. For the downregulated genes, cell differentiation, cytokine-mediated signaling and innate immune response were the BPs predominantly altered. The GO annotation results (overlapping genes, p value and FDR) of the DEGs after PIR silencing are summarized in Table S1, and the top ten pathways in the pathway-enrichment analysis included: cytokine-cytokine receptor interaction, PI3K-AKT signaling pathway, focal adhesion, MAPK-signaling pathway, and JAK-STAT-signaling pathway (Fig. 4e). PI3K-AKT and cytokine-cytokine receptor interaction pathways were the most significantly enriched terms, both with an FDR = 9.00e−26.
JARID1B a transcription factor putatively regulated by PIR. Due to the fact that Pirin acts as a transcriptional co-regulator with NF1 and NFkB 25 , we set out to determine whether Pirin modulated the expression of genes targeted by other transcription factors in melanocytes. A computational analysis of transcription factors likely to be associated with Pirin and its target genes was performed using the Biojupies Interactive Notebook, which significantly predicted that JARID1B could be one such transcriptional regulator working between PIR and its downstream modulated genes. JARID1B was recently described to be a slow cycling gene involved in the epigenetic regulation and malignant reprogramming of melanoma cells [26][27][28] . In our data, more than 100 DEGs were putative targets of this protein and significantly, the most strongly enriched GO terms for this sub- www.nature.com/scientificreports/ set of genes included regulatory processes related to transcription, proliferation, metabolism, morphogenesis, communication and differentiation ( Fig. 5a and see a summary of the details of these GO entries in Table S2). It is also interesting to note that an association between JARID1B regulation and the PI3K-AKT pathway has been described, which was enriched in response to PIR silencing 29 . Therefore, we wondered whether PIR could somehow modulate JARID1B expression. To address this possibility, melanoma cells were co-transfected with a PIR-coding plasmid (pPIR) or an empty plasmid (pCTR) and a plasmid in which eGFP expression was driven by the JARID1B promoter. Activation of the JARID1B promoter was determined through the proportion of cells expressing GFP and significantly, an increase in Pirin expression in MeWO and A2058 melanoma cells led to a significant decrease in JARID1B promoter activation in both metastatic melanoma cell lines (Fig. 5b). www.nature.com/scientificreports/ Based on our findings regarding the DEGs detected following PIR silencing in melanocytes, and the significant increase in E2F1 (p value 5 × 10−4, FDR 0.022) and c-MYC (p value 7.8 × 10-5, FDR 0.006) expression observed, and given that these are target genes of JARID1B 30,31 , their expression was evaluated in Pirin overexpressing melanoma cells by RT-qPCR. Significantly weaker expression of the JARID1B, E2F1 and c-MYC genes was clearly seen in melanoma cells that overexpressed Pirin relative to the control cells (Fig. 5c). Accordingly, it should be noted that the THY-1 and LUM genes (whose expression is shown in Fig. 4c) are also targeted by JARID1B, and both genes are down-regulated in the metastatic melanoma cells overexpressing PIR. Moreover, the rates of proliferation assayed (Fig. 3c) were in agreement with this molecular analysis, since both c-MYC and E2F1 are proteins involved in cell proliferation and cell cycle regulation. These data are consistent with the downregulation of JARID1B gene expression. Together, these results suggest that PIR might modulate melanoma proliferation by targeting the slow-cycling transcriptional regulator JARID1B.

Discussion
Metastasis represents a significant problem in melanoma, even in early-stage (i.e. stage I-II) patients, since 10% of them relapse within five years after their initial diagnosis 32,33 . This is especially relevant for stage IIB and IIC patients (AJCC 8th edition) who often show worse disease progression than those diagnosed in stage III, in which sentinel node involvement, metastasis in transit or metastasis through microsatellites are detected 34 . In this context, efforts are now being directed towards the early characterization of patients, which is essential in a scenario where the absence of effective therapies for metastatic melanoma requires prognostic markers capable of discerning patients at greater risk of metastasis at early stages of the disease 34,35 . This search for early prognostic factors in melanoma is one of our main research objectives [36][37][38][39][40] , and here we focused on studying Pirin expression as an early diagnostic and prognostic marker of melanoma, as well as its role in the cellular mechanisms driving melanoma progression.
Pirin belongs to the superfamily of cupin proteins and it participates in the regulation of different cellular processes, acting as a protein kinase inhibitor, antioxidant or putative transcriptional co-factor 6,25,[41][42][43] . Increased expression of Pirin has been associated with lymph node metastasis in breast cancer 44 and enhanced expression of this protein has been also observed in colorectal tumors 17 . In melanocyte context, a homogenously high Pirin expression in benign melanocytes from nevi was observed, in contrast to the heterogeneous expression in malignant melanomas. Focusing on stage I and II melanoma patients, we explored whether Pirin expression might be related to a more aggressive phenotype. The patients were divided into disease-free patients and patients who developed metastasis during follow-up. Based on multivariate analyses by Logistic Regression and Cox models, including age and the Breslow index as co-variates, strong Pirin expression was significantly associated with a risk of metastasis, suggesting its importance as a prognostic marker.
Regarding the functional role of Pirin, this protein has been shown to promote cell proliferation and migration in breast cancer models 45 and to support cell migration favoring EMT process in cervix carcinoma cells 14 . Moreover, inhibition or decrease of Pirin in melanoma has been shown to diminish cellular migration capacity of melanoma cells 13,18,25 . Nevertheless, neither proliferation nor migration were affected in colon carcinoma cells 17 . We studied the proliferation and migration of metastatic melanoma cells in which Pirin was overexpressed and we found this upregulation did not modify migration but rather, it did induce a significant decrease in the proliferation rate of both the melanoma cell lines studied. Discrepancies among results may be linked to the tumor type or origin as role of Pirin in EMT-linked increase on migration capacity has been observed in tumors of epithelial origin while melanoma arises from the malignant transformation of neural crest-derived non-epithelial cells 46 . In fact, in our settings, Pirin overexpression did not lead to an alteration on Vimentin or N-Cadherin levels. Lack of effect on migration may also be linked to the methodological approach as migration has been always being affected by the decrease or inhibition of Pirin instead of overexpression-based settings 14,16 . In addition, Pirin has been shown to interact with NF-1, BCL-3 and NF-kB, forming an activated complex able to promote previously mentioned EMT and cell migration 13 . Here, BCL-3 expression in nevi and melanoma biopsies was analyzed in parallel to Pirin expression, yet no relationship between them was found (Supplementary file, Fig. S5). Hence, other cellular signaling pathways would appear to be involved in the metastatic activity of Pirin in melanoma cells. On the other hand, the pro-proliferative role of Pirin has been based on its capacity to support E2F1 expression 45 while our results, in both melanoma cell lines, indicate the opposite effect suggesting a cell-type dependent effect.
The cellular activity of Pirin has mainly been studied in terms of extracellular matrix tumorigenicity 16,25,47 and thus, despite its broad distribution, there is little information regarding the role of Pirin in non-transformed cells and tissues 6,48 . To better understand the role of Pirin in a melanocytic context, we firstly studied Pirin expression in primary melanocytes and melanoma cell lines by RT-qPCR and in Western Blots, demonstrating that primary melanocytes exhibit generally stronger and more homogenous Pirin expression than melanoma cell lines, which had significantly lower expression and more heterogeneity among the different cell lines analyzed.
In accordance with the antiproliferative activity observed when Pirin is overexpressed in melanoma cells, the transcriptomic analysis following PIR-silencing in primary melanocytes here revealed an enrichment of genes involved in the negative regulation of cell proliferation, the G1/S transition and extracellular matrix organization and positive regulation of cell migration (Table S1). Furthermore, deleterious mutations in the PIR gene were recently identified in breast cancers that could affect protein structure, stability and function 45 . These results could explain the discrepancies found when studying different cancers or different tumor cell lines. On the other hand, melanoma heterogeneity was recently proposed to be due to the co-existence of different melanoma cell phenotypes and adaptive phenotype plasticity given that transcriptional reprogramming could drive melanoma progression 49 . Transcriptional reprogramming has been detected at different stages of melanoma, with enhanced mesenchymal traits in circulating melanoma cells and proliferative features in metastatic tumors 49  www.nature.com/scientificreports/ Hence, cells with different phenotypes may interact in a cooperative manner and contribute to successful metastatic progression 50,51 .
In recent years, attention is paid to epigenetic regulation in melanoma 28 , which led to the description of JARID1B as an epigenetic regulator implicated in the transcriptional reprogramming of several tumor cells and in tumor heterogeneity 52 . Although JARID1B expressing melanoma cells represent only a small proportion of the cells in the primary and metastatic melanoma populations 53 , the RNA-seq dataset and the transcription factor enrichment analysis found that JARID1B could target more than 100 of the DEGs identified. Furthermore, cotransfection experiments showed a decrease of JARID1B promoter activation after Pirin overexpression, pointing to a functional relationship between Pirin and JARID1B expression. In addition, we demonstrated that the overexpression of Pirin in both the metastatic melanoma cell lines studied led to a significant decrease in JARID1B gene expression, and that of its target genes E2F1 and c-MYC 30,31 . These results may explain the antiproliferative effect of Pirin observed in melanoma cell lines. Indeed, in canine oral melanoma cell lines JARID1-inhibitors drove anti-proliferative activity and overcame cisplatin resistance 54 .
From our data, we believe that in melanocytes Pirin expression could regulate the rate of proliferation through JARID1B and E2F1 pathway 39,55-57 . Moreover, as JARID1B is involved in the neural differentiation process 58 , we hypothesize that it could be maintaining the undifferentiation state of the low proliferative Pirin expressing melanoma cells and potentially a more aggressive phenotype. Indeed, tumor cells with a slow-cycling phenotype may be metabolically active and highly aggressive, with increased potential to grow and metastasize 27,59,60 .
In summary, our study supports the use of Pirin staining along with the Breslow index as a prognostic marker at early stages (I-II) of melanoma. High Pirin protein levels are associated with a more significant probability of metastasis, as well as a shorter time until this clinical end-point. In addition, we propose that Pirin could play an important role in modulating the differentiation and proliferative state of melanoma cells by regulating JARID1B gene expression. However, further research will be necessary to better understand the mechanisms underlying this phenomenon, which could shed light on useful therapeutic strategies for these tumors. After surgery of the primary tumor, clinical check-ups were scheduled every three months for the first two years of the follow-up, and every six months thereafter, until a five-year follow-up had been completed. Annual revisions were then scheduled up to the tenth year post-surgery. The patients who developed metastasis during the follow-up period were again examined every three months for two years after their metastasis had been diagnosed. The presence or absence of metastasis was assessed in all patients by physical examination, as well as through laboratory and radiological testing (X-rays and/or computed tomography -CT-scanning). Metastases were detected in 92 of the 239 melanoma patients studied (38%), including those in whom the disease had already spread at the moment of diagnosis. Disease stages were classified according to the AJCC 8th edition 19 , and the clinical and diagnostic data for each patient was collected retrospectively from centralized electronic and/or paper medical records. For the statistical prediction analysis, only melanoma patients at early disease stages (I and II) were considered (159 subjects in total), and inclusion in the "disease-free" group required a minimum tracking of 2 years.

Patients.
The study was carried out in accordance with the Helsinki Declaration and it was approved by the Drug Research Ethics Committee of the Department of Health of the Basque Government (reference 16-99) (http:// www. euska di. eus/ comite-etico-inves tigac ion-con-medic ament os/). Written informed consent was obtained from all the subjects and the melanoma biopsies collected were stored at the Basque Biobank until use (https:// www. bioba ncova sco. org/). For more details about the informed consent see supplementary material. The clustergrammer displays the gene expression in for each sample in the RNA-seq dataset. Every row of the heatmap represents a sample and each column represents a gene. Prior to displaying the heatmap, the raw gene counts were normalized using the logCPM method, filtered by selecting the 100 genes with most variable expression and finally transformed using the Z-score method. (b) The figure contains a scatter plot which displays the log2-fold changes and statistical significance of each gene calculated by performing a differential gene expression analysis. Every point in the plot represents a gene. Red points indicate significantly up-regulated genes and blue points indicate down-regulated genes. Gene fold changes were transformed using log2 and displayed on the x axis; P-values were corrected using the Benjamini-Hochberg method, transformed using log10, and displayed on the y axis. (c) Validation of RNA Sequencing data using RT-qPCR o PIR, LUM and THY-1 gene expression in HEMn-LP primary melanocytes with dampened Pirin expression, and in A2058 and MeWO metastatic cell lines transfected with a plasmid overexpressing Pirin. The average of three independent assays is shown and the red line highlights the normalized control expression: *p value < 0.05. (d) Bar charts displaying the results of the Gene Ontology enrichment analysis of Biological process generated using EnrichR. The x axis indicates the gene number for each term, the blue bars represent the terms enriched from the downregulated gene subset and the red bars represent the terms enriched from the upregulated genes subset. (e) Canonical pathways significantly modulated by PIR silencing. Every row of the figure represents one enriched process with an FDR cut-off of < 0.05, and each bar represents the number of genes included in each pathway.  blocks were subjected to antigen retrieval in citrate buffer (pH 6.1) and steam for 105 min, and then analyzed by Pirin immunohistochemistry with an anti-Pirin antibody (PA5-29777: Thermo Fisher Scientific, Waltham, MA, USA) and using the EnvisionTMG|2 Sistema/AP Kit (Dako Corporation, Denmark). The slides were counterstained with hematoxylin and images were obtained using a NanoZoomer S210 Digital slide scanner (Hamamatsu C13239-01). The staining intensity was evaluated as negative, low or high expression after independent examination by two observers. Discordant assessments were reviewed jointly to obtain a conclusive consensus evaluation.
Statistical analysis. Bivariate statistical testing of categorical parameters (e.g. sex and melanoma prognosis) was performed with Pearson's chi-squared tests, whereas the Cochran-Armitage Test was employed when Pirin was involved. Statistical differences between two or more groups were estimated with a two-tailed Kruskal-Wallis test or by ANOVA if normality was assumed based on the Shapiro-Wilk test. A logistic regression model was employed to test the association between Pirin and prognosis, adjusting for possible sources of co-variation. For this analysis, marginal effects were also computed 61 . Finally, a Cox Proportional-Hazard model was established to test the association of Pirin expression with the onset of metastasis (measured in months). All p-values were computed non-parametrically using the R package "coin" 62 , controlling for the false discovery rate (FDR) in the case of multiple testing 63 , and statistical evidence of the effects was quantified using the BIC (Bayesian Information Criterion) and Bayes factors 64 . The Bayes Factors indicate the probability of the alternative hypothesis (in our case, there is an association of Pirin expression with metastasis) relative to the null hypothesis (no association of Pirin expression with metastasis) and vice versa, such that: BF01 = Probability (Null Hypothesis)/ Probability (Alternative Hypothesis); and BF10 = 1/BF01 (the inverse). All the analyses were carried out using R version 3.6.1.

Cell lines, and the proliferation and migration assays.
In the present work, eleven melanoma cell lines and three primary melanocyte cell lines were studied. The primary human melanocytes were purchased from Invitrogen (Carlsbad, CA, USA), and they correspond to lightly (HEMn-LP, C-002-5C), moderately (HEMn-MP, C-102-5C) and darkly (HEMn-DP, C-202-5C) pigmented neonatal foreskin. All primary human The melanoma cell lines were cultured in an appropriate medium supplemented with 10% fetal bovine serum (FBS), L-glutamine (2 mM) and antibiotics, according to the supplier's indications. All the primary human melanocytes and melanoma cell lines were cultured at 37 °C in and atmosphere of 5% CO 2 and 95% humidity. For monitoring cell proliferation, the control melanoma cells and the melanoma cells with stable PIR gene transfected clones were labeled with CFSE, a dye which is retained in the viable cells and gradually halves as the cells division. Briefly, desired density grown cells were incubated in pre-warmed PBS containing CFSE working solution (1 μM). After incubation, cells were detached and fix in paraformaldehyde 4% and analyzed by flow cytometry. The reduction in CFSE fluorescence signal was visualized by computer analysis using WEA-SEL 3.7.1-Flow Cytometry Data Analysis software. Moreover, we carried out the 2,3-bis-(2-methoxy-4-nitro-5-sulfophenyl)-2H-tetrazolium-5-carboxanilide (XTT) assay, following the manufacturer's instructions. Absorbance at 490 nm was measured at 24 and 48 h, and cell viability was calculated in relation to untreated control cells as: (experimental absorbance/untreated control absorbance) × 100.
In migration assays, after seeding in 24-well plates the cell monolayers were incubated with Mitomycin C (0.5 µg/mL) for 2 hours and they were then scraped with a sterile plastic micropipette tip. Wound closure was observed over 48 h, taking photos at 0, 24 and 48 h under a light microscope. . After 15 minutes on ice in lysis buffer, the lysates were cleared by centrifugation at 10,000 g for 5 min and the concentration of the total protein recovered was determined using the bicinchoninic acid assay.

Sequencing analysis.
For RNA Sequencing, the RNA concentration was first measured in a Qubit 2.0 RNA assay Kit (Invitrogen, Carlsbad, CA, USA, #Q32855). The quality of all the RNAs assayed was characterized on an Agilent 2100 Bioanalyser using an Agilent 6000 Nano Chip #5067-1511, with an optimal RNA Integrity Number (RIN) between 8 and 10. Sequencing libraries were prepared following the "TruSeq Stranded mRNA Sample Preparation Guide (Part # 15031058 Rev. E)" with the corresponding kit [Illumina Inc., #RS-122-2101 or RS-122-2102 (Set A or B, respectively)]. Briefly, mRNA was purified from the total RNA (400 ng), fragmented and primed for cDNA synthesis. The first strand was synthesized for 10 min at 25  The raw data have been deposited in the NCBI Sequence Read Archive (SRA) through the Gene Expression Omnibus (BioProject accession number PRJNA843921). The reads were initially mapped to the reference genome from the UCSC genome browser (Human GRCh38/hg38 version) using Hisat2. The reads were assembled into transcripts using StringTie and the transcript read count was quantified using the feature counts software of the Subread package in the R statistical environment. Differentially expressed genes (DEGs) were identified by means of the edgeR package in R-studio, implementing a negative binomial distribution to assess statistical significance. Normalization was performed using the trimmed mean of M values (TMM) method. For each gene, the reads per kilobase of exon model (RPKM) per million mapped reads were transformed to log2-fold changes using R version 3.1.0. The analysis of variance (p < 0.05) and false discovery rate (FDR < 0.05) tests were performed using the R program (version 3.1.0) to select genes exhibiting significantly different expression patterns. The upregulated and downregulated DEG sets were generated by extracting the 500 genes with the respectively highest and lowest values from the gene expression signature. Spearman's rank correlation coefficient analysis (SRCCA) was used to determine the strength of association between samples. Enrichment results were generated by analyzing the upregulated and downregulated gene sets using BioJupies 65 and EnrichR 66 , which are freely available at https:// amp. pharm. mssm. edu/ bioju pies and http:// amp. pharm. mssm. edu/ Enric hr, respectively. A KEGG pathway enrichment analysis of the DEGs was performed using the KOBAS online database (available online: http:// kobas. cbi. pku. edu. cn/).

RT-qPCR analysis.
Total RNA was isolated from the cultured cells using the RNeasy Mini kit (Qiagen Inc, Hilden, Germany) and cDNA was synthesized from the total RNA (1 µg) from each sample using the iScriptTM cDNA Synthesis kit (Bio-Rad, Hercules, CA, USA), according to the manufacturer's instruction. Quantitative real-time RT-PCR (RT-qPCR) assays were carried out on an iCycler PCR platform (Bio-Rad, Hercules, CA, USA). The reaction mixture contained cDNA (0.1 µl) from the reverse transcription reaction, together with forward and reverse specific primers (see Table 2), and the iQTM SYBR ® Green Supermix (Bio-Rad, Hercules, CA, USA) in a final reaction volume of 20 µl.
The PCR began by heating the reaction at 95 °C for 10 min, followed by 45 cycles of denaturation at 95 °C for 30 s, annealing at the corresponding temperature for each gene (56-61 °C) for 20 s and extension at 72 °C for 30 s. Each assay included a negative control (a sample without cDNA). The expression data were generated from two amplification reactions with the samples and controls run in triplicate. Optical data obtained by RT-PCR were analyzed using the MyiQ Single-Color Real-Time PCR Detection System Software v.1.0 (Bio-Rad). Analysis of the melting curve for each PCR assay and 1.5% agarose gel electrophoresis analysis of randomly selected samples was performed to confirm the specificity of the amplification products. The expression of three different housekeeping genes (ACTB, GAPDH and RPS15) was also analyzed to normalize the expression data using the Gene Expression Macro Software Version 1.1 (Bio-Rad Laboratories, Hercules), where the relative expression values were calculated by the comparative Ct method 67,68 .