Identication of Crucial Methylation Genes Associated With Immune Cell Inltration in Twin-to-twin Transfusion Syndrome Using Co-methylation Analysis

Background: DNA methylation (DNAm), is an important transcriptional regulation mechanism, relevant to various diseases. Twin-to-twin transfusion syndrome (TTTS) is a complication in twin pregnancies resulting from disproportionate blood circulation. Survivors of TTTS show a high risk of neurodevelopmental abnormalities, particularly in the hippocampus, which is important in learning and memory. Here, we investigate gene expression and DNAm in hippocampus tissues of TTTS specimens. Methods: DNAm and gene expression levels were compared among the three groups: 10 recipients, 10 donors, and 10 matched control, using methylation microarray. We further explored the immune inltration of six immune cell sub-populations using EpiDISH analysis. The methylated sites related to immune cell inltration were identied using the WGCNA package. We explored the core methylation genes in the protein-protein interaction network using the MCODE plugin in Cytoscape software. Results: There were 188 differential methylation sites among three groups. Based on WGCNA, we found that the turquoise module containing 174 CpG sites is signicantly related to the immune inltration level. And four hub genes correlated with immune inltration level, namely, PTPRJ, FYN, LYN, and AKT1, and were identied using gene sub-network analysis. Conclusions: We identify the four hub methylation genes related to immune inltration in the TTTS. The molecular function of hub genes is still explored in the future research.

Background TTTS occurs in 10 ~ 15% of monochorionic (MC) twin pregnancies during the second trimester [1,2]. MC twins are monozygotic and result from the cleavage of a single zygote, resulting in two fetuses sharing the same genotype [3]. TTTS is thought to occur secondary to an imbalance in the transfer of blood between two monochorionic twins through placental vascular anastomoses [4]. The fetus that outputs blood is the donor, while the fetus that receives blood is the recipient. The donor develops features of anemia and hypovolemia, including oliguria and growth restriction, while the recipient shows signs of hypervolemia with polyuria, polyhydramnios, and visceromegaly [5]. These pathophysiological changes may in uence the growth and development of fetuses, especially neurodevelopment. Long-term neurodevelopmental disorders (NDIs) are present in 30% TTTs survivors, severe brain damage in 8-18%, and cerebral palsy in 3-12% [6]. Therefore, studying the mechanisms of TTTS on fetal neurodevelopment will be helpful for the diagnosis and treatment of fetuses with TTTS.
Learning and memory are important nervous system functions and important indices of neurodevelopment. In mammals, the hippocampus is the core region related to learning and memory function, and the number of nerve cells and neurogenesis in the hippocampus have been proven to be closely related to the learning and memory function [7]. DNA methylation (DNAm) is the most widely investigated epigenetic modi cation. DNAm is the mechanism whereby a methyl group is covalently added to cytosine residues in a Cytosine-phosphate-Guanine context (CpG) [8]. DNAm is essential for numerous biological functions such as ensuring chromosomal stability and genomic imprinting [9]. The precise regulation of DNAm is essential for normal cognitive function. The pathogenesis of many complex neurological and psychiatric diseases involve DNAm dysregulation, including autism and schizophrenia [10].
In this study, we analyzed three groups of hippocampus tissues using methylation array, including TTTS recipients, TTTS donors, and matched controls. We screened hub methylation genes related to immune in ltration based on the co-methylation analysis in TTTS.

Recruitment and sample collection
The donors (N = 10) and recipients from TTTS pregnancies (N = 10) and matched controls (N = 10) were recruited from the Shengjing Hospital of China Medical University from June 2018 to June 2019. TTTS diagnosis and staging were determined with ultrasound criteria according to Quintero [11,12]. TTTS cases were selected from cases of inevitable abortion before 24 weeks. In the control group, all were singletons who were spontaneously aborted before 24 weeks. There was no signi cant difference in the time to termination of pregnancy between TTTS and matched controls (P >0.05). The gravidae approved the use of hippocampus tissue for scienti c purposes by signing informed consent forms. The study was approved by the Human Ethics Committee of the Shengjing Hospital of China Medical University (2018PS360K).

