DNA Methylation Signature in Monozygotic Twins Discordant for Psoriatic Disease

Background: Psoriatic disease is a multifactorial inflammatory condition spanning from skin and nail psoriasis (Pso) to spine and joint involvement characterizing psoriatic arthritis (PsA). Monozygotic twins provide a model to investigate genetic, early life environmental exposure and stochastic influences to complex diseases, mainly mediated by epigenetics. Methods: We performed a genome-wide DNA methylation study on whole blood of monozygotic twins from 7 pairs discordant for Pso/PsA using the Infinium Methylation EPIC array (Illumina). MeDiP—qPCR was used to confirm specific signals. Data were replicated in an independent cohort of seven patients with Pso/PsA and 3 healthy controls. Transcriptomic profiling was performed by RNAsequence on the same 7 monozygotic twin pairs. Results: We identified 2,564 differentially methylated positions between psoriatic disease and controls, corresponding to 1,703 genes, 59% within gene bodies. There were 19 regions with at least two DMPs within 1 kb of distance and significant within-pair Δβ-values (p < 0.005), among them SNX25, BRG1 and SMAD3 genes, all involved in TGF-β signaling pathway, were identified. Co-expression analyses on transcriptome data identified IL-6/JAK/STAT3 and TNF-α pathways as important signaling axes involved in the disease, and they also suggested an altered glucose metabolism in patients’ immune cells, characteristic of pro-inflammatory T lymphocytes. Conclusion: The study suggests the presence of an epigenetic signature in affected individuals, pointing to genes involved in immunological and inflammatory responses. This result is also supported by transcriptome data, that altogether suggest a higher activation state of the immune system, that could promote the disease status.


INTRODUCTION
Psoriatic arthritis (PsA) and psoriasis (Pso) represent different manifestations of the same disease spectrum, i.e., psoriatic disease, a systemic chronic inflammatory condition, where genetic factors contribute to disease susceptibility and where biomarkers are virtually absent (Ritchlin et al., 2017). Monozygotic twins (MZ) exhibit a variable degree of concordance for several complex disorders, including autoimmune diseases (Ceribelli and Selmi, 2020), and twin studies contributed to identify susceptibility genes. MZ twins show concordance rates as high as 64% for psoriasis, while in the case of PsA these are almost equal in MZ and dizygotic (DZ) twins although these rates are difficult to compare, based on the bias related to the coexistence of Pso and PsA (Pedersen et al., 2008). Nonetheless, these observations suggest that epigenetic factors might play a major role (Ceribelli and Selmi, 2020), as supported by the experimental evidence in healthy MZ twins (Brodin et al., 2015;Generali et al., 2017) and by the phenotype changes according to age in MZ twins (Fraga et al., 2005). Further, DNA methylation patterns are associated with MZ twins' discordance in the imprinted regions of monogenic syndromes (Ollikainen and Craig, 2011) and in inflammatory diseases with largely incomplete MZ twin concordance, such as systemic lupus erythematosus, type-1 diabetes, rheumatoid arthritis, and systemic sclerosis (Absher et al., 2013;Coit et al., 2013;Liu et al., 2013;Ramos et al., 2019). Limited data from epigenome-wide association studies (EWAS) are available for Pso and PsA and a deregulated epigenetic machinery has been demonstrated in psoriatic skin and whole blood (Zhou et al., 2016).
We hypothesize that epigenetic changes, in specific DNA methylation, may contribute to the onset of Pso/PsA in discordant MZ twins. (Pai et al., 2015;Schübeler, 2015). We performed a peripheral blood EWAS on a cohort of 14 (seven pairs) MZ twins clinically discordant for Pso/PsA and identified DNA methylation changes with specific differentially methylated positions (DMPs) associated with psoriatic disorder. The characterization of a possible circulating psoriatic diseaseassociated methylome would identify specific candidate loci to be further investigated as molecular biomarkers.

