Epigenetic consequences of hormonal interactions between opposite‐sex twin fetuses

Abstract Previous studies reported inconsistent evidence about some phenotypic traits of females in human opposite‐sex twins (opposite‐sex females [OSF]) being distinct from females in same‐sex twins (SSF). Comparatively, less evidence showed significant differences between males in OS twins (opposite‐sex males [OSM]) and males in same‐sex twins (SSM). The twin testosterone transfer hypothesis suggests that prenatal exposure of testosterone in utero may be a possible explanation for the differential traits in OSF; however, the underlying mechanism is unknown. Here, we investigated the potential epigenetic effects of hormone interactions and their correlation to the observed phenotypic traits. In the study, DNA methylomic data from 54 newborn twins and histone modification data (H3K4me3, H3K4me1, H3K27me3, and H3K27ac) from 14 newborn twins, including same‐sex females (SSF), OS twins, and same‐sex males (SSM) were generated. We found that OSF were clearly distinguishable from SSF by DNA methylome, while OSM were distinguishable from SSM by H3K4me1 and H3K4me3. To be more specific, compared to SSF, OSF showed a stronger correlation to males (OSM and SSM) in genome‐wide DNA methylation. Further, the DNA methylomic differences between OSF and SSF were linked to the process involving cognitive functions and nervous system regulation. The differential H3K4me3 between OSM and SSM was linked to immune responses. These findings provide epigenetic evidence for the twin testosterone transfer hypothesis and offer novel insights on how prenatal hormone exposure in utero may be linked to the reported differential traits of OS twins.

findings provide epigenetic evidence for the twin testosterone transfer hypothesis and offer novel insights on how prenatal hormone exposure in utero may be linked to the reported differential traits of OS twins.

K E Y W O R D S
DNA methylation, histone modification, opposite-sex twins, twin testosterone transfer hypothesis

INTRODUCTION
The incidence of spontaneous dizygotic twinning occurs between 1% and 4%, which might be increased by advanced maternal age and the application of assisted reproductive technologies. 1 Among the total dizygotic twins, approximately 40% are opposite-sex (OS). 2 Androgens, particularly testosterone, serve important functions in early embryonic development. 3 Testosterone production in fetuses increases from 8 to 24 weeks of gestation and reaches maximal level between 10 and 15 weeks. 4,5 It overlaps with the period of rapid brain development when microglial cells colonize the cerebrum during 4-24 gestational weeks. 6 Twin testosterone transfer hypothesis states that a female having a male co-twin is exposed to higher levels of prenatal testosterone. 7 Potential differences in physiological, cognitive, and behavioral traits have been reported between females in OS twins (OSF) and females in SS twins (SSF). 5 Compared to SSF, OSF may show more masculine traits, such as conservatism, 8 breaking rules, 9 poorer performance in high school and college, less likely to marry, reduced fertility, and lower income. 10 Some studies also suggested no differences between OSF and SSF. [11][12][13] In contrast to OSF, males in OS twins (OSM) show weak and inconsistent evidence of behavioral and biological changes. 10,14,15 Animal studies showed that testosterone triggered a DNA demethylation event in mouse embryonic neural stem cells and could affect the global acetylation pattern of histone H3. 16 Cisternas et al. found that inhibiting DNA methylation in neonatal mice disrupted the testosterone-dependent masculinization of neurochemical phenotypes. 17 In addition, human germline cells undergo global DNA demethylation from 7 to 19 weeks, coincident with testosterone production. 18 Therefore, epigenetic modifications may be a candidate mechanism to explain the more masculine traits in OSF.
To explore the potential influences of prenatal hormone exposure on fetuses, we generated epigenomic data from OS and SS newborn twins. We excluded effects of postnatal environment and socialization so that the epigenetic consequences of prenatal hormones interactions could be investigated exclusively. By comparing DNA methylation and histone modification data, we aim to address how epigenetics may be altered by hormone exposure and may lead to differentiated traits between OSF and SSF, as well as OSM and SSM.