Tissue processing and DNA extraction
After the termination of pregnancy, the brain tissues were immediately removed, and the hippocampus tissues were dissected on ice and washed using ice-cold phosphate-buffered saline (PBS). Then, the hippocampus tissues were snap-frozen in liquid nitrogen for 10 min and stored at −80 ℃ until processing. Genomic DNA was extracted from hippocampal tissues using the AllPrep DNA/RNA MiniKit (80204, Qiagen, Hilden, Germany) according to the manufacturer's protocol.

Illumina HumanMethylation850K bead chip
Genomic DNA from each sample was prepared for sodium bisul te conversion using the EZ DNA methylation Gold Kit (Zymo Research, Irvine, CA, USA). We further assessed genome-wide DNAm based on Illumina In nium HumanMethylation850K BeadChip (Illumina Inc, San Diego, CA, USA) microarray.
The array data were analyzed using ChAMP package in R for deriving the methylation level [13]. We have denoted the methylation status of all probes as the β value, which is the ratio of the methylated probe intensity to the overall probe intensity.

Differential methylated CpG sites analysis
The Limma package was used to identify differential methylated CpG sites [14]. We used |Δβ| ≥0.20 and adjusted P value ≤0.05 as differentially methylated site cut-off points. The CpG sites with Δβ ≥0.20 were considered hypermethylated, and the CpG sites with Δβ ≤−0.20 were considered hypomethylated.
Function enrichment analysis WEB-based Gene Set Analysis Toolkit (WebGestalt) is a functional enrichment analysis web tool, which has on average 26,000 unique users from 144 countries and territories per year according to Google Analytics. A total of three well-established and complementary methods for enrichment analysis in the WebGestalt: Over-Representation Analysis (ORA), Gene Set Enrichment Analysis (GSEA), and Network Topology-based Analysis (NTA) [15]. We performed the gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the identi ed differential methylation probes (DMP) based on the WebGestalt. P value < 0.05 was used as statistical signi cance.

Epigenetic Dissection of Intra-Sample Heterogeneity
The EpiDISH (Epigenetic Dissection of Intra-Sample Heterogeneity) package provides tools to infer the fractions of a priori known cell subtypes present in samples. The EpiDISH package provides a function called CellDMC, which allows the identi cation of differentially methylated cell types in epigenome-wide association studies (EWAS). The EpiDISH package contains four references: two whole blood subtype references, one generic epithelial reference of total immune cells, and one breast tissue reference. [16-18].

Co-methylation network analysis
To analyze the intercorrelation among the identi ed DMPs, we performed the co-methylation network analysis using the WGCNA package [19]. WGCNA is a network analysis tool that uses hierarchical clustering of correlated methylation states between DMPs to construct weighted co-methylation modules. WGCNA can transform the adjacency matrix into a topological overlap matrix. DMPs can be divided into different co-methylation modules based on the topological overlap matrix dissimilarity measurements. Here, we set the soft-threshold power as 16 to identify key co-methylation modules. The modules with the highest correlation with immune in ltration were selected for further analysis.

Protein-protein interaction network analysis
We constructed the protein-protein interaction (PPI) network using the Search Tool for the Retrieval of Interacting Genes (STRING) database. We upload the differential methylation genes to the STRING database and de ned the cut-off as the interaction score of 0.4 (median con dence) [20]. The PPI network was visualized using Cytoscape software [21]. Molecular Complex Detection (MCODE) analysis was used to screen the key gene modules within the PPI network. A k score ≥3 was selected as the cutoff to de ne the key gene module. We further explored the expression level of the hub differential methylation genes among the three groups.

Differential methylated CpG sites analysis
We select the probes with Δβ ≥0.2 and adjust P value ≤0.05 for further analysis. We found 1676 CpG site DMPs between the recipient group and the donor group. Of these, 1204 probes were hypermethylated in the recipient group and 472 probes were hypomethylated (Fig.1A). A total of 120,334 probes were differentially methylated between the recipient group and the matched control group. Among these, 776 probes were hypermethylated in the recipient group and 119,558 were hypomethylated (Fig.1B). Compared to that in the matched control group, there were 143,807 gene probes differentially methylated in the donor group; 316 probes were hypermethylated, whereas 143,491 probes were hypomethylated (Fig.1C). A total of 188 probes were differentially methylated in three groups (Fig.1D).

GO and KEGG pathway functional enrichment analysis
We performed GO functional annotation and KEGG pathway enrichment analyses on the nearest neighbor genes of differentially methylated CpG sites among the three groups. The results of GO annotation revealed that changes in biological processes are signi cantly enriched in spinal cord development, regulation of nervous system development, positive regulation of neurogenesis, positive regulation of neuron differentiation, and positive regulation of neuron projection development ( Fig.2A).
Changes in molecular function are mainly enriched in glutamate receptor binding, cytoskeletal protein binding, ionotropic glutamate receptor binding, T cell receptor binding, and Ras guanyl-nucleotide exchange factor activity ( Fig.2A). Changes in cell components are mainly enriched in the cell junction, cell projection part, plasma membrane-bound cell projection part, cell-cell junction, and synapse ( Fig.2A). The results of KEGG pathway analysis show that pathways are mainly enriched in the Fc epsilon RI signaling pathway, adherens junction, Fc gamma R-mediated phagocytosis, human papillomavirus infection, and platelet activation (Fig.2B).

Immune in ltration by EpiDISH analysis
We rst explored the difference in immune in ltration in hippocampus tissues of the recipient group, the donor group, and the matched control group from TTTS in six subpopulations of immune cells based on EpiDISH analysis [16-18]. We visualized the results obtained from 10 recipients, 10 donors, and 10 matched control hippocampus tissues (Fig.3A&B) [22]. There were no differences in the proportion of B cells, NK cells, CD4T cells, and CD8T cells. Compared with the matched control and donor tissues, the hippocampal tissues from the receipt group were found to contain a higher proportion of monocytes (Fig.3C).

Results of co-methylation network analysis
To further investigate CpG sites related to immune in ltration, we performed co-methylation analysis using the WGCNA package in the R environment (Fig.4A&B) [19]. Based on WGCNA analysis, highly interconnected CpG sites were clustered in the same modules. In each module, CpG sites were co-methylated. We identi ed a total of three modules. We generated the heatmap of module-trait correlations by setting the soft-threshold power as 16. We identi ed that the turquoise module was most highly correlated with monocyte in ltration (Fig.4C) with an especially high correlation coe cient of 0.44 (P = 0.02). The pink module contained a total of 174 CpG sites as shown in Fig.4D. We performed further analysis for the turquoise module CpG sites.

Gene network analysis
We uploaded the genes related to the CpG sites in the turquoise module to the STRING database and constructed the PPI network with 40 nodes and 47 edges using Cytoscape software. Among the 84 nodes, 40 central node genes were identi ed with medium con dence (0.4) (Fig.5A). To determine the core sub-network, we set the k score as 3 using the MCODE plug-in. The core module contains 4 genes, namely, PTPRJ, FYN, LYN, and AKT1 (Fig.5B).

Discussion
Monochorionic diamniotic (MCDA) twins comprise 20-30% of spontaneous and 4-5% of iatrogenic twin gestations. MCDA pregnancies are at higher risk for perinatal complications than dichorionic twins and singleton pregnancies [23]. This is attributed to the presence of vascular anastomoses connecting the fetal circulations. About 10-15% of MCDA twins may develop TTTS [1]. TTTS is a heterogeneous clinical entity and may display a wide spectrum of severities. It is di cult to completely understand the reasons for the range of clinical variations observed [24]. Research on long-term neurodevelopmental outcomes found that the survivors of TTTS often have neurodevelopmental injury. Both recipients and donors are equally at risk for adverse neurodevelopmental outcomes [25]. In this study, we investigated the hub methylation genes using methylation microarray to explore the pathogenesis of neurological impairment in TTTS.
As research has shown that changes in the immune microenvironment play an important role in neurodevelopment impairment [16-18], we further explored the difference in immune in ltration among three groups in six sub-populations of immune cells based on EpiDISH analysis. We found no difference in the proportions of B cells, NK cells, CD4T cells, and CD8T cells. However, the proportion of monocytes in the recipient group is higher than the donor and matched control groups. Gene sub-network analysis found the four hub genes: PTPRJ, FYN, LYN, and AKT1. PTPRJ, also called DEP-1 or CD148 is a widely expressed receptor-like protein tyrosine phosphatase comprising an extracellular domain consisting of eight bronectin-type repeats, a transmembrane domain, and a cytoplasmic phosphatase domain [26,27]. In the human immune system, PTPRJ is widely expressed on mature thymocytes, peripheral T cells, B cells, NK cells, granulocytes, macrophages, and DC populations [28]. Schneble et al. reported that the protein-tyrosine phosphatase PTPRJ is a positive regulator of microglial cell migration and phagocytosis, and PTPRJ may exert its effects through negative regulation of the SRC-family kinase Fyn [29]. Microglia that can mature into distinct populations of the central nervous system monocytes are considered the most susceptible sensors of brain pathology [30].
FYN, belonging to the Src family kinases, is a non-receptor or cytoplasmic tyrosine kinase that is expressed in the neurons and glia in the nervous system [31]. FYN is primarily in several transduction pathways in the central nervous system, such as axon guidance, myelination, and oligodendrocytes formation. FYN regulates mechanisms related to learning and memory processes and is attributed to the development of several neurological disorders including Alzheimer's disease and multiple sclerosis [32]. Grant et al. found that fyn mutant mice showed impaired spatial learning, which further illustrates that FYN may involve the core process of the growth of neurons in the developing hippocampus [33].
LYN is also a member of the Src family. Its activation relies on the dephosphorylation of the C-terminal inhibitory tyrosine by membrane-bound receptor-type phosphatases [34]. Using single-cell RNA-seq, Roberts et al. found that LYN is a key regulator of monocytes homeostasis, signaling, and gene expression [35]. LYN may be involved in the neuroprotective effect during ischemic injury in the hippocampal CA1 neurons. Zhang et al. reported that positive modulation of alpha-amino-3-hydroxy-5methyl-4-isoxazolepropionic acid receptors (AMPARs) could exert neuroprotective effects via Lyn-ERK1/2-CREB signaling [36]. AKT1, which is highly expressed in the brain, is involved in memory formation and synaptic plasticity. AKT1 contributes to several cellular functions such as cell survival, cell growth, and metabolism [36]. AKT1 is signi cantly linked to schizophrenia, a neurological disorder characterized by brain dysconnectivity, with a molecular basis in aberrant synaptic plasticity [37]. As one of the effectors in the pathway from neuregulin 1 (the NRG1-ERBB4-PI3K-AKT1 pathway), lower AKT1 protein levels are expressed in schizophrenia [38]. Hippocampus formation is prominently implicated in the pathogenesis of schizophrenia [39].
In this study, we focused on immune microenvironment-related genes involved in the pathogenesis of hippocampal impairment in TTTS using a methylation microarray. Collectively, we put forward the hypothesis that hyperperfusion/hypoperfusion may change methylation levels of the core genes PTPRJ, FYN, LYN, and AKT1 and further lead to changes of immune in ltration level that are related to neurological impairment. Declarations YSW and LCX conceived and designed the study. LC analyzed the data and wrote the paper. YSW revised the paper for important intellectual content. All authors read and approved the nal manuscript.

Availability of data and materials
The data used to support the ndings of this study are available from the corresponding author upon request.
Ethics approval and consent to participate All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee, and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Informed consent was obtained from all individual participants included in the study.

Consent for publication
Not Applicable.