biomedicines Article The Histone Variant MacroH2A1 Impacts Circadian Gene Expression and Cell Phenotype in an In Vitro Model of Hepatocellular Carcinoma Annalucia Carbone 1,† , Elisabetta De Santis 2,† , Olga Cela 3 , Vincenzo Giambra 2 , Luca Miele 4 , Giuseppe Marrone 4 , Antonio Grieco 4 , Marcus Buschbeck 5 , Nazzareno Capitanio 3,‡ , Tommaso Mazza 6,‡ and Gianluigi Mazzoccoli 1, *,‡ 1 Department of Medical Sciences, Division of Internal Medicine and Chronobiology Laboratory, Fondazione IRCCS Casa Sollievo della Sofferenza, 71013 San Giovanni Rotondo, Italy; annalucia.carbone@gmail.com 2 Institute for Stem Cell Biology, Regenerative Medicine and Innovative Therapies (ISBReMIT), Fondazione IRCCS Casa Sollievo della Sofferenza, 71013 San Giovanni Rotondo, Italy; e.desantis@operapadrepio.it (E.D.S.); v.giambra@operapadrepio.it (V.G.) 3 Department of Clinical and Experimental Medicine, University of Foggia, 71100 Foggia, Italy; olga.cela@unifg.it (O.C.); nazzareno.capitanio@unifg.it (N.C.) 4 Fondazione Policlinico Universitario A. Gemelli-IRCCS, Catholic University of the Sacred Heart, 00168 Rome, Italy; luca.miele@policlinicogemelli.it (L.M.); giuseppe.marrone@policlinicogemelli.it (G.M.); antonio.grieco@unicatt.it (A.G.) 5 Josep Carreras Leukaemia Research Institute, IJC Building, Can Ruti Campus Ctra de Can Ruti, Camí de les   Escoles s/n, 08916 Badalona, Spain; mbuschbeck@carrerasresearch.org 6 Bioinformatics Unit, Fondazione IRCCS Casa Sollievo della Sofferenza, 71013 San Giovanni Rotondo, Italy; Citation: Carbone, A.; De Santis, E.; t.mazza@css-mendel.it Cela, O.; Giambra, V.; Miele, L.; * Correspondence: g.mazzoccoli@operapadrepio.it; Tel./Fax: +39-(0882)-410-255 Marrone, G.; Grieco, A.; Buschbeck, † These younger authors contributed equally to this article. M.; Capitanio, N.; Mazza, T.; et al. ‡ These senior authors contributed equally to this article. The Histone Variant MacroH2A1 Impacts Circadian Gene Expression Abstract: Hepatocellular carcinoma (HCC) is a leading cause of cancer-related death worldwide. and Cell Phenotype in an In Vitro A foremost risk factor for HCC is obesity/metabolic syndrome-related non-alcoholic fatty liver Model of Hepatocellular Carcinoma. disease (NAFLD) and non-alcoholic steatohepatitis (NASH), which is prompted by remarkable Biomedicines 2021, 9, 1057. https:// changes in transcription patterns of genes enriching metabolic, immune/inflammatory, and circadian doi.org/10.3390/biomedicines9081057 pathways. Epigenetic mechanisms play a role in NAFLD-associated HCC, and macroH2A1, a variant Academic Editor: Massimo Moro of histone H2A, is involved in the pathogenesis modulating the expression of oncogenes and/or tumor suppressor genes and interacting with SIRT1, which crucially impacts the circadian clock Received: 4 August 2021 circuitry. Hence, we aimed to appraise if and how macroH2A1 regulated the expression patterns of Accepted: 19 August 2021 circadian genes in the setting of NAFLD-associated HCC. We took advantage of an in vitro model Published: 20 August 2021 of liver cancer represented by HepG2 (human hepatocarcinoma) cells stably knocked down for macroH2A1 and conducted whole transcriptome profiling and deep phenotyping analysis. We Publisher’s Note: MDPI stays neutral found up-regulation of PER1 along with several deregulated circadian genes, enriching several with regard to jurisdictional claims in important pathways and functions related to cancer onset and progression, such as epithelial- published maps and institutional affil- to-mesenchymal transition, cell cycle deregulation, and DNA damage. PER1 silencing partially iations. mitigated the malignant phenotype induced by the loss of macroH2A1 in HCC cells. In conclusion, our findings suggest a modulatory role for the core circadian protein PER1 in liver carcinogenesis in the context of a lack of the macroH2A1 epigenetic and transcriptional landscape. Copyright: © 2021 by the authors. Keywords: circadian; cancer; biological clock; HCC; PER1 Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons 1. Introduction Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ Liver cancer, the most common form of which is hepatocellular carcinoma (HCC), is 4.0/). a leading cause of death worldwide [1,2]. Therapeutic strategies to treat advanced HCC Biomedicines 2021, 9, 1057. https://doi.org/10.3390/biomedicines9081057 https://www.mdpi.com/journal/biomedicines Biomedicines 2021, 9, 1057 2 of 21 provide only marginal effects [3]. Obesity/metabolic syndrome, currently a pandemic, is a risk factor mainly for HCC among different types of cancer. In fact, obese individuals with a high body mass index (BMI 35–40) had the greatest increase in cancer risk (the risk of HCC was increased by ~5 times) [4]. Obesity is accompanied in more than 90% of cases by non-alcoholic fatty liver disease (NAFLD). NAFLD is the result of an imbalance between the availability of lipids (deriving from the uptake of circulating lipids or from de-novo lipogenesis) and the elimination of lipids (through the oxidation of free fatty acids or the secretion of lipoproteins rich in triglycerides). Consequently, visible lipid droplets accumulate in the hepatocytes. The natural history of NAFLD in humans is not easy to establish since most studies are retrospective; however, there is strong and convincing evidence that patients affected by NAFLD have a higher incidence of HCC. Subjects with NAFLD progress to non-alcoholic steatohepatitis (NASH) in 10% of the cases and eventually to HCC [5]. The high amount of fat in the liver can induce the development of HCC through inflammation, activating specific signal pathways such as those of NF-kB, STAT3/5, JNK, and SIRT1, and consequently the expression of a set of growth factors and cytokines (IL-6, IL-10, TNFα, and lymphotoxins) [3,4]. Indeed, alterations in hepatocyte metabolism and proliferation during steatosis and HCC are triggered by dramatic changes in gene tran- scription patterns. Metabolic pathways, inflammatory processes, and immune response are rhythmically driven by the biological clocks, light-entrainable intracellular mechanisms endowed in all cells and tissues, capable of self-sustaining oscillation, and committed to driving behavioral and physiological rhythms in the organism through hierarchical net- working [6]. These molecular clockworks comprise transcriptional/translational feedback loops operated by core circadian genes and proteins (Arntl, Clock, Per1-2, Cry1-2), which drive rhythmic expression changes in thousands of downstream genes enriching pathways that control critical cell functions [7–9]. The biological clock plays a critical role in NAFLD onset and progression to NASH and HCC, impacting basic cellular processes, such as autophagy, cell cycle, DNA dam- age response, xenobiotic detoxification, reactive oxygen species production, and inac- tivation [10,11]; however, the epigenetic mechanisms involved in HCC and associated with obesity/metabolic syndrome/NAFLD have not been investigated in detail. DNA is wrapped around histone octamer’s backbone to the chromatin structure. The packaging of chromatin in the nucleus is regulated at different levels to allow transcriptional plasticity. One of these levels is the replacement of canonical histones (H2A, H2B, H3, and H4) with histone variants [12,13]. A variant of histone H2A is represented by macroH2A1, which is encoded by H2AFY, a gene present in both sex chromosomes and autosomes. Chromatin domains containing the macroH2A1 variant can be hundreds of kilobases long and occupy approximately a quarter of the human genome [14–16]. MacroH2A1 is thought to act as a strong transcriptional modulator that can repress transcription and activate some genes in response to growth signals that are not yet well defined [17,18]. Mice genetically ablated for macroH2A1 are hallmarked by liver steatosis and glucose and lipid metabolism alterations related to the ge- netic background [19,20]. MacroH2A1 is therefore considered a new factor involved in the pathogenesis of NAFLD. Interestingly, wild-type mice fed a lipogenic diet poor in methyl groups (methionine/choline-deficient diet), which induces liver steatosis and inflamma- tion, show an increase in the total liver content of macroH2A1 [21]. Hence, the lack of macroH2A1 associated with chromatin impacts the complex interaction between enzymes and signal cascades that regulate the intracellular lipid turnover of hepatocytes [20]. MacroH2A1 is present in two isoforms, macroH2A1.1 and macroH2A1.2, generated following alternative RNA splicing. The expression of the macroH2A1.1 and macroH2A1.2 isoforms is predictive of the onset of lung cancer, breast cancer, and gastrointestinal cancers identifying this histone as a valuable marker for risk stratification in cancer patients [22]. The macroH2A1.1 isoform is expressed at high levels in cells that undergo senescence, a Biomedicines 2021, 9, 1057 3 of 21 debated anticancer mechanism, suggesting that macroH2A1.1 may be a useful biomarker of senescent tumor cells [13,18]. The regulation of the expression of oncogenes and/or tumor suppressors by macroH2A1 in hepatocytes could be particularly relevant for the development of HCC associated with metabolic syndrome since the activity of these genes often mechanically connects NAFLD to the onset of HCC; therefore, macroH2A1 isoforms, in addition to being proposed as markers for lung, breast, and skin cancer, could play a key role in the pathogenesis of liver cancer as well [22–24]. As observed in macroH2A1 KO mice, NAFLD is associated with the development of a spectrum of liver damage, including carcinogenesis, and steatosis- associated HCC is histologically featured by immunopositivity for histone macroH2A1 isoforms marks [25]. Importantly, histone macroH2A1 deletion in HCC cells sets out CD4(+)CD25(+)FoxP3(+) regulatory T cells activation and upholds paracrine-mediated chemoresistance [26]. MacroH2A1.1 binds a small metabolite called O-acetyl-ADP-ribose (OAADPR) to high affinity, which macroH2A1.2 does not bind. OAADPR is generated by an enzymatic reac- tion catalyzed by SIRT1, a nicotinamide-adenine-dinucleotide (NAD)-dependent histone and protein deacetylase, whose activation is considered protective against cell metabolic and oxidative stress, along with being considered protective against aging [27]. SIRT1 is frequently found over-expressed in HCC and in this cancer type preserves liver cancer stem cells and their tumorigenic and self-renewal features, propping up tumor recurrence, metastatic spreading, and resistance to therapy [28,29]. Interestingly, liver-specific SIRT1 transgenic mice are protected against the development of HCC associated with metabolic syndrome [30], and SIRT1 is crucial for the functioning of the biological clock [31–33]. On these premises, and considering that macroH2A1 is a remarkable transcriptional modifier, we aimed to evaluate if and how this histone variant affected the expression pattern of cir- cadian genes in the setting of HCC. We addressed this issue, taking advantage of an in vitro model of malignant hepatic tumor represented by HepG2 (human hepatocarcinoma) cells stably knocked down for macroH2A1 through short hairpin RNA and analyzed by RNA- Seq for whole transcriptome profiling and by deep cellular phenotyping. Several circadian genes were found deregulated, in particular, the core clock gene PER1 was up-regulated, with enrichment of several important pathways and functions related to carcinogenesis, such as epithelial-to-mesenchymal transition (EMT), altered cell cycle control, and DNA damage. PER1 silencing to some extent ameliorated HCC cell phenotype in the epigenetic and transcriptional context spawned by loss of macroH2A1 histone variant. 2. Materials and Methods 2.1. Primary Datasets Bioinformatics analyses were performed on two publicly available genomic data: GSE117459 and GSE11923. The GSE117459 dataset renders the RNA-Seq profiling of the whole transcriptome of HepG2 cells knocked down for macroH2A1 [24]. Among the gene symbols in the dataset, we pointed out the genes included in the GSE11923 dataset, rendering high-temporal resolution profiling of C57BL/6J mouse liver samples through Affymetrix arrays. Circadian (24 ± 4 h) transcripts were categorized using Fisher’s G-test at a false-discovery rate of <0.05 and COSOPT [34]. Array probe IDs/nucleotide sequences of 24-h oscillating genes were listed using BioDBnet (https://biodbnet-abcc.ncifcrf.gov/db/db2db.php, accessed on 4 June 2021) and 1842 Ensembl Transcript IDs were recovered from the latter dataset. 2.2. Cell Culture The HepG2 parental cell line, acquired from the American Type Culture Collection (ATCC, Rockville, MD, USA) and the HepG2 cell line with macroH2A1 knock out (KO) were cultured in Dulbecco’s Modified Eagle Medium (DMEM) high glucose (Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS) (Thermo Fisher Scientific, Waltham, MA, USA) and 1% penicillin/streptomycin (P/S) (Thermo Biomedicines 2021, 9, 1057 4 of 21 Fisher Scientific, Waltham, MA, USA). Cell cultures were grown in a 5% CO2 incubator at 37 ◦ C. 2.3. Generation of Lentiviral Vectors As previously described, lentiviral production and transduction of HepG2 cells were performed to stably suppress macroH2A1 expression through shRNAs [23,24,26]. 2.4. RNA-Seq For RNA-Seq, total RNA was extracted from control and macroH2A1 KD HepG2 cell lines with TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA). Indexed libraries were prepared from 2 mg/ea purified RNA with the TruSeq Total Stranded RNA Sample Prep Kit (Illumina, Cambridge, UK) according to the manufacturer’s instructions. Libraries were quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and pooled so that each index-tagged sample was present in equimolar amounts; the final concentration of the pooled samples was 2 nmol/L. Pooled samples were then subjected to cluster generation and sequencing using an Illumina HiSeq 2500 System (Illumina, Cambridge, UK) in a 2 × 100 paired-end format at a final concentration of 8 pmol/L. Short reads were aligned against the hg19 genome assembly using STAR (ver. 2.5.1a). Piled-up reads were counted with the htseq-count tool. Normalization of reads counts and their comparisons were performed using the EdgeR-based R pipeline. Regarding the assessment of differentially expressed circadian genes for macroH2A1 knockout (KO) versus control (CTL) cells, statistically significant differences between groups in gene expression were considered using a cut-off of >±1.5 for the fold change and <0.05 for the adjusted p-value ≤ 0.05. 2.5. Analysis for Functional Categorization of the Most Significant Pathways Enriched by Deregulated Circadian Genes Functional and pathway analyses were conducted using Ingenuity Pathway Anal- ysis (IPA; QIAGEN, Redwood City, CA, USA; www.qiagen.com/ingenuity accessed on 4 June 2021) and specifically considering 1842 Ensembl Transcript IDs recovered from the GSE117459 dataset rendering the RNA-Seq profiling of the whole transcriptome of HepG2 cells knocked down for macroH2A1 and categorized as circadian (oscillating with 24 ± 4 h periodicity) using Fisher’s G-test at a false-discovery rate of <0.05 and COSOPT [34]. Ar- ray probe IDs/nucleotide sequences of 24 h oscillating genes were listed using BioDBnet (https://biodbnet-abcc.ncifcrf.gov/db/db2db.php accessed on 4 June 2021). Both enrich- ment analyses were based on the prior calculation of the activation z-scores, which infer the activation states of predicted biological diseases/functions and pathways. An enrichment score (Fisher’s exact test, p-value) was calculated to measure the overlap between observed and predicted regulated gene sets. We considered p-values < 0.05, z-scores > 1 (minimum activation threshold), and z-scores < −1 (minimum inhibition threshold) as significant. 2.6. Knockdown of PER1 by RNA Interference Short interference RNA (siRNA) against PER1 was synthesized by (Sigma-Aldrich, St. Louis, MO, USA). The sequences of the sense and anti-sense PER1 siRNA were 50 -r (CUUUCCAACAGAUCUGUAA) d (TT)-30 and 50 -r (UUACAGAUCUGUUGGAAAG) d (TT)-30 , respectively. The negative control (scrambled) siRNA sequences were as follows: 50 -r (UUCUCCGAACGUGUCACGU) d (TT)-30 and 50 -r (ACGUGACACGUUCGGAGAA) d (TT)-30 . For the siRNA transfection experiments, control and macroH2A1 KD HepG2 cell lines were seeded in 6-well plates (5 × 105 cells per well). After 24 h, cells were transiently transfected with the PER1-specific siRNA diluted in Opti-MEM (Thermo Fisher Scientific, Waltham, MA, USA) using Lipofectamine™ RNAiMAX Transfection Reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol. After transfection, the cells were incubated for 72 h and subjected to various analyses. Biomedicines 2021, 9, 1057 5 of 21 2.7. Western Blotting Cells transfected with siRNA were lysed using RIPA buffer supplemented with phos- phatase and protease inhibitors (Roche, Basel, Switzerland). The protein amount was determined using a Pierce BCA Protein Assay kit (Thermo Fisher Scientific, Waltham, MA, USA). The obtained proteins were subjected to SDS–PAGE and then transferred to PVDF membranes, that were incubated overnight (ON) at 4 ◦ C with specific antibodies against: PER1 (1:500), PER2 (1:500), Timeless (1:2000), TIPIN (1:1000), and BNIP3 (1:500), which were purchased from Abcam (Cambridge, UK); CRY1 (1:1000), CRY2 (1:1000), which were purchased from NOVUS (Centennial, CO, USA); ARNTL1 (1:1000), CLOCK (1:1000), SirT1 (1:1000), WEE1 (1:1000), c-MYC (1:1000), CDH1 (1:1000), ZEB1 (1:1000), LC3B (1:1000), Atg7 (1:1000), ULK1 (1:1000), and Beclin-1 (1:1000), which were purchased from Cell Signaling (Beverly, MA, USA); GABARAPL1 (1:1000) (Thermo Fisher Scientific, Waltham, MA, USA); C/EBPβ (1:500), β-Actin (1:1000), which were purchased from Santa Cruz Biotechnology (Dallas, TX, USA); p53 (1:1000) (Dako Agilent Pathology Solutions, Carpinteria, CA, USA). Membranes were then incubated for 1 h at RT with goat anti-mouse IgG or goat anti-rabbit IgG secondary antibodies (BIO-RAD, Hercules, CA, USA). Proteins were revealed with Pierce ECL Western Blotting Substrate (Thermo Fisher Scientific, Waltham, MA, USA) and images were acquired using ChemiDoc™ Imaging system (BIO-RAD, Hercules, CA, USA). Protein levels were quantified using the Image Lab software (version 2.1.0.35.deb) (BIO-RAD, Hercules, CA, USA) and normalized to actin. 2.8. Cell Proliferation Assay To evaluate the proliferation ability, HepG2 cells, control, and macroH2A1 KD, trans- fected with negative control and PER1 siRNA, were seeded into 96-well plates in 90 µL complete DMEM medium, at a density of 1 × 104 cells/well, and incubated overnight at 37 ◦ C in 5% CO2 atmosphere to enable cell adhesion. Cell number and density of viable cells were determined at 0, 24, 48, 72, and 96 h using PrestoBlue™ Cell Viability assay (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instruc- tions. Two biological replicates were prepared, and each condition was assayed in triplicate; the results are expressed as mean ± standard deviation (SD). 2.9. Cell Migration Assay Migration was assessed using an in vitro scratch wound healing assay. HepG2 cells with the different conditions were seeded in a 6-well plate at a density of 1 × 106 cells/well for 24 h to achieve 100% confluence followed by starvation in serum-free DMEM medium for 2 h to completely inhibit cell proliferation. After generating a linear scratch wound in the cell monolayer with a sterile 200 µL pipette tip, the culture medium containing detached cells and cellular debris was removed by washing and replaced with a fresh medium. The scratched area was captured at 0 and 72 h after scratching using an Olympus Inverted Phase Contrast Microscope with a 10× phase objective. The width of the gap between the invasion fronts was measured using Fiji software (ImageJ) to calculate the rate of wound closure. The wound closure values refer to the initial position. All experiments were performed at least two times in triplicates. 2.10. Flow Cytometry Analysis 2.10.1. Cell Cycle To analyze cell cycle phase distribution, 105 /mL cells were cultured as previously described in wells of a 12-multiwell plate for each condition to evaluate the incorporation of BrdU after 2 h of in vitro growth. The incorporation of BrdU was carried out with 10 µL of BrdU 10 mg/mL by using APC BrdU Flow Kit (Becton Dickinson, BD, Franklin Lakes, NJ, USA). Subsequently, treated cells were fixed, permeabilized, and treated with DNase to expose incorporated BrdU. Cells were then stained with APC-conjugated anti-BrdU antibody and resuspended in a 7-AAD solution (1×). All FACS data were acquired on Biomedicines 2021, 9, 1057 6 of 21 FACS Calibur, Canto2 (Becton Dickinson, BD, Franklin Lakes, NJ, USA) and analyzed using FlowJo software (Becton Dickinson, BD, Franklin Lakes, NJ, USA). 2.10.2. Half Maximal Inhibitory Concentration (IC50) of Cisplatin To determine cisplatin’s IC50, 105 cells/500 uL were cultured as previously described in 24-multiwell plate and incubated in controlled atmosphere (37 ◦ C, 5% O2 , 5% CO2 , RH > 95%) for 72 h. Cisplatin (Sigma-Aldrich, S. Louis, MO, USA) tested concentrations were 0.1, 0.5, 1, 5, 10, and 50 µmol/L. After 72 h of in vitro growth, cells were detached and washed once in PBS1X 3% FBS. The DAPI fluorescent DNA dye (1:1000 dilution; Sigma-Aldrich) was used to identify live cells. 2.10.3. DNA Damage and Apoptosis Assay with and without Cisplatin Cells for each condition were treated in vitro with 1 µmol/L cisplatin for 72 h to assess DNA damage. Subsequently, treated cells and negative controls were fixed in PFA 4% for 15 min at RT and then permeabilized in ice-cold 90% methanol for 10 min in ice. Cells were washed once in BD PERM/WASH buffer 1× (Becton Dickinson, BD, Franklin Lakes, NJ, USA) to remove methanol and then resuspended in 100 µL of 1:100 diluted Anti-Phospho-Histone H2AX (Ser139) (Cell Signaling, Beverly, MA, USA). After 1 h of primary antibody incubation, cells were stained for 30 min with AlexaFluor 488 conjugated secondary antibody (Thermo Fisher Scientific, Waltham, MA, USA). Cells were then washed once in BD PERM/WASH buffer 1X and resuspended in 400 µL of BD PERM/WASH buffer 1X for flow cytometric analysis. To assess apoptosis, cisplatin-treated cells as reported above and controls were stained following the kit protocol PE Annexin V Apoptosis Detection Kit I (RUO) (Becton Dickinson, BD, Franklin Lakes, NJ, USA). 2.10.4. Mitochondria and Cellular Respiration Measurements Oxygen consumption rates (OCR) were assessed in intact HepG2 cells by Clark electrode (Oxygraph, Hansatech Instruments Ltd., Narborough Rd, Pentney, King’s Lynn PE32 1JL, UK). Briefly, 2–3 × 106 cells were suspended in 0.5 mL of DMEM, and the attained endogenous basal OCR was followed for 5 min; after that, 1.0 µg/mL of the FoF1 ATP-synthase inhibitor oligomycin (OCROligom.) was added, followed by the addition of 1.0 µmol/L of the uncoupler valinomycin (OCRValinom.). Finally, 2 µmol/L antimycin A + 2 µmol/L rotenone (inhibitors of the mitochondrial respiratory chain complex III and complex I, respectively) were added and the attained residual OCR subtracted from the previously recorded OCRs to estimate the mitochondrial-dependent respiratory activity expressed as nmol of O2 consumed/min/106 cells. The OCRValinom. is the maximal respiratory capacity; the difference between basal-OCR and OCROligom. is a measure of the ATP-generating respiration; the OCROligom./basal-OCR ratio is a measure of the proton leak (i.e., the capacity of the mitochondrial inner membrane to dissipate the respiration-mediated electrochemical proton gradient). Where indicated, 4 µM etomoxir, an inhibitor of carnitinepalmitoyltransferase 1a (CPT-1a), was present in the medium from the beginning of the assay. 2.10.5. Statistical Analysis Data are shown as mean ± standard error of the mean (SEM). Group comparison was assessed by one-way ANOVA and Tukey’s honestly significant difference test, using R (R Core Team, 2020). p-values were adjusted using the Benjamini–Hochberg procedure. Statistical significance was declared if p-value ≤ 0.05. 3. Results 3.1. Loss of MacroH2A1 Histone Variants Changes the Expression of Circadian Genes in HCC Cells Using a cut-off of >±1.5 for the fold change and <0.05 for the adjusted p-value, the assessment of differentially expressed circadian genes for macroH2A1 knockout (KO) Among the most up-regulated genes, HSD17B2 (fold change 263.67) encodes 17β- hydroxysteroid dehydrogenase type 2, an enzyme crucial for steroid hormones metabo- lism and for intratumoral estrogens and androgens metabolism (intracrinology) impact- ing tumor progression and patient’s prognosis in sex-steroid-dependent cancers as well Biomedicines 2021, 9, 1057 7 of 21 as in estrogen and androgen receptor-positive liver cancer [35–39]. MPDZ (fold change 76.83) encodes phosphatidylinositol 3-monophosphate (PtdIns3P)-binding protein that dynamically handles key intracellular signaling and operates intracellular second mes- sengers engaged versus control in numerous (CTL) cellular cells showed processes a slight such as cell transcriptional proliferation, divergence betweenmigration, and the matched survival [40]. Deregulation of PI3K-dependent pathways is straightforwardly engaged HepG2 cell lines with 69 differentially expressed genes, 24 up-regulated and 45 down- in carcinogenesis through abnormal activation regulated (Figure 1, Supplementary Table S1). of the Figure 1. Figure Histonevariant 1. Histone variantmacroH2A1 macroH2A1knockdown knockdownaffects affectsthe theexpression expressionlevel levelof ofcircadian circadiangenes. genes. (A) Heat-map rendering mRNA expression levels of 24 h oscillating genes upon histone (A) Heat-map rendering mRNA expression levels of 24 h oscillating genes upon histone variant variant macroH2A1knock macroH2A1 knockdown downinin HepG2 HepG2 cells. cells. GFP1 GFP1andand GFP3 GFP3 indicate indicate HepG2 HepG2 control control cells, cells, KO1 KO1 and and KO2 indicate macroH2A1 knocked-down HepG2 cells (biological replicates). (B) Volcano plot of differentially expressed circadian genes between HepG2 KD and HepG2 control cells. Indicated in orange are the genes with statistically significant different expression, i.e., p-value < 0.05 and abs (FC) > 2. (C) Enrichment analysis for functional categorization of the most significant pathways. Functional and pathway analyses were conducted using Ingenuity Pathway Analysis (IPA; QIAGEN, Redwood City, CA, USA; www.qiagen.com/ingenuity accessed on 4 June 2021). Among the most up-regulated genes, HSD17B2 (fold change 263.67) encodes 17β- hydroxysteroid dehydrogenase type 2, an enzyme crucial for steroid hormones metabolism and for intratumoral estrogens and androgens metabolism (intracrinology) impacting tumor progression and patient’s prognosis in sex-steroid-dependent cancers as well as in estrogen and androgen receptor-positive liver cancer [35–39]. MPDZ (fold change 76.83) encodes phosphatidylinositol 3-monophosphate (PtdIns3P)-binding protein that dynam- ically handles key intracellular signaling and operates intracellular second messengers engaged in numerous cellular processes such as cell proliferation, migration, and sur- vival [40]. Deregulation of PI3K-dependent pathways is straightforwardly engaged in Biomedicines 2021, 9, 1057 8 of 21 carcinogenesis through abnormal activation of the PI3K/AKT/mTORC1 axis [41] and when highly expressed at the protein level in HCC tissue associates with a more favorable prognosis in patients affected by liver cancer [40,42]. MPDZ modulates DLL4-induced Notch signaling in the vasculature through physical interaction with DLL1 and DLL4, assisting their interaction with the adherents junction protein Nectin-2 and improving Notch signaling activity with tumor angiogenesis drop off [43]. The protein encoded by ABCG5 (fold change 16.35) is one of the ATP-binding cassette (ABC) transporters superfamily members, involved in extra and intracellular membrane transport of different molecules and, in particular, of neutral sterol in hepatobiliary and transintestinal cholesterol excretion and its inactivating mutations cause a rare genetic disorder, Sitosterolemia [44]. SGTB (fold change 3.09) encodes small glutamine-rich tetratri- copeptide repeat containing beta, an isoform of SGT, a co-chaperone of heat shock cognate protein of 70 kDa (Hsc70), and its down-regulation by miR-365b enhances HCC cell mi- gration and invasion in vitro and low SGTB tissue expression was found correlated with more favorable HCC progression and prognosis [45,46]. The protein-coding gene SER- PINH1 (fold change 2.80) encodes a member of the serine protease inhibitor superfamily. This serpine is co-expressed with SNHG6 (small nucleolar RNA host gene 6) through the miRNA-139-5p interplay and enhances the in vitro malignant phenotype of HCC cell lines (HepG2, Hep3b, HLE, and Huh-7) [47]. The most down-regulated gene, GNG12 (fold change −115.27), encodes guanine nucleotide-binding protein subunit gamma-12, a G protein family member entailed in the modulation of the inflammatory signaling cascade, which when highly expressed in tumor tissue, enhances cancer cell growth and predicts poor prognosis in pancreatic cancer [48] as well as in liver cancer, in this case in cooperation with GNA12 (fold change −1.79) [42]. Another highly down-regulated gene, CDO1 (FC −60.54), encodes cysteine dioxygenase type 1 and is a tumor suppressor found silenced by promoter methylation in different cancer types, such as non-small cell lung cancer and gastrointestinal cancers, including liver cancer [49–52]. SPARC (secreted protein acidic and rich in cysteine, fold change −36.12) encodes a glycoprotein involved in extracellular matrix remodeling and influencing cell growth and hepatic fibrogenesis in chronic inflammation as well as cell proliferation, migration and invasion, angiogenesis, and tumor progression in HCC, where its decreased expression in tumor tissue predicts poorer prognosis [53–55]. CYP7A1 (fold change −31.38) encodes a member of the cytochrome P450 superfamily of rate-limiting monooxygenases engaged in drug metabolism, lipids/steroids/cholesterol synthesis, and cholesterol to bile acids conversion. Interestingly, bile acids homeostasis disorders along with augmented hepatic bile acids toxic levels associates with higher HCC incidence [56]. RHBDD1 (fold change −29.99) encodes an intramembrane-cleaving serine protease that, similar to other rhomboid proteins, is involved in human cancer progression with a growth-promoting role negatively impacting patients’ survival in different cancer types (e.g., colorectal cancer and breast cancer) [57–59]. Among the core clock genes encoding the circadian proteins that hardwire the molec- ular clockwork, we found up-regulation of PER1 (period circadian regulator 1, fold change 1.5), which is a candidate tumor suppressor gene and, when highly expressed in tumor tissue, is a favorable prognostic biomarker in liver cancer [42,60]. PER1 drives a key con- nection amid the molecular clockwork, the cell cycle machinery, and the DNA damage response (DDR) cooperating with ATM and Chk2 checkpoint proteins [61,62]. When over- expressed, PER1 was shown to set off apoptosis associated with DNA damage and altered expression of critical cell cycle controllers, such as WEE-1, Cyclin B1, and CDC2 in various human breast (MDAMB-231), colon (SW48), lung (NCI-H460), and endometrial (Ishikawa) cancer cell lines [61,62]. 3.2. Circadian Genes Deregulated upon MacroH2A1 Knockout Enrich Important Pathways in HCC Cells Performing functional and pathway analyses with IPA exclusively on 1842 circadian genes sorted out from the GSE117459 dataset (see Materials and Methods), we found sig- nificant enrichment in seven down-regulated pathways (Sirtuin signaling pathway, hepatic Biomedicines 2021, 9, 1057 9 of 21 fibrosis signaling pathway, autophagy, GP6 signaling pathway, ErbB2-ErbB3 signaling, stearate biosynthesis I, fatty acid β-oxidation I, type 2 diabetes mellitus signaling), and three up-regulated pathways (NRF2-mediated oxidative stress response, PTEN signaling, superpathway of methionine degradation). (Figure 1, Supplementary Table S2). Furthermore, we found significantly decreased activation for several functions related to carcinogenesis in the hepatobiliary and gastrointestinal system as well as lipid molecular transport and lipid and steroid synthesis/metabolism. On the contrary, we found increased activation of functions related to cancer aggressiveness and progression, such as epithelial- mesenchymal transition, cell cycle (interphase), DNA fragmentation, and recombination (Supplementary Table S3). Among the significantly down-regulated pathways, the sirtuin signaling cascade is particularly relevant, which interplays both with macroH2A1 epigenetic regulation and circadian clock circuitry functioning. In the liver, SIRT1 mainly controls glucose and lipid metabolism, fine-tuning insulin sensitivity, gluconeogenesis/glycolysis, and fatty acid oxidation/cholesterol metabolism. SIRT1 operates as an intracellular metabolic sensor and, through deacetylation and activation of PGC1-α and PPARα, enhances glucose homeostasis and fatty acids β-oxidation [63]. Further, adjusting the transcriptional activity of CLOCK/BMAL1 heterodimer according to nutrient sensing through intracellular NAD+ levels, SIRT1 takes part in the functioning of the biological clock, and deacetylating PER2 regulates the expression of oscillating tissue-specific genes, joining cellular metabolism to the molecular clockwork [64]. 3.3. MacroH2A1 Histone Knockdown in HCC Alters the Expression of Proteins Entailed in Key Cell Functions Loss of macroH2A1 histone variants determines variations in the expression of dif- ferent circadian proteins and silencing of PER1 seems to unveil that its deregulation may account, at least in part, for these alterations (Figure 2). In particular, the core transcrip- tion factor BMAL1 is increased in macroH2A1 KO cells and further increases upon PER1 silencing, while the expression of CLOCK does not seem to be influenced by the knock- down of macroH2A1 and subsequently of PER1. According to transcriptomics data, the expression of PER1 is increased following knockdown of the macroH2A1 histone variants, and its silencing is accompanied by a striking reduction in the expression of CRY1, CRY2, and PER2, as well as SIRT1, at the protein level. The reduction in the expression of the circadian proteins operating the negative arm of the transcriptional–translational feedback hardwiring of the biological clock could be related to the extremely low level of expression of SIRT1, which plays a fundamental role in the correct functioning of the circadian tran- scriptional process [65]. Knockdown of macroH2A1 histone decreases WEE-1 and, even more, c-MYC protein expression, with further reduction following PER1 silencing. These two proteins link the biological clock to cell cycle control (WEE-1 blocking G2/M transition and c-MYC upholding G0/G1 transition) and their alteration, in the context of macroH2A1 KO-related epigenetic landscape, should be deciphered considering other protein changes to explain the greater replication of macroH2A1 KO cells, in particular in comparison to cells with PER1 silencing [66]. The biological clock intervenes in the coordination of cell cycle progression, DNA replication, and DDR through CRY2 and the interaction of TIMELESS/TIPIN with the fork protection complex on the one hand and with the cell cycle checkpoint proteins Chk1 and the ATR-ATRIP complex on the other hand to organize mitotic kinase activation with DNA replication termination and DNA damage checkpoint response to protect fork integrity and genome stability [67–69]. Protein levels of TIMELESS were slightly increased upon knockdown for macroH2A1 histone, contributing probably to lower DNA damage-induced apoptosis upon cisplatin challenge, an effect that is reverted by PER1 silencing. Opposite patterns were shown for the p53 protein, which is greatly decreased in macroH2A1 KO cells, with partial rescue upon PER1 silencing, and could greatly contribute to the phenotypic changes observed in these cells, particularly regarding the cell cycle and apoptosis. P53 represses the transcription of PER1 and, for its part, PER1, through physical interaction, reduces the stability and transcriptional activity of p53 [70]. tional activity of p53 [70]. The biological clock also drives the circadian oscillation of the proteins involved in the autophagy process by controlling the oscillation of the transcription factors C/EBPβ, an important transactivator of autophagy proteins [71], as well as TFEB/TFE3, driving au- tophagy in relation to nutrient sensing and metabolism [72]. Interestingly, among the eval- Biomedicines 2021, 9, 1057 10 of 21 uated proteins, ULK1 and GABARAPL1 were previously described to oscillate with the circadian pattern, while the expression of ATG7 and BECLIN-1 was not rhythmic [71]. Figure Figure2.2. Western Westernblot blot analysis analysisofof the the expression expression levels levels of of circadian circadian proteins proteinsandandregulatory regulatoryproteins proteins and transcription factors involved in key cell processes in the considered cell conditions. and transcription factors involved in key cell processes in the considered cell conditions. (A)(A) Western West- ern blot blot analysis analysis of theof expression the expression levels levels of proteins of proteins involved involved in the in the molecular molecular clockwork. clockwork. (B) (B) West- Western ern analysis blot blot analysis of theofexpression the expression levels levels of proteins of proteins involvedinvolved in cell processes. in cell processes. (C) blot (C) Western Western blot analysis analysis of of the expression the expression levels of levels of proteins proteins involved in involved in theprocess. the autophagic autophagic process. See the See and Materials the Materials Methods and Methods section for moresection for RDM details. more siRNA details.NEG RDM=siRNA siRNANEG = siRNA scramble-treated scramble-treated HepG2 cells; RDMHepG2 cells; siRNA RDM siRNA PER1 = siRNA PER1-treated HepG2 cells; M1 siRNA NEG = siRNA scramble-treated PER1 = siRNA PER1-treated HepG2 cells; M1 siRNA NEG = siRNA scramble-treated macroH2A1 macroH2A1 knocked down HepG2 cells; M1 siRNA PER1 = siRNA PER1-treated macroH2A1 knocked down HepG2 cells; M1 siRNA PER1 = siRNA PER1-treated macroH2A1 knocked down knocked down HepG2 cells; bars, standard error of mean (SEM); three biological replicates were HepG2 cells; bars, each assayed standard in triplicate anderror of mean results were (SEM); expressedthree as biological replicates mean ± standard errorwere each assayed of mean (SEM). in triplicate and results were expressed as mean ± standard error of mean (SEM). In our hands, some autophagy proteins, such as GABARPL1 and BNIP3, were in- The biological clock also drives the circadian oscillation of the proteins involved in creased, while LC3B, Beclin-1, and in particular ULK and ATG7, were decreased upon the autophagy process by controlling the oscillation of the transcription factors C/EBPβ, knockdown of macroH2A1 histone, with further reduction upon PER1 silencing. These an important transactivator of autophagy proteins [71], as well as TFEB/TFE3, driving latter protein changes suggest that the two different autophagic pathways, the conven- autophagy in relation to nutrient sensing and metabolism [72]. Interestingly, among the tional Atg5-dependent pathway and the alternative Atg5-independent pathway, are re- evaluated proteins, ULK1 and GABARAPL1 were previously described to oscillate with strained. The interpretation of this pattern is difficult, considering that autophagy is a the circadian pattern, while the expression of ATG7 and BECLIN-1 was not rhythmic [71]. double-edged sword in carcinogenesis, on one side providing tumor suppression by re- In our hands, some autophagy proteins, such as GABARPL1 and BNIP3, were in- straining cancer-cell survival and inducing cell death, on the other side supporting the creased, while LC3B, Beclin-1, and in particular ULK and ATG7, were decreased upon proliferationof knockdown ofmacroH2A1 cancer cells, histone, aiding them withinfurther enduring the hostile reduction upon tumor PER1microenvironment silencing. These latter protein changes suggest that the two different autophagic pathways, the conventional Atg5-dependent pathway and the alternative Atg5-independent pathway, are restrained. The interpretation of this pattern is difficult, considering that autophagy is a double-edged sword in carcinogenesis, on one side providing tumor suppression by restraining cancer- cell survival and inducing cell death, on the other side supporting the proliferation of cancer cells, aiding them in enduring the hostile tumor microenvironment and improving growth and aggressiveness [73]. The enhancement of malignant phenotype observed in macroH2A1 knocked down cells could also be related to a decrease in E-cadherin (epithelial, CDH1) in the presence of a high level of the transcription factor ZEB1 (Zinc Finger E-Box Binding Homeobox 1), which hallmarks the activation of the epithelial-to- mesenchymal transition (EMT) program, a long-lasting and reversible process involved in cancer progression and featured by loss of cell polarity and cell–cell adhesion, achieve- ment of migratory and invasive features together with stem cell properties and reduced sensitivity to chemotherapy [74]. Interestingly, upon PER1 silencing, CDH1 increased and ZEB1 decreased in macroH2A1 knocked down cells, suggesting reasonable EMT revert and hinting that this core clock protein plays a key role in determining cell phenotype and fate in the epigenetic background triggered by knockdown of macroH2A1 histone. ment of migratory and invasive features together with stem cell properties and reduced sensitivity to chemotherapy [74]. Interestingly, upon PER1 silencing, CDH1 increased and ZEB1 decreased in macroH2A1 knocked down cells, suggesting reasonable EMT revert and hinting that this core clock protein plays a key role in determining cell phenotype and Biomedicines 2021, 9, 1057 fate in the epigenetic background triggered by knockdown of macroH2A1 histone.11 of 21 3.4. Knockdown of MacroH2A1 Histone Modulates Proliferative and Migratory Capacity of HCC 3.4. Cells Knockdown of MacroH2A1 Histone Modulates Proliferative and Migratory Capacity of HCC Cells The evaluation of cell proliferation did not reveal significant differences between con- trol cells; however, cells The evaluation knocked of cell down for proliferation macroH2A1 did not revealhistone variants significant show slower differences pro- between liferation at 48 and 72 h. Following PER1 silencing, these cells control cells; however, cells knocked down for macroH2A1 histone variants show slowerseem to acquire greater proliferative at proliferation capacity, 48 and 72 although not statistically h. Following significant. PER1 silencing, As regards these cells seem tothe abilitygreater acquire to mi- grate, the picture seems to be the opposite, as following knockdown of macroH2A1 proliferative capacity, although not statistically significant. As regards the ability to migrate, his- tonepicture the variants, the cells seems to beshow greater ability the opposite, to migrate as following when evaluated knockdown as a percentage of macroH2A1 histone of wound area variants, closure the cells showafter 48 h ability greater (p-valueto