Study population and sample collection
A total of 56 families having twins were recruited and then divided into four main categories: OSF, OSM, SSF, and SSM. The twins were born at the Third Hospital of Peking University; an accurate initial medical history confirmed that none suffered birth defects. On the delivery day, ∼1 mL of umbilical cord blood was collected from each twin and processed immediately. Whole-blood aliquots of 300 µL were used for DNA extraction, and the remaining blood was stored at −80 • C for backup. We also collected ∼20 mL blood and extracted cord blood mononuclear cells (CBMCs) as soon as possible.

DNA extraction
The QIAmp DNA Blood Mini Kit (Qiagen, cat.: 51106) was used to extract genomic DNA (gDNA) from the stored whole blood samples.

CBMC extraction and fixation
CBMCs were isolated by Ficoll density-gradient centrifugation (TBD, cat.: LDS1075). Freshly-prepared 1% formaldehyde in 1 × phosphatebuffered saline (PBS) was used to fix CMBCs for 10 minutes at room temperature. Glycine buffer (125 mmol/L) was added to quench the crosslinking reaction, and the solution was incubated for another 5 minutes at room temperature. The fixed CBMC pellet was cryopreserved at −80 • C after washing once with cold 1 × PBS for subsequent chromatin immunoprecipitation sequencing (ChIP-seq).
Bisulfite conversion was conducted with a MethylCode Bisulfite Conversion Kit (Thermo Scientific, cat.: MECOV-50). The converted samples were size-selected by excising gel slices containing 160-700 bp DNA (2% TAE gel). A gel DNA recovery kit (Vistech, cat.: PC0313) was used to recover and then amplify the converted DNA fragments through PCR with a maximum of 12 cycles with Kapa HiFi U+ Master Mix (Kapa Biosystems, cat.: KK2801). Barcodes were introduced during this process.
A Fragment Analyzer Automated CE System (Analysis Kit: cat.: DNF-474-0500) was used to assess the fragment distributions, and a Library Quant Kit for Illumina (NEB, cat.: E7630L) was used to measure molar concentrations. The final qualified libraries were sequenced using the PE150 strategy on an Illumina X Ten sequencer.

ChIP-seq
ChIP-seq was performed as previously described. 20

RRBS data analysis
Before read mapping, adapters and low-quality bases in the reads were removed with Trim Galore to optimize pairedend alignments. The kept reads were aligned to the Homo sapiens reference genome (human GRCh38/hg38) using Bismark (version 0.18.1) with the default parameters. The uniquely mapped reads with <2% mismatch were retained for further analyses. Differentially methylated cytosines (DMCs) between any two groups were identified using methylKit with a q-value threshold of 0.001 and a mean methylation difference of 20%.
In the RRBS samples, the 17.65% (3/17) of SSM and 33.33% (6/18) of SSF are monozygotic twins. In this study, we only focused on the epigenomic characteristics of OS twins (OSF and OSM) by comparing them with SS twins (SSF and SSM). Prenatal hormone exposure of twins is related to the genders of twins, but not zygosity. Thus, from the perspective of hormonal interactions, there are no differences between monozygotic SSF and dizygotic SSF, as well as monozygotic SSM and dizygotic SSM.

ChIP-seq data analysis
To optimize paired-end alignment, Trimmomatic was used for trimming before read mapping. 21 The kept reads were aligned to the Homo sapiens reference genome (human GRCh38/hg38) using the Burrows-Wheeler Aligner with the default parameters in paired-end mode. 22 The nonredundant reads with ≤2% mismatches were retained for further analyses. The MACS2 peak caller 23 (version 2.1.0) with the parameter setting "-keep-dup = 1, -broad" was used to identify peaks of histone modification. The R/Bioconductor package DiffBind was used, 24 which uses DESeq2 or edgeR to identify differential peak (DPs) with q-values <.05 and fold changes >2 for H3K4me1, H3K4me3, H3K27ac, and H3K27me3. The shared outputs of the DESeq2 and edgeR methods were taken as the final DPs.