Study Participants
We enrolled pairs of MZ twins previously identified through the collaboration with the Italian Twin Registry-Istituto Superiore di Sanità and residing in the area surrounding the Humanitas Research Hospital in Rozzano (Milan, Italy). Participants were evaluated jointly by a rheumatologist and a dermatologist to define the diagnosis of Pso according to the physical examination and/or PsA according to the CASPAR classification criteria (Taylor et al., 2006). The study was approved by the local IRB and, following the signing of an informed consent, all participants provided peripheral whole blood samples at the time of the joint rheumatological and dermatological evaluations. Whole blood and peripheral blood mononuclear cells (PBMCs) were obtained from both twins at the same time using heparinised vacutainers and processed in parallel. In total, 7 MZ twin pairs discordant for Pso or PsA (the presence of either being defined psoriatic disease) were included. Epigenetic data were validated in an independent cohort of seven patients with Pso/PsA and 3 healthy controls. Disease subtype, current and/or past drug history of the affected twins, comorbidities, family history, smoking status and smoking behaviour, diet, alcohol consumption and allergies of both twins were recorded.

DNA Extraction and Methylation Analysis
Genomic DNA was extracted from whole blood using an automated DNA extractor (Chemagic Star workstation; Hamilton, ON, Canada), following manufacturer's instructions. For each sample, 500 ng of genomic DNA was bisulfite converted using the EZ DNA Methylation Kit (ZymoResearch, Irvine, CA, United States). Bisulfite-treated samples were then processed using the EPIC array (Illumina, San Diego, CA, United States), according to manufacturer's instructions. The iScan system (Illumina) was used to scan the arrays and obtain the raw intensity data files. DNA methylation data were processed using the R package minfi (Aryee et al., 2014), to obtain DMPs. First, a quality control was performed to evaluate samplespecific methylation parameters, exploiting the presence of control probes on the EPIC array. Then, the mean detection p-value was evaluated across all samples, and those probes with a p-value > 0.01 in one or more samples were discarded from the dataset. We also removed those probes where SNPs may affect the CpG, and underperforming probes, as suggested by Zhou and others (Zhou et al., 2017). Finally, we also filtered out probes localized on the X and Y chromosome, to avoid sex-specific effects. After this QC step, the final dataset included 762,451 probes, which accounts for 88% of the EPIC probes. Finally, data were normalized using the minfi quantile normalization algorithm. DMPs were calculated using the R-package limma (Ritchie et al., 2015), the function lmFit, and eBayes, based on a paired statistics, and correcting for batch effects with the sva package (Leek et al., 2020).