Comparisons in different groups
To study the characteristics of OSF, we compared the Pearson's correlation coefficients of OSF with males (OSM or SSM

Clustering, functional enrichment analysis, and statistics
Hierarchical clustering was used for unsupervised clustering of the epigenomes of OS and SS samples with the R function "hclust," in the R base package stats. Kmeans clustering was used for unsupervised clustering of DMCs using the R function kcca with the parameter "family = kccaFamily (which = NULL, dist = distCor)" in the flexclust package. To calculate distance, we used the dist function, and for similarity matrices we used the cor function of the stats package in R with the parameter "method = correlation." Potential biological functions of genomic regions, such as DMRs and DPs, were determined by GREAT with a p-value threshold of 10 −5 . For violin boxplots, the center dot indicated the average; the significance of differences between two groups was determined by the Wilcoxon rank sum test.

OS twins exhibited distinct DNA methylomic characteristics
To systematically study the epigenetic differences between OS twins and SS twins, we performed RRBS on cord blood collected from 18 SSF,17 SSM, and 19 OS twins (Tables S1 and S2). We conducted hierarchical clustering on OSF versus SSF and OSM versus SSM based on the DNA methylation level of top 3% CpG sites with the highest standard deviation. The results showed that samples in OSF versus SSF or OSM versus SSM were obviously clustered into corresponding subgroups ( Figure 1A,B), revealing the different DNA methylation features of OS twins and SS twins with same genders.
To further identify the distinction between OS twins and SS twins, we performed inter-group correlation analysis for the four different groups (OSF, OSM, SSF, and SSM), and the comparisons in siblings were excluded. The correlation between OSF and males (OSM or SSM) was significantly higher than the correlation between SSF and males. By contrast, there were no significant differences between the correlation coefficient of "OSM versus females (OSF or SSF)" and that of "SSM versus females" (Figures 1C  and 2). Meanwhile, OSF was more strongly correlated with OSM than with SSM ( Figure 1C; Table S3). Furthermore, the intra-group correlation coefficient of "OSF versus OSF" (twin pairs were excluded) was highest among females versus females while "OSM versus OSM" did not show this characteristic (Figure 2). This indicated that the DNA methylome of OSF became more homogeneous among individuals, potentially shaped by common processes occurring in utero.

Different DNA methylation patterns are related to their potential distinct functions in OSF and OSM
In order to study the potential biological functions of hormone-induced DNA methylomic changes, we identified DMCs with q-values <0.001 and methylation differences >20% in "OSF versus SSF," and "OSM versus SSM." The results showed that there were ∼1000 DMCs in these two comparisons. Between the two comparisons, only 10% (108) of hypo methylated cytosines (hypo-DMCs) and 10% (108) of hyper methylated cytosines (hyper-DMCs) were overlapped ( Figure 3A and B; Table S4). To further determine whether OSF and OSM induced same methylation changes, we performed cluster analysis based on DNA methylation levels of merged hyper-DMCs and hypo-DMCs in all samples. The data showed that for both hyper-DMCs and hypo-DMCs, no consistent patterns of DNA methylation levels were present in any cluster of OSF and OSM ( Figure 4A). In addition, the overlapping proportion of DMCs associated genes was less than 20% ( Figure 3C); the genomic distances between DMCs in OSF and OSM were mostly more than 100 kb ( Figure 3D). These results suggested that the induced changes in DNA To identify the potential impacts of those specific alterations in OSF and OSM, we performed gene ontology enrichment analysis and found that hyper-DMCs of "OSF versus SSF" were significantly enriched in nervous system and catabolic process, such as neuropeptide catabolic process, hormone catabolic process, and collagen catabolic process. Comparatively, hypo-DMCs of "OSF versus SSF" were also mainly involved in nervous sys-tem regulation, such as negative regulation of neuron projection regeneration, negative regulation of axon regeneration, and regulation of integral component of postsynaptic membrane (Figures 3E and 4B,C; Tables S5 and  S6). Therefore, the differential DNA methylation in OSF might mainly affect regulation and catabolic processes of the nervous system. Comparatively, DNA methylation changes in OSM affected RNA catabolic process and neuromuscular process controlling posture ( Figure 3F,G; Tables S7 and S8).