Regulatory Annotations
The eFORGE v2.0 (https://eforge.altiusinstitute.org/) (Breeze et al., 2016) tool was used to identify if the associated CpGs (top 1,000) were enriched in cell-specific regulatory elements, such as DNase I hypersensitive sites (DHSs) (active regulatory regions) and loci with overlapping histone modifications (H3K4me1, H3K4me3, H3K9me3, H3K27me3, and H3K36me3) across available cell lines and tissues from the Roadmap Epigenomics Project, BLUEPRINT Epigenome, and ENCODE (Encyclopedia of DNA Elements) consortia data. Chromatin Hidden Markov Models (ChromHMMs) bioinformatic approach was adopted to annotate Pso/PsAassociated DMPs through the combination of epigenetic maps across multiple cell types (Chen et al., 2013).

Gene Ontology Analysis
Enrichment analysis for the host genes of the differentially methylated positions was performed using the "topGO" Bioconductor tool, the "biological processes," the "molecular functions" database, and the elim algorithm (Alexa et al., 2006). Gene ontology enrichment was performed on those genes containing at least two DMPs.

MeDIP-qPCR
MeDIP was performed to detect immunoprecipitated methylated DNA with an anti-5′-methyl-cytosine antibody. The MagMeDIP qPCR kit (Diagenode, Belgium) was used following manufacturer's instructions. Briefly, genomic DNA was sheared using a BioRuptor sonicator (Diagenode, Belgium) to produce 400 bp fragments, which were checked by gel electrophoresis. The obtained fragments were immunocaptured with a monoclonal antibody specific for 5-methyl cytosine (supplied by the kit). Methylated DNA was then washed and purified from beads with a DNA Isolation Buffer (DIB) and proteinase K provided by the kit. Following an incubation at 55°C and at 100°C (for 15 min each), the supernatant containing the DNA was used for qPCR analysis to evaluate enrichment.

RNA Sequencing
Whole blood was collected in a PAXgene Blood RNA Tube (PreAnalytiX, Switzerland) and RNA extraction performed using the Maxwell simplyRNA Blood Kit (Promega, United States), following the manufacturer's instructions. RNA quality was assessed by LabChip GX Touch (PerkinElmer, Waltham, MA). Libraries were prepared starting from 500 ng of total RNA, using the TruSeq Stranded mRNA Library Prep Kit (Illumina) and following the manufacturer's instructions. Samples underwent a paired-end 75 bp sequencing using a NextSeq 500 platform (Illumina).

Data Analysis
Sequencing reads were mapped to the human genome (hg19) using STAR (version 2.5.2). Transcript quantification from mapped reads was performed using HTSeq-count (version 0.6.1p1) and the human transcripts annotations from Ensembl database (GRCh37 version). Differential expression analysis was carried out using R and the DESeq2 package, considering an FDR<0.1 as threshold. The enrichment analysis was performed using the Enrichr online tool (Chen et al., 2013), the MSigDB Hallmark pathway, and the PheWeb databases. The co-expression of genes was evaluated starting from transcriptome data and the Coseq R package (Rau and Maugis-Rabusseau, 2018). We applied a centered log ratio (CLR) transformation to count data before fitting a Gaussian Mixture Model, in order to identify gene clusters. For the enrichment analysis we selected only those genes that were attributed to a cluster with a conditional probability>0.9. For each cluster, the average expression profile of all the genes belonging to the cluster was compared between affected and non-affected MZ twins by the Wilcoxon rank sum test.
The estimation of the cellular composition of the samples was performed using the immunedeconv R package, and the quantiseq method (Sturm et al., 2020).

Quantitative Real-Time Polymerase Chain Reaction
Quantitative RT-PCR was performed following the protocol previously described (Vecellio et al., 2016). RT-PCR was performed in triplicate and the 2−ΔCt method was used to calculate the expression of SMAD3, BRG1 and SNX25 relative to β-Actin used as housekeeping gene (ID Assay qHsaCED0036269, Bio-Rad Laboratories, Kidlington, United Kingdom).

Statistical Analysis
We performed one sample t-test to determine statistically significance in RT-PCR and MeDIP followed by qPCR experiments. Significance values were set to p < 0.05. In Figure 6 and related text, data are represented as mean ± SEM.

Patient and Public Involvement
It was not appropriate or possible to involve patients or the public in the design, or conduct, or reporting, or dissemination plans of our research.

Pso/PsA-Associated Whole Blood Methylome Study
The genome-wide DNA methylation profile was investigated using the Illumina Inflnium MethylationEPIC BeadChips (EPIC arrays). The general demographic and clinical characteristics of the MZ twins for which methylation data were available are illustrated in Table 1. After performing the quality check and the exclusion of low-quality and unreliable probes, sites that overlap with SNPs, probes with a very low signal and those located on the X and Y chromosomes, methylation data for 762,451 sites were obtained. The experimental workflow and the analyses performed are shown in Figure 1.

Identification of DMPs Between Psoriatic Disease Affected and Non-affected Twins
The principal component analysis (PCA) highlighted the high similarity between MZ twins (PC1 vs PC2, Figure 2A). However, when higher components were considered, a clear cluster separation between affected and non-affected twins emerged, suggesting a different methylation pattern related to the disease status (PC7 vs PC8, Figure 2B). The statistical analysis did not demonstrate genome-wide significant probes, but when a more stringent filter was applied to the analysis (p < 0.005), 2,564 DMPs were observed ( Figure 2C), mapping to 1,703 genes ( Figure 2D), in the majority of cases at the level of gene Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 778677 bodies (56%), less frequently in intergenic (27%) and promoter regions (16%), in line with the probe distribution in the array, where 51% of probes map to gene bodies, 29% to intergenic and 20% to promoter regions. When evaluating whether the Pso/PsA-associated DMPs reside in genomic regulatory regions involved in transcriptional regulation (such as gene promoters, enhancers, transcriptional start sites, specific histone modifications) in different cell types and tissues we observed that the Pso/PsA-associated CpGs were enriched in H3K36Me3 specific histone modifications, especially in blood cells, such as monocytes, lymphocytes, NK cells and hematopoietic stem cells (FDR<0.01, Supplementary Figure  S1A). To predict disease-relevant cell types, we investigated the ChromHMM state, which was suggestive of either strong or weak gene transcription overlapping associated DMPs in blood cells ( Figure 3) and epithelial fibroblasts (Supplementary Figure S1B) among others (FDR<0.05, p < 10 −10 ).
We found 19 regions with at least two DMPs within 1 kb, significant within-pair Δβ-values and a suggestive p < 5 × 10 −3 , as shown in Figure 2C. The 19 regions associated with the identified DMPs are listed in the Supplementary Table S1 and

Psoriatic Disease MZ Twins Have a Distinct Transcriptomic Profile
The RNA sequencing analysis, performed on the same cohort of MZ twins didn't highlight differentially expressed genes surviving the multiple testing correction. Nonetheless, we performed an enrichment analysis selecting those genes characterized by a significant unadjusted p-value (N 662). Very interestingly, when considering deregulated genes, we observed an enrichment in the pathways related to "oxidative phosphorylation," "inflammatory response," and "MYC targets" (Supplementary Figure S2, Supplementary Table S2).
As a next step, we searched for gene clusters that were coexpressed with a different profile between discordant MZ twins. The analysis identified 20 clusters ( Figure 5A), the majority showing a statistically significant difference in the average expression profiles between the two groups. We ranked the clusters based on the average expression ratio between affected patients and controls, and we selected the most dysregulated three (clusters 14, 19, 13) to perform an enrichment analysis considering the genes included in each cluster. Interestingly, cluster 19 genes, characterized by a higher average expression profile in affected MZ twins (average expression profile ratio: 1.14, Wilcoxon p < 2.2e-16), were associated to "Psoriasis vulgaris" and "Psoriasis" in the PheWeb database ( Figure 5B). Moreover, the same set of genes showed a trend for enrichment in the "IL-6/JAK/ STAT3" pathway ( Figure 5C), a signaling axis that has been implicated in psoriasis (Andrés et al., 2013). Cluster 13, also showing a higher expression in patients compared to controls (average expression profile ratio: 1.1, Wilcoxon p < 2.2e-16), confirmed the enrichment in genes involved in the "IL-6/JAK/STAT3" pathway, as well as in "inflammatory response" and "TNF-α signaling pathway" ( Figure 5D). Cluster 14, characterized by a higher expression in controls (average expression profile ratio: 0.86, Wilcoxon p < 2.2e-16), was instead mainly enriched in "oxidative phosphorylation" Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 778677 5 pathway ( Figure 5E). Finally, we evaluated through a deconvolution method, the cellular composition of the RNAseq samples, but no significant differences were observed between affected and non-affected twins (data not shown).

Evidence for SNX25 Deregulation in Affected MZ Twins
Among the differentially methylated regions, only SNX25 resulted also nominally deregulated in MZ twins' transcriptome data. This gene is involved in the TGF-β pathway and potentially relevant in Pso and PsA pathophysiology. The observed downregulation (logFC -0.48, p 0.02) is concordant with the methylation profile, with affected MZ twins showing higher methylation levels ( Figure 4B). To confirm this result, a quantitative PCR specific for SNX25 promoter was carried out after immunoprecipitation of the methylated DNA in PBMCs of MZ twins. The analysis showed the same trend, although not significant, observed in the methylation array ( Figure 6A). No significant differences were observed in the PBMCs' SNX25 expression of affected and healthy twins ( Figure 6B). However, when SNX25 expression was tested in PBMCs of an independent cohort of 7 patients with Pso/PsA and 3 healthy controls the same trend observed in RNAseq data was present ( Figure 6C).

SMARCA4/BRG1 and SMAD3 Show Methylation Differences with Discordant β-Value Variation
A Gene Ontology (GO) analysis of DMPs associated to the 19 regions showed an enrichment in transcription factor binding, transcription corepressor and transcription coactivator activity, SMAD binding and histone -lysine-N-methyltransferase activity (p < 0.005) ( Figure 7A).
Among the 19 genes identified with the Infinium Methylation EPIC array, we found significant interactions for SMARCA4/BRG1 (transcription activator BRG1, brahma-related gene-1) and SMAD3 (mothers against decapentaplegic homolog 3), by performing proteinprotein interaction analysis (String v11.2 database; www.string-db. org) ( Figure 7B). The BRG1/SMAD3 interaction was indeed already demonstrated in the literature: immunoprecipitation from cells overexpressing FLAG epitope-tagged BRG1 and HA-tagged SMAD proteins showed that BRG1 is a direct interactor of SMAD3 (Xi et al., 2008). Based on this observation, we further analyzed the methylation status of these two genes. Methylation scatterplots showed a significantly different methylation level of the probes (p < 0.005) localized in the SMAD3 and SMARCA4 genes, in affected and healthy  Figures S3A,B) although with discordant β-value variation in the CG probes. We interrogated the promoter of SMARCA4/BRG1 and SMAD3 to evaluate any change in the level of methylation in whole PBMCs from affected and healthy twins. Immunoprecipitation of 5′-methylcytosine antibody followed by quantitative PCR was performed with specific primers for SMARCA4/BRG1 and SMAD3, but we could not identify significant differences between affected and healthy twins (n 4 twin couples, Supplementary Figure S3C). Further, no changes were found in SMARCA4/BRG1 and SMAD3 mRNA expression in PBMCs (Supplementary Figure S3D). We also evaluated the mRNA expression of SMAD3 and SMARCA4/BRG1 in the selected control cohort of 7 patients with Pso/PsA and 3 healthy controls, confirming no significant changes (Supplementary Figures  S3E,F).

DISCUSSION
Epigenetic changes have widely been proposed to be crucial in the pathogenesis of complex diseases (Petronis, 2010) and, more specifically, DNA methylation is linked to different immunemediated disorders, from rheumatoid arthritis to multiple sclerosis (Liu et al., 2013;Dick et al., 2014). While epigenetic changes are ideal mechanisms to determine disease discordance in MZ twins, translating EWAS data into a quantifiable epigenetic difference and mechanistic understanding for complex diseases like Pso and PsA (herein cumulatively identified as psoriatic disease) is challenging, particularly in the absence of disease biomarkers. We performed the first genome-wide integrative approach to identify differential DNA methylation in whole blood from a selected cohort of MZ twins discordant for psoriatic disease, similar to what has been performed in other complex diseases (Javierre et al., 2010;Rakyan et al., 2011). Our EWAS approach identified distinct DNA methylation patterns at specific DMPs differentiating affected and healthy twins, mapping to the gene body of 1,703 genes while we did not observe genome-wide significant differences [p < 3.6 × 10 −8 (Saffari et al., 2018)]. The ENCODE and Blueprint regulatory annotations showed an enrichment for H3K36me3 histone mark overlapping the significant DMPs. DNA promoters methylation is largely known for its repressive effect on transcription initiation of protein-coding genes, non-coding RNAs, or transposon repeats (Teissandier and Bourc'his, 2017), while the function of intragenic methylation remains unclear (Baubec et al., 2015) and it may be involved in the modulation of alternative splicing, in exon definition and recognition (Maunakea et al., 2013). The commonly accepted dogma where methylation negatively correlates with gene expression as and vice versa, not always hold as, in different cancers, within promoter regions it has been observed a positive correlation between methylation and gene expression (Chatterjee and Vinson, 2012;Kim et al., 2016;Spainhour et al., 2019). H3K36me3 is required for preventing ). An asterisk (*) next to a p-value indicates the term also has a significant adjusted p-value (<0.05).
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 778677 8 spurious transcription initiation stemming from transposon fragments or cryptic promoters, preserving transcription process to canonical promoters (Verma et al., 2018). It has been recently demonstrated that H3K36me3 may act as a docking site for DNMT3B (DNA Methyltransferase 3 Beta), and together with reduced H3K4me3, altogether they repress the default activity of intragenic CpG islands and act as transcriptional promoters (Baubec et al., 2015). Further, the enrichment of methylated cytosines in regulatory regions that we identified in blood cell subsets suggests a dysregulation in psoriatic disease, which might be cell-type specific.
The enrichment analysis performed on transcriptome data suggested a significant role of "oxidative phosphorylation," "inflammatory response" and "myc targets" pathways in disease susceptibility. Indeed, these pathways have been already associated to autoimmune disorders, and also to Pso/ PsA. Regarding the "oxidative phosphorylation" pathway, it was observed that there is a direct link between dysregulated glucose metabolism in lymphocytes and autoimmunity (Bantug et al., 2018). In fact, enhanced glycolytic activity is a common feature of pro-inflammatory effector CD4 lymphocytes such as T helper (Th) 1 and Th17, cells characterized by a Warburg-like glycolytic metabolism, since they require aerobic glycolysis for their differentiation and effector functions. On the other hand, Regulatory T (Treg) cells, which suppress inflammatory responses and promote tolerance, after the initial activation, are characterized by an enhanced oxidative metabolism (Kornberg, 2020). The co-expression analysis again evidenced the dysregulation of this pathway, since we found an enrichment of the "oxidative respiration" pathway in cluster 14, characterized by a higher gene expression in non-affected MZ twins. It is worth mentioning that one of the top dysregulated genes in the transcriptome data is PDK4, coding for Pyruvate Dehydrogenase Kinase 4, a protein implicated in glucose metabolism. The gene was upregulated in affected MZ twins (logFC 0.64, unadjusted p 1.1*10 −4 ). A recent study evaluated the effect of PDK4-deficiency on the development of experimental autoimmune encephalomyelitis, the mouse model of multiple sclerosis, another autoimmune disorder, induced by pathogenic Th17 cells. Very interestingly, knock-out mice developed a less severe disease and showed a decrease in Th17 cells and an increased infiltration of Foxp3+ Tregs in the central nervous system (Allon Wagner et al., 2020). All together these results suggest that in affected MZ twins there might be a higher activation state of pro-inflammatory, and hence diseasesupporting, Th lymphocytes. Also MYC is one of the key players that coordinate metabolic reprogramming and activity in immune cells (Gnanaprakasam and Wang, 2017), which play a pivotal role in the development of inflammation and autoimmunity.
The finding that seven genomic regions are characterized by a concordant β-value variation in all the probes distributed, either with increased or decreased methylation level between affected and healthy twins, is of primary importance in order to define a specific methylation signature associated with psoriatic disease. Among these seven genes, SNX25, a negative regulator of TGFβ pathway and intracellular trafficking, was previously found associated among the differentially methylated sites between psoriatic and normal skin (Verma et al., 2018). Interestingly, transcriptome data pointed to a deregulation of this gene in affected MZ twins. Further analysis performed to replicate the results showed a trend, although not significant, both in methylation status and in expression level of the gene. Additional analyses are needed on a larger cohort of patients to confirm these results.
A small number of the 19 genes with at least two DMPs within 1 kb of distance and significant within-pair Δβ-values was found significantly associated with psoriatic disease. Among those, SMAD3, which is involved in the downstream signaling pathway of TGF-β, is associated with BRG1 to mediate TGFβ-induced transcriptional regulation at multiple genes, by means of chromatin remodeling and gene expression regulation. Although the methylation level of the CpG probes located in SMAD3 and SMARCA4/BRG1 was significant, this did not follow a consistent trend. Further, despite no significant changes in RNA expression were observed, a possible contribution of the differential methylation pattern in the definition of the alternative splicing profile of the gene cannot be ruled out.
A few studies evaluated the role of DNA methylation in the disease susceptibility specifically in psoriasis (Gervin et al., 2012;Li et al., 2020) at genome-wide level in immune cells. The most comprehensive one (Gervin et al., 2012), focused on lymphocytes, didn't disclose any differentially methylated regions between cotwins; however, when gene expression was also considered, different genes were identified. Among the 50 top hits, 7 of them (14%) were found differentially methylated in our dataset. Moreover, GO analysis revealed enrichment in processes associated with "immune response" and "psoriasis." This is in accordance with our results, since in the co-expression analysis we highlighted an enrichment in genes associated with "inflammatory response," "IL6/ STAT/JAK signaling," as well as involved in psoriasis. The role of inflammatory response was also evidenced by another genome-wide transcriptome analysis (Rawat et al., 2020) where a consistent, though small, pattern of changes was observed for a set of genes associated with "neutrophil activation" and "inflammation." Even though all the data point to the involvement of the same pathways in disease susceptibility, the lack of a high degree of overlap in the identified genes can be due to the use of different starting samples.
We are aware of the strengths and limitations of our study. Among the former, the chosen model of MZ twins is a unique setting to investigate complex diseases, as hereditary factors are not the major determinants of immunological functions, as proven in healthy MZ twins (Brodin et al., 2015). Further, using twins can improve the statistical power of a genetic study as it reduces the amount of genetic and/or environmental variability. Conversely, the results obtained from twin studies cannot be universally correlated to the general population, due to lack of randomization (Sahu and Prasuna, 2016). We are aware that our twin's cohort is heterogenous (see Table 1 for patients' characteristics), as it includes patients affected by psoriasis, psoriatic arthritis and both manifestations simultaneously. Further, the multitier approach with a small independent cohort of unrelated patients and controls also represents an advantage for our study. Among the study limitations, DNA methylation was measured in peripheral blood cells while epigenetic changes associated with psoriatic disease may be present in specific tissues such as the skin, but the collection of these samples through biopsies is invasive, and the collection of synovial fluid may not be feasible in most of PsA cases when synovitis is not the predominant domain. For this reason, peripheral blood is the best accessible alternative that reflects FIGURE 7 | Pathway enrichment and protein-protein interaction analysis of psoriatic disease MZ twins. (A) Top 25 differentially enriched GO terms identified by the "topGO" Bioconductor tool, using the "molecular functions" database, ordered by elimFisher p value. The width of the dots indicates the percentage of the enriched genes out of the total number of genes belonging to each term. (B) protein-protein interaction analysis evaluated with stringdb (https://string-db.org).
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 778677 different pathophysiological pathways and there is high potential clinical utility for any identified blood-derived epigenetic disease marker. Nevertheless, we are aware that the results we obtained in whole blood, could represent a limitation, as the observed changes cannot be associated with a specific cell type (i.e., T-cells, monocytes).
In conclusion, our genome-wide study on the unique model of MZ twins demonstrates high similarity in whole blood-based methylomes of psoriatic disease-discordant samples. However, we identified DMPs specific for psoriatic disease and few candidate loci that warrant to be further evaluated. The evaluation of the causal relation between genetic variants, biomarker levels and DNA methylation is becoming of undisputed importance. Their interplay is relevant to define plausible biological pathways in psoriatic disease that are regulated by epigenetic mechanisms.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GSE186713 for methylation data, GSE186724 for RNA-seq data.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Study-specific ethical approval was provided by the local ethics committee. All study participants provided written informed consent. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
All authors contributed to manuscript preparation. AC, EP, NI, FM, MV, MD, and CS actively participated to rheumatological evaluations and experiments, and supervised the whole project from the initial planning to the present manuscript; EP, GC, MR, RA, and SD performed and supervised the epigenetics and transcriptomics experimental parts; MS and SB enrolled twin couples involved in the study; FS and AC evaluated the enrolled patients from the dermatological point of view; MV, EP, AC, SD, and CS drafted the manuscript, and all the authors revised the final version prior submission.

FUNDING
AC received funding from Fondazione Eli Lilly per la ricerca medica ONLUS (Grant number ICH 2205). MV received funding from Versus Arthritis UK (Grant number 21428).