3.3
Differential histone modifications are more observed in "OSM versus SSM" group than "OSF versus SSF" group In addition to DNA methylome, histone modifications may also be altered by the environment in utero. 25,26 Therefore, we used ChIP-seq in core blood mononuclear cells collected from seven OS twins, four SSF twins, and three SSM twins to examine changes in histone modifications, including H3K4me1, H3K4me3, H3K27ac, and H3K27me3. Based on genome-wide histone modifications, neither the "OSF versus SSF" group nor the "OSM versus SSM" group could be clearly identified and separated in hierarchical clustering analysis (Figure 5A-D; Tables S9 and S10). Differential peak (DP) analysis revealed that more than 1000 GAIN H3K4me1 DPs and GAIN F I G U R E 4 Differentially methylated cytosines (DMCs) for OSF versus SSF, and OSM versus SSM. A, K-means clustering for the all merged hyper-and hypo-DMCs. The hyper-DMCs or hypo-DMCs for OSF versus SSF, and OSM versus SSM were merged if they had 1-bp common region at least, respectively. The hyper-DMRs (hypo-DMRs) were classified into six clusters by k-means clustering on row scaled DNA methylation level. B and C, Enrichment terms of gene ontology (GO, biological process) with adjusted p-value <10 −5 were shown for the hyper-DMCs (B) and hypo-DMCs (C) of OSF versus SSF. The terms associated with nervous system were highlighted by red.
H3K4me3 DPs were identified in "OSM versus SSM," while only a few GAIN/LOSS DPs could be identified in "OSF versus SSF" (Figure 5E). Through GAIN DPs of H3K4me1, OSM, and SSM could be clearly distinguished ( Figure 5F). Furthermore, functional enrichment analysis showed that GAIN DPs of H3K4me1 between "OSM versus SSM" were enriched in the processes of gland morphogenesis, DNA replication, and transport ( Figure 5G; Table S11). In addition to H3K4me1, OSM and SSM could also be separated by GAIN DPs of H3K4me3 ( Figure 6A), which were highly enriched in immune response ( Figure 6B; Table  S11). We also focused on further analyzing the relationship between DMCs and DPs. Almost no overlaps in genomic regions between DMCs and DPs were observed, along with only a few DMC-and DP-associated genes (including GAIN H3K4me1 and GAIN H3K4me3) overlapped among all groups ( Figure 7A,B). The genomic distances between DMCs and DPs (including GAIN H3K4me1 and GAIN H3K4me3) were mostly more than than 100 kb ( Figure 7C).

DISCUSSION
Prenatal hormone interactions in OS twins may explain their different phenotypic traits in some clinical cohorts, while the underlying mechanism is still unclear. In this study, we recruited 56 pairs of twins (OS and SS twins) to exclude potential confounders such as postnatal environmental and social factors. For the first time, the epigenetic consequences of prenatal hormone interactions were investigated. We systematically explored the differences in genome-wide DNA methylation changes and histone modifications between OS and SS twins. From DNA methylome data, OSF presented different characteristics compared to SSF. "OSF versus males (OSM or SSM)" showed stronger correlation than "SSF versus males," suggesting that OSF tended to be more masculine and might become closer to their OSM brothers in DNA methylome. In contrast, OSM and SSM did not have significant differences when comparing their correlations with OSF and SSF. Therefore, OSM may be less influenced by their female co-twins according to the DNA methylation data.
We further investigated the potential functions of DNA methylome changes in OSF. Analysis of DMCs showed methylation changes were associated with nervous system development and regulation, which might explain the reported traits about differential perceptual and cognitive in OSF. In two cohort studies, OSF between 20 and 38 months tended to be somewhat weak in expressive vocabulary and scored lower on the MacArthur Communicative Developmental Inventory than SSF. 27,28 OSF also performed better than SSF twins in visuospatial cognition assessed by the Mental Rotation Test. 29,30 Although DMC-related genes were only correlated with nervous system functions in OSF, both OSF and OSM presented a pronounced enrichment of DMC-related genes involved in catabolic processes. A Swedish cohort study of twins (>60 years old) reported that OSF had a moderately higher body mass index, body weight, and rate of dyslipidemia than SSF. 31 The phenotypes and associated potential risks need further attention.
In addition to DNA methylation, epigenetic histone modification may regulate chromatin status and gene expression. 32 Previous studies were mainly focused on the epigenetics of monozygotic twins to study the cause for some types of phenotypic discordance between monozygotic twin pairs. 33 Our study offers some insights on histone modification differences between OSM and SSM. We found that DPs of H3K4me1 and H3K4me3 could be used to distinguish OSM from SSM; yet none of the four histone modifications could distinguish OSF from SSF. Interestingly, we found that differential H3K4me3 modifications in "OSM versus SSM" might be related to immune responses. A Danish cohort study showed lower early-life mortality risks in OSM than SSM. 34 When evaluating the effects of epigenetic modifications in OS twins, age may be a critical factor. [35][36][37][38] An Australian study using two cohorts (younger cohort and older cohort) reported that masculine characteristics (OSF scoring higher on rule-breaking than SSF) only occurred in the younger cohort (mean = 23.2), while the older cohort (mean = 41.2) showed no significant differences. 9 Another study from Denmark and Sweden showed no significant differences in the incidence rate ratios of cancers between OS and SS twins. 39 Previous studies were mainly focused on adult twins and may therefore obtain inconsistent results. Our study used newborns, which excluded postnatal environmental and age as potential confounders, so that the hormone effects of OS twins in utero could be investigated exclusively. However, our study also has some limitations. The sample size was not very large compared to clinical cohort studies. We would need more clinical samples to verify our results. In the same-sex twins, we didn't differentiate their zygosity (including monozygotic twins and dizygotic twins). It may be an influence factor when we performed inter-group comparing analysis between SS twins and OS twins. The specific biases and differences introduced by monozygotic or dizygotic twins need further investigation in a larger cohort considering the small portion of monozygotic twins in our study. The postnatal follow-up data should also be analyzed in the future to further confirm the epigenetic changes and to investigate the effects of age. However, for those traits that may only appear long after, the phenotypes could also be strongly impacted by environmental factors which may also be a confounding factor.
In conclusion, we provide an epigenetic basis to the understanding of hormonal interactions in OS twins. Our results highlighted that DNA methylation of OSF was significantly influenced by their male co-twins in newborns with potential impacts on nervous system development and regulation. DMCs involved in catabolic processes are prone to be affected in both OSF and OSM, and immune responses in OSM might be affected through histone modifications. Long-term follow-up is necessary to provide more solid and insightful evidence.

C O N F L I C T O F I N T E R E S T
The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

E T H I C A L A P P R O VA L
All blood samples were obtained after participants have given written informed consent and were fully anonymized (each individual was given a code). The Reproductive Study Ethics Committee of Peking University Third Hospital approved the study (approved protocol no. 201752-044 in 2017/06/20). All relevant ethical regulations were followed.

A U T H O R C O N T R I B U T I O N S
Siming Kong, Yong Peng, Wei Chen, and Liying Yan wrote the manuscript. Wei Chen and Xinyi Ma collected the study materials, samples, and patient data, and performed the experiments. Yong Peng developed the analysis methods and performed bioinformatics analysis. Jie Qiao, Liying Yan, and Wei Chen developed the experimental concept and design. All the authors read and approved the final manuscript.

C O D E AVA I L A B I L I T Y
The bioinformatics pipelines and scripts used for our analysis are available at https://github.com/CTLife/Opposite-Sex-Twins_Epigenome.

D ATA AVA I L A B I L I T Y S TAT E M E N T
All the raw Data included in this study have been uploaded to the Sequence Read Archive database (https://trace.ncbi. nlm.nih.gov/Traces/home/) and are available for download via accession number GSE136849.