Genomic map of candidate human imprint control regions: the imprintome

ABSTRACT Imprinted genes – critical for growth, metabolism, and neuronal function – are expressed from one parental allele. Parent-of-origin-dependent CpG methylation regulates this expression at imprint control regions (ICRs). Since ICRs are established before tissue specification, these methylation marks are similar across cell types. Thus, they are attractive for investigating the developmental origins of adult diseases using accessible tissues, but remain unknown. We determined genome-wide candidate ICRs in humans by performing whole-genome bisulphite sequencing (WGBS) of DNA derived from the three germ layers and from gametes. We identified 1,488 hemi-methylated candidate ICRs, including 19 of 25 previously characterized ICRs (https://humanicr.org/). Gamete methylation approached 0% or 100% in 332 ICRs (178 paternally and 154 maternally methylated), supporting parent-of-origin-specific methylation, and 65% were in well-described CTCF-binding or DNaseI hypersensitive regions. This draft of the human imprintome will allow for the systematic determination of the role of early-acquired imprinting dysregulation in the pathogenesis of human diseases and developmental and behavioural disorders.


Introduction
Data from model systems and humans demonstrate that environmentally induced epigenetic modifications occurring early in development can cause longterm gene expression changes in important mechanistic pathways involved in disease pathogenesis. Affected diseases include neurological disorders [1,2], cardio-and cerebrovascular diseases [3], and cancers [4][5][6], and their major risk factors, such as obesity and the attendant dysfunction in metabolism, nutrient acquisition, fat deposition, appetite, and satiety [7,8].
Covalent DNA methylation of cytosines in CpG dinucleotide sites is the most studied epigenetic modification and is hypothesized to link environmental exposures to these diseases. Nevertheless, data generated from case-control or crosssectional studies, the most cost-efficient epidemiologic study designs, have been difficult to interpret. This is because methylation marks measured in easily accessible peripheral cell types from otherwise healthy individuals do not always reflect methylation in inaccessible cell types, tissues, and organs involved in the formation or progression of a chronic disease. Moreover, DNA methylation at susceptible loci can adaptively change throughout life in response to environmental exposures or disease. Indeed, methylation levels can also naturally diverge with normal cell differentiation, ageing, and environmental influences [9][10][11][12].
Known exceptions are CpG methylation marks that are stochastically established before specification that control metastable epiallele expression [13] and imprinting control regions (ICRs) that regulate the monoallelic expression of imprinted genes [6,14]. CpG methylation of metastable epialleles and ICRs is established before gastrulation and is mitotically heritable. Thus, these epigenetic marks are normally similar across tissues and cell types throughout an individual's life. Unlike metastable epialleles however, ICRs are defined by parent-of-origin specific methylation marks that are important gene dosage regulators based on the allele's parental origin. Consequently, in contrast to epigenetic marks controlling metastable epiallele expression, methylation marks regulating imprinted genes are similar across individuals [15,16]. Importantly, changes in the methylation patterns in ICRs are implicated in adult-onset diseases suspected to have foetal origins, including neurological disorders, cancers, and metabolic diseases stemming from abnormal growth and nutrient acquisition disorders [17,18].
Together, these features make ICRs attractive targets for dissecting disease aetiology, particularly since imprinted genes comprise an estimated 1-6% of the human genome, are over-selected for growth regulators, and are critical in early embryonic development [3,19]. The stability of these methylation marks with age also makes them longterm 'records' of early exposures that are difficult to obtain through questionnaires or other exposure assessment assays [6]. Yet, despite their biological and clinical relevance, only 25 of the ICRs regulating the 100-150 identified human imprinted genes, and of the 300-1,000 genes predicted to be imprinted [19], are currently known [14].
Leveraging existing specimens, recent advances in genome sequencing, and computational capabilities, we sought to comprehensively characterize human ICRs using DNA methylsequencing of tissues representing the three germ layers as well as the gametes, to create a catalogue of the set of imprint regulatory DNA methylation marks in humans, the 'human imprintome' [14]. This characterization of human ICRs should enable a more detailed understanding of the epigenetic basis of numerous pathophysiologies with foetal origins, and with that, advance the ability to diagnose, prevent, and treat a number of developmental disorders and diseases.

Materials and subjects
Sperm from three adults and tissues from twelve 65 to 95-day-old embryos of both sexes (confirmed by sex-linked marker genotyping) and of African and European descent were selected for bisulphite DNA sequencing. These tissues were obtained from the National Institutes of Health funded Laboratory of Human Embryology at the University of Washington, Seattle, WA; they were snap frozen to preserve DNA/RNA integrity (NCSU Institutional Review Board #3565). Embryonic tissues were used for identifying ICRs because the gametic and somatic imprint marks are intact, and because monoallelic gene expression of imprinted genes occurs primarily during embryonic development [20][21][22][23]. Sperm DNA was from the TIEGER study at Duke [23].

Whole-genome bisulphite sequencing (WGBS)
Libraries for NextSeq sequencing using HiSeq 2500 were prepared from bisulphite converted DNA derived from tissues representing the three germ layers using previously described methods [21]. Thirty (three from sperm obtained as previously described [23,24] and 27 from somatic tissues) of the 36 samples passed quality control standards for sequencing by Illumina NextSeq with 12-15X coverage. The sequenced somatic tissue libraries were 8 kidney (mesoderm), 8 liver (endoderm), and 11 brain (ectoderm). Libraries were index-tagged for separating reads after multiplex sequencing and pooled into groups of nine, with each group split for sequencing into three separate lanes. Splitting samples across lanes ensured that no single sample was disproportionally affected by technical variability specific to an individual sequencing lane (e.g., low read numbers or low read quality). If a problem persisted, samples exhibiting consistently low quality across lanes were rerun or removed from analysis. These sequence data were supplemented with publicly available oocyte sequence data (JGAS00000000006) [25] and part of control sperm data from PRJNA754049 (control males in this study were otherwise healthy individuals who did not use cannabis) [26]. Details are in the Supplementary Materials and Methods.

Bioinformatic approaches to identify ICRs
Samples were separated by index sequences and aligned to a reference in silico bisulphite-converted genome. Reads without unique alignment to the reference sequence due to either repetitive sequence or loss of information because of cytosine conversion were eliminated, as were duplicate reads, indicative of clonally amplified original random DNA fragments. From these reads, methylation fractions and read counts were calculated for all CpG sites in the genome. We developed a candidate ICR identification pipeline (putICR) application using a Ruffus framework in Python. The workflow is described in detail in Supplementary Materials and Methods and outlined in Supplementary Figure S7.
This application was designed to scan the genome and identify regions of allelic differential methylation based on four criteria that define genomic imprinting: 1) ≥5 consecutive CpG sites, consistent with known cis-acting imprinted control regions [14]; 2) methylation levels of 50% ±15%, supportive of monoallelic methylation (i.e., approximately 100% methylation on one parental allele and 0% on the other); 3) similarity of methylation levels in tissues from the three germ layers (i.e., brain, liver, and kidney), as expected for methylation marks established before tissue specification; and 4) similarity of methylation across individuals, indicative of sequence regions playing critical roles in regulating imprinted gene dosage, which should not generally vary by sex, ethnicity, developmental age, or from person-to-person. Based on the methylation levels, we developed an online tool, the imprintome browser, which is linked to standard genome browsers to visualize the methylation level for each CG site in the 1,488 candidate ICR regions (https://humanicr.org/).
To assess the reproducibility of the putICR pipeline, the raw sequence data and the same four criteria were provided to an independent bioinformatics group (Sciome, Inc, Research Triangle Park, NC) (Supplementary Materials and Methods). Additionally, fully methylated or unmethylated regions from gametic sequences were compared, as these are the original inherited parent-of-origin specific regions.
We used pyrosequencing (sequencing primers in Supplementary Table S3) to determine if methylation patterns at previously identified ICRs in the three germ layers in embryos were similar to those obtained from DNA from accessible adult tissues frequently used in epidemiologic studies, such as peripheral blood components and human umbilical vein endothelial cells (HUVECs) (Figure 3).

WGBS library preparation
WGBS libraries were prepared from three human tissues -brain (ectoderm), kidney (mesoderm), and liver (endoderm) -representing the three germ layers from 12 embryos (6 male, 6 female), and three sperm samples, resulting in 39 total libraries. Twenty-seven of the tissue libraries (i.e., 8 kidney (mesoderm), 8 liver (endoderm), and 11 brain (ectoderm)), and the three sperm libraries passed quality checks for Illumina HiSeq2500 (Illumina, San Diego, California, USA), resulting in 30 total libraries sequenced (27 somatic and 3 gametic). Gametic sequence data were augmented with publicly available human oocyte sequence data (accession number JGAS00000000006 [25]) and control sperm sequence data from PRJNA754049 [26]. The average number of per-sample reads was 153 million (range 74-231 million), covering an average of 23.1 billion bases per sample (range 11.2-34.9 billion). Approximately 80% of reads uniquely aligned to the in silico bisulphite-converted human genome (hg38). Of the 29.2 million CpG sites in the human genome [27], an average of 26.6 million (91%, range 86-94%) were covered by aligned reads for the set of 30 samples.

Characteristics of candidate ICRs from tissues representing the three germ layers
Most of the sequences obtained from the brain, kidney, and liver (75%) had a sequence coverage of greater than 20X. The sequence coverage for oocytes (accession number JGAS00000000006) was lower due to the decreased availability of DNA from this source [25] while the sperm coverage, by incorporating data Figure 1. Detection of ICRs genome-wide. A bioinformatics selection algorithm, putICR, was used to identify methylation fractions using whole-genome bisulphite sequencing (WGBS) of DNA from kidney, liver, and brain embryonic tissues and gametes (i.e., sperm and oocyte). (a) Genome coverage was >20X for the somatic tissues in more than 75% of the fraction base sampled. The genome coverage was lower in the available oocyte gametic sequence data (accession number JGAS00000000006), but higher in sperm because our data were supplemented with control sperm data from PRJNA754049 [25,26]. (b) The size range of candidate imprint control regions (ICRs) averaged 248 bp. (c) This is similar to previously identified ICRs (Avg: 375 bp) [14].
from Schrott et al. [26], was higher (Figure 1 (a)). Using the percent methylation fractions calculated for each CpG site, candidate ICRs were defined based on the following criteria: five or more consecutive CpG sites within a 300 bp region having methylation levels of approximately 50 ± 10-20% in tissues from all three germ layers. These criteria are consistent with an ICR being a genomic series of cis-acting CpG sites that are established in embryonic stem cells, resulting in one parental allele being near fully methylated (~100%) while the other is unmethylated (~0%) for ≥80% of the sites.
Using the most relaxed criteria (50 ± 20%), we identified 7,559 candidate ICRs, including 21 of the 25 known ICRs [14]. A more stringent methylation criterion of 50 ± 15% decreased the number of candidate ICRs to 1,488 while still detecting 19 of 25 known ICRs (Table 1  and Supplementary Table S1). Further, restricting the window to 50% ±10% decreased the number of candidate ICRs to 127, including 15 of the known ICRs.
The 1,488 novel ICRs methylated at 50% ±15% ranged from 10 to ~4,000 bp long, with a median length of 248 bp (Figure 1(b)). This is similar to the median size of the known ICRs (375 bp) (Figure 1(c)), but with a tailing distribution of much longer candidates. As expected for ICRs, differentiated tissues derived from the three embryonic germ layers exhibited similar methylation levels, consistent with the establishment of these methylation marks in the stem cells before tissue specification. Importantly, these regions also showed similar methylation fractions across individuals, consistent with ICRs controlling gene dosage. We have developed a corresponding imprintome browser depicting methylation fractions for each CpG site in the 1,488 candidate ICRs identified for each embryonic germ layer (https://humanicr.org/).
As a sensitivity analysis, Sciome, Inc. (Research Triangle Park, NC, USA) used an independent calling pipeline with similar procedures for adapter trimming and alignment to hg38 and the four ICR call criteria, but without pooling aligned reads. This generated 1,225 ICRs, including 19 of the 25 characterized ICRs. Of the 1,225 ICRs discovered by Sciome, Inc. (Research Triangle Park, NC, USA), 900 (62%) were found in our initial analysis (Supplementary  Table  S1 and Supplementary Table S2).

Parent-of-origin methylation patterns in gametes
Parent-of-origin methylation patterns are a property of authentic inherited ICRs and can be discerned from the methylation patterns of gametes. We examined the parent-of-origin methylation patterns of the novel ICRs using WGBS results from sperm and oocyte sequence data. Oocyte data was obtained from public databases (accession number JGAS00000000006) [25]. Our sperm sequence data was supplemented with partial control sperm samples from PRJNA754049 [26]. We compared the embryonic and gametic methylation patterns of the 1,488 candidate ICRs seeking to identify those either fully unmethylated (<10%) or fully methylated (>90%) in gametes. This methylation pattern indicates control of genomic imprinting, which once established in the gametes and inherited persists throughout development. As proof of concept, we first examined whether this strategy could successfully identify known ICRs, and 19 of the 25 known ICRs [14] were captured using this strategy (Table 1 and  Supplementary Table S1). Despite the low readdepth of the oocyte sequence data, we found that 332 (Table 1) of the 1,488 candidate ICRs (Supplementary Table 1) had methylation patterns consistent with originating in either the sperm or egg (visualize in https://humanicr.org/).
The ICRs for most known imprinted human genes (https://www.geneimprint.com/site/genesby-species) are unknown, but most are now  [14]. They were also identified in our whole-genome bisulphite sequencing (WGBS) screen for candidate ICRs (red boxes); visualize at https://humanicr.org/. The ICRs defined herein overlap and extend beyond the currently known ICRs. (c) PTCHD3 (ICR_643) and (d) MCTS2P/HM13 (ICR_1191) contain only candidate ICRs (red boxes) determined in this study, and potentially control novel imprinted genes. The methylation levels for ICRs is 50 ± 15% for the candidate ICRs. DNA methylation levels determined by WGBS in sperm and oocytes are also provided. These sequence data were supplemented with publicly available gametic sequence (accession number JGAS00000000006) and part of the control sperm data from PRJNA754049 [25,26]. Dots indicate hemi-methylation (green), hypomethylation (blue), and hypermethylation (yellow) based on WGBS. The vertical dashed red lines delineate the candidate ICR regions.         Figure S2B). This is similar to what has been observed in the dog [32], but different from that observed in the mouse, which also has an ICR in the promoter region [33]. We also identified novel candidate ICRs mapping near PTCHD3 (ICR_643, Figure 2(c), Table 1) and MCTS2P/HM13 (ICR_1191, (Figure 2(d)), (Table 1). Other novel candidate ICRs include those mapping near WNT10A (ICR_177, Supplementary Figure S3A, Table 1) and ADNP2/PARD6G (ICR_1052, 1053, Supplementary Figure S3B, Table 1). Such findings indicate the presence of additional imprinted genes in the human genome.

Similarity of methylation marks across accessible tissues
For diagnostics and public health screening, it is critically important that DNA methylation marks be reproducible using other sequencing technologies and be similar to that in otherwise healthy human's accessible tissues (e.g., peripheral blood, maternal or foetal tissues discarded at birth (e.g., placenta), or HUVECs). These tissues also often serve as controls in epidemiologic studies, and act as surrogates for affected inaccessible tissues.
To assess reproducibility, we selected two wellcharacterized ICRs regulating the imprinted expression of PEG3/ZIM2 (ICR_1142, (Figure 3 (a)), (Table 1) and Supplementary Table S1) and PEG10 (ICR_475, Figure 3(b), Table 1 and  Supplementary Table S1) for which we were also able to developed pyrosequencing assays using the PCR primer sets provided in Supplementary Table  S3 to determine DNA methylation levels. We used pyrosequencing to test whether CpG methylation averaged 50% in these canonical ICRs samples. The DNA methylation patterns determined with pyrosequencing of embryonic brain, kidney, and liver fell between 35% and 65% and averaged approximately 50%. This was comparable to that found with WGBS in these well-characterized ICRs (Figure 3(c,d)).
Samples of internal tissues are generally unavailable for research, diagnostics or public health screening purposes. To determine if methylation marks at known ICRs are similar in accessible tissues from healthy individuals, we again characterized ICRs regulating the imprinted expression of PEG3/ZIM2 (ICR_1142, (Figure 3(a)), (Table 1) and Supplementary Table S1) and PEG10 (ICR_475, Figure  3(b), Table  1 and Supplementary Table S1) in CD14 − monocytes from newborn cord blood and HUVECs. The methylation patterns determined with pyrosequencing were comparable to those found with WGBS Regions of differential methylation meeting criteria reported herein for ICRs that overlap with regions of gamete-specific methylation. Parental allele methylated indicates which gamete source showed hypermethylation -Paternal (P -sperm) or Maternal (M -oocyte) based upon both sperm and oocyte methylation data. Previously reported ICRs that do not have gamete specific methylation are also included and labelled as having Somatic (S) methylation. ICRs overlapping ENCODE annotated regions of CTCF binding and DNase I hypersensitivity are indicated by * and ^, respectively. ICRs overlapping previously published ICRs of imprinted genes are indicated by #. + Criteria for an ICR: 1) ≥5 consecutive CpG sites, 2) methylation levels of 50% ±15%, 3) similarity of methylation levels in tissues from the three germ layers (i.e., brain, liver, and kidney), and 4) similarity of methylation across individuals.
in these well-characterized ICRs in the CD14 − monocytes and HUVECs (Figure 3(e,f)), falling between 35% and 65% methylation and averaging 50% in these accessible tissues as expected for bona fide ICRs. These results indicate that accessible peripheral tissues can be used to assess the effects with previously known (blue boxes) and candidate ICRs (red boxes); visualize at https://humanicr.org/. DNA methylation levels determined by WGBS in sperm and oocytes are also provided. These sequence data were supplemented with publicly available gametic sequence (accession number JGAS00000000006) and part of control sperm data from PRJNA754049 [25,26]. Dots indicate hemi-methylation (green), hypomethylation (blue) and hypermethylation (yellow) based on (WGBS on the human imprintome from early life exposures to chemicals, physical agents (e.g., radiation, blunt force trauma, heat stress, etc.), and other adverse physiological conditions.

Novel ICR methylation and human pathology
Many of the novel ICRs identified are in regions previously implicated in the pathogenesis of human diseases. For example, four candidate differentially methylated ICRs are located within the Down syndrome (DS) critical region at chromosome 21q22 [34]. ICR_1317 resides in the promoter of HLCS (Supplementary Figure S4A, Supplementary Table S1); ICR_1318 is in intron 1 of RIPPLY3 (Supplementary Figure S4A, Supplementary Table S1); ICR_1319 is close to the 3' UTR of KCNJ6 (Supplementary Figure  S4B, Supplementary Table S1); and ICR_1320 is in intron 1 of GET1 (Supplementary Figure S4C, Supplementary Table S1). The functions of these genes have developmental implications, playing roles in the metabolism required for infant growth and development (HLCS); transcriptional regulation controlling pharyngeal development (RIPPLY3); cell membrane potential regulating G-protein coupled receptors in both cardiac and neuronal signalling, with mutations connected to developmental delay with facial abnormalities and intellectual phenotypes (KCNJ6); and intracellular transport and positioning of proteins involved in signal transduction pathways connected to retinal deterioration and nystagmus (GET1).
Interestingly, ICR_1319 and the highly restricted Down syndrome critical region (HR-DSCR) of only 34 kb [35] flank KCNJ6 on distal 21q22.13, suggesting that not only gene duplication but also imprinting dysregulation may be involved in DS. Moreover, the first intron of the GET1 gene harbours a single nucleotide polymorphism (rs2244352, chr21:39,386,047) (Supplementary Figure S4C) implicated in comitant esotropia, a condition resulting in poor binocular vision that affects ~20% of Down syndrome cases. The SNP rs2244352 is associated with risk in a parent-of-origin dependent manner, with significant association between incidence and paternal inheritance of the minor 'T' allele [36,37]. The major allele for rs2244352 is G. Thus, carrying the minor T allele would eliminate the formation of a CG dinucleotide able to be methylated; however, loss of methylation on the paternal chromosome is likely not the causative mechanism.
It has been shown that this region is differentially methylated, but that it is the maternal allele that is methylated [36]. Our gametic data is also consistent with maternal methylation, showing hypomethylation in sperm and hypermethylation in oocytes (Supplementary Figure S4C). With no loss of paternal methylation possible to cause dysregulation of this putative ICR, we hypothesize an alternate mechanism, in which repressor binding occurs on unmethylated CG sites in this region. This would be comparable to regulation at the H19/Igf2 ICR, by CTCF binding of the unmethylated allele [14]. Thus, if 'CG' to 'TG' substitution alters the recognition sequence for repressor binding, paternal GET1 expression could be activated, adding to the expression from the maternal copy protected from repression by methylation of the putative ICR.
Similarly, four candidate ICRs are located within the DiGeorge syndrome critical region at chromosome location 22q11.2 (Supplementary Figure S5 A, C) [38,39]. Hemizygous microdeletions in this crucial region result in developmental disorders, such as velopharyngeal insufficiency, variable conotruncal heart defects, and cognitive and behavioural disorders (e.g., schizophrenia, bipolar disorder, and autism) [38]. Moreover, differential brain effects have been reported with the maternal or paternal deletion of 22q11.2, suggesting a role of imprinted genes in the aetiology of these psychiatric abnormalities [38]. The ICRs identified are ICR_1355 in the promoter of DGCR5 (Supplementary Figure  S5A, Supplementary Table S1); ICR_1356 is in the promoter region of TBX1 and ICR_1357 is intronic (Supplementary Figure S5B, Table 1 and  Supplementary Table S1). ICR_1358, which overlaps pseudogene ABHD17AP4 (Supplementary Figure S5C, Supplementary Table S1), is upstream of SERPIND1 and SNAP29 (Supplementary Table  S1). Functionally, DGCR5 is a long non-coding RNA that acts as a regulator of apoptosis and proliferation; TBX1 is a transcription factor involved in the regulation of developmental processes, with its deletion directly linked to physical malformations seen in DiGeorge syndrome; SERPIND1 is a protease inhibitor with critical functions in blood clotting; and SNAP29 is a member of a family of synaptic vesicle trafficking regulators associated with developmental disorders of the CNS.
We previously used a computer algorithm to predict the genome-wide imprint status of human genes from sequence features [19]. DLGAP2 was predicted to be imprinted and demonstrated to be paternally expressed. It is a membrane-associated protein that plays a role in synapse organization and signalling in neuronal cells. It was subsequently implicated in autism based upon CNV analyses [40]. In this study, we identified an ICR_502 in intron 1 (Supplementary Figure S6A, Supplementary Table S1); this gene also contains ICR_503-506 (Table 1 and  Supplementary Table S1). Thus, the function of DLGAP2 could potentially be altered both genetically and epigenetically in the formation of autism. Interestingly, of the 102 annotated genes we predicted with the use of computer algorithms to be imprinted [19], we have now identified ICRs for 35 of them (34%) (Supplementary Table S1 and  Supplementary Table S4), providing additional supporting evidence that they are imprinted.
PRDM16 functions as a transcription coregulator in the development of brown adipocytes and increased expression may protect against obesity [41]. We previously predicted it to be imprinted and paternally expressed in both mice [42] and humans [19]. Interestingly, ICR_11 and ICR_12 (Supplementary Table S1), while overlapping TTC34, are also within 500kb of PRDM16. Although the oocyte methylation data are sparse, it appears that the paternal allele at these ICR regions is sparsely or entirely unmethylated while the maternal allele is methylated, suggesting that PRDM16 is paternally expressed as predicted by Luedi et al. [19]. Using less stringent criteria (i.e., 50 + 20% DNA methylation), there is also an ICR in the promoter region of PRDM16 (ICR_3070_54) and three additional intronic ICRs (i.e., ICR_3070_55-57; Supplementary Figure S6B). Should this gene be experimentally confirmed to be imprinted, genomic imprinting would be involved in the maturation of both white [43] and brown fat cells, making the role of imprinted gene expression in metabolism and obesity more extensive than previously appreciated.
Imprinted genes have been identified on the X-chromosome in mice [44], but not in humans. Nevertheless, we have identified 98 putative ICRs on the X-chromosome (ICR_1391-1488, Supplementary Table S1), indicating that imprinted genes are also present on the human X chromosome. Gamete DNA methylation was adequate to determine parent-of-origin methylation for 11 of these ICRs (6 maternally methylated and 5 paternally methylated, Table 1). Only DHRSX (ICR_1394), a secretory protein associated with starvation-induced autophagy [45], is located in a pseudoautosomal region -PAR1 ( Figure 4). The role of these candidate imprinted genes in the genesis of parental-dependent behavioural disorders, such as those observed in Turner syndrome [46] needs further investigation.

Functional significance of ICR methylation
Bona fide ICRs typically regulate gene expression by controlling access to transcriptional regulatory sites within the imprinted genes. We examined the overlap of our 1,488 candidate ICRs with DNase I hypersensitive sites, regulatory motifs, and transcription factor-binding sites for nearby genes. Of the 332 ICRs that were either unmethylated or fully methylated in gametes, 200 (60%) overlapped with DNase I hypersensitive regions. Of those 332 ICRs, 178 were hypermethylated in sperm DNA sequences, and 154 were hypermethylated in oocytes.
Moreover, for each ICR, we looked for regulatory motifs within ±5,000 bp using Analysis of Motif Enrichment (AME) [47] against the Homo Sapiens Comprehensive Model Collection (HOCOMOCO). This consists of motifs for 680 known transcription factors. Approximately onethird of the 680 known transcription factorbinding sites are within or near candidate ICRs, supporting the functional significance of these ICRs (Supplementary Table S5).
The limited number of ICRs previously characterized [14] are over-selected for growth effectors, . Candidate ICR for DHRSX (ICR_1394) in the pseudo autosomal region, PAR1, on the X chromosome; visualize at https:// humanicr.org/. Whole genome bisulphite sequencing (WGBS) identified a candidate ICR (red box) in embryonic brain, kidney, and liver tissue for a candidate imprinted gene associated with starvation induced autophagy [45]. DNA methylation determined by WGBS in sperm and oocytes is also shown; these sequence data were supplemented with publicly available gametic sequence (accession number JGAS00000000006) and part of control sperm data from PRJNA754049 [25,26]. Dots indicate hemi-methylation (green), hypomethylation (blue) and hypermethylation (yellow) based on WGBS. The vertical dashed red lines delineate the candidate ICR regions.
with dysregulation associated with a wide range of conditions, including metabolic disorders, cancers, neurological diseases, language development deficits, schizophrenia, and bipolar affective disorders [48,49]. Thus, we examined the predicted functions of the 914 genes associated with the 332 candidate ICRs with gamete-specific methylation (Supplementary Table S6) by performing gene ontology (GO) analysis using Gorilla (http://cblgorilla.cs.technion.ac.il/) (Supplementary Table  S7). The most significant high-level biological processes identified by Gorilla were cAMP biosynthetic processes, cellular response to glucagon, and adenylate cyclase-activating G proteincoupled receptor signalling pathway (Supplementary Table S7). Functional analysis of candidate regions using the Comparative Toxicogenomics Database (CTD) identified protein-binding and DNA-binding as enriched GO terms in the functional annotations (Supplementary Table S8). Of the top five pathways identified by Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., https://www.qiagenbioinfor matics.com/products/ingenuity-pathway-analysis) , three were important neurological pathways: gonadotropin-releasing hormone (GnRH) signalling, endocannabinoid neuronal synapse pathway, and g-aminobutyric acid receptor signalling (Supplementary Table S9). Five hundred eight genes were associated with neurological disease, while 168 were related to nervous system development and function. Of the top networks identified by IPA, RNA post-transcriptional modification, cell cycle and DNA replication, recombination, and repair had the highest score, with auditory disease, cellular compromise, and neurological disease coming in second. Results of KEGG and REACTOME analysis of enriched pathways included neuronal system, transmission across chemical synapses, signal transduction, pathways in cancers, circadian entrainment, axon guidance, cholinergic synapse, glutamatergic synapse, and calcium signalling pathways (Supplementary  Table S10) [50].
We also explored the role of the candidate imprinted genes in human disease using the Mouse Genome Informatics (MGI) database and the CTD and their curation of OMIM disease categories and inferred diseases via chemical-gene chemical-disease associations [51]. Consistent with the preceding analyses, several neurological disorders correlated with a subset of the genes in MGI, including Alzheimer's disease (n = 2), autism spectrum disorder (n = 6), various ataxias, and cranial developmental disorders. When inferred disease associations were factored into the CTD analysis, the number of gene-disease associations restricted to neurological disorders expanded greatly, e.g., 242 genes were associated with Alzheimer's disease and 243 with Parkinson's disease. For Alzheimer's disease, 16  ; the remaining genes were shared by both disorders, indicating overlapping chemical mechanisms at the molecular level in the aetiology of these diseases.

Discussion
Epigenetic mechanisms, such as DNA methylation, are believed to link adverse intrauterine exposures to adult disease susceptibility; however, supporting empirical evidence from humans has remained scarce because of a lack of recognizable and archivable patterns of early epigenetic effects that can be detected and quantified in the epigenome [52][53][54]. This is especially important for human studies that rely on DNA from sample types accessible in otherwise healthy human populations (e.g., saliva or peripheral blood). While sequence regions controlling the monoallelic expression of imprinted genes have been previously proposed as targets for such studies [6], until now, only 25 of these regions have been described [14], with potentially hundreds more unknown.
Herein, we determined that the human imprintome is comprised of 1,488 regions with characteristics typically observed in known ICRs. Interestingly, the overlap between the 850,000 CpG sites in the Illumina Infinium Methylation EPIC microarray (Illumina, Inc., San Diego, CA, USA) and the 22,279 CpG sites in the human imprintome is only 7%. Using genomic DNA obtained from three tissue types isolated from embryos of both sexes and at least two ethnicities, these regions include most of the previously identified ICRs [14]. Of the 1,488 regions, 332 also exhibited parent-of-origin specific methylationthe hallmark of an ICR. Of those ICRs, 209 overlapped with DNase I hypersensitive regions. These novel ICRs have a median CpG dinucleotide content of 248, similar to the 25 previously characterized ICRs with a size range of 10 to ~4,000 bp. When an overlap was identified between a novel and a known ICR, the novel sequence typically extended beyond the boundaries of the prior reported sequence. Chromosomes frequently contain clusters of genes that are controlled by a single ICR (i.e., imprinted domains). Thus, our 200 strongest candidate ICRs, with overlapping DNase I hypersensitive sites and gamete-specific methylation, could regulate as many as 400 imprinted genes if each ICR controlled only two genes, which is commonly observed. Expanding this to include all 332 candidate ICRs with gametespecific methylation increases the number of potentially imprinted genes to a little more than 500 or ~3% of the genome.
Many known imprinted genes have developmental functions involving the regulation of cell function and growth such as pruning of synapses and adipocyte accrual affecting the life course. Thus, ICRs are of particular interest in studying the early origins of a wide range of common chronic diseases, including neurological disorders and metabolic diseases such as obesity, type 2 diabetes, and cancers. The majority (90%) of our candidate ICRs fall within 5,000 bp of genes involved in fundamental processes, including cAMP biosynthetic processes, cellular responses to glucagon and adenylate cyclase-activating G protein-coupled receptor signalling. Notably, a large proportion of the candidate ICRs are near genes involved in metabolic and neurological diseases, consistent with the observations of previously characterized ICRs. Characterizing and experimentally confirming the complete repertoire of ICRs could lead to the development of dietary interventions and pharmacological targeting of specific regions to advance the precision of nutritional and chemical therapies.
For example, if the candidate ICR proximal to the zinc finger transcription factor gene PRDM16 is experimentally confirmed to be imprinted, altered methylation can be detected in any tissue as early as birth. As PRDM16 controls the bidirectional fate decision between brown adipocytes and myoblasts [55], early detection of alteration in an imprinted ICR could provide opportunities for pharmacological or dietary manipulation to enhance the expression of PRDM16, preventing metabolic dysfunction, including obesity. Similarly, disposition to hypoplasia of the thymus and parathyroid glands and conotruncal cardiac malformations and schizophrenia that potentially result from imprinting disorders in the DiGeorge syndrome critical region could be identified earlier, and potentially even during gestation when interventions may be more effective.
ICRs have the unique features of the early establishment of DNA methylation marks, similarity across cell types and tissues, and stability over the lifespan. These characteristics facilitate their broad use as stable archives of early developmental exposures that may alter metabolic and other developmental and behavioural processes in adulthood [56][57][58]. Thus, ICRs are logical targets for evaluating the early origins of disease using accessible cell types obtained at variable ages [59][60][61][62]. A predetermined reference panel of ICRs could be vital in identifying early exposures. Their effects on DNA methylation have a wide range of potential uses, particularly in chronic disease epidemiology where relevant past exposures may be difficult to quantify. Such a predetermined reference panel could also be key in unravelling epigenetic responses to early life chemical and non-chemical stressors that result in molecular 'wear and tear' as reflected in methylation changes in accessible tissue, and may improve the precision and usefulness of epigenetic clocks [63,64]. The expectation that the baseline methylation fraction of ICRs is approximately 50% and is stable over time, with little epigenetic drift, also supports the use of ICR methylation as biomarkers of adult disease susceptibility at any time during the life course. Therefore, having the complete repertoire of bona fide ICRs -the imprintome -should improve our understanding of the early origins of adult diseases.
While the use of methylation marks as exposure proxies has long been advocated, and methylation patterns of ICRs provide a rare 'epigenetic responsive' window to early exposures, our results should be interpreted in the context of the study limitations. Firstly, our algorithm used to identify candidate genome-wide ICRs was tested at three methylation thresholds: relaxed 30-70% (50 ± 20%), moderate 35-65% (50 ± 15%), and stringent 45-55% (50 ± 5%). We selected the moderate set for further analysis. While this approach is pragmatic, it does not accommodate the possibility that methylation patterns defining ICR boundaries may be more fluid, requiring additional statistical and experimental interrogation of ICRs to refine the boundaries and, with that, the imprinted genes they regulate. Bioinformatic approaches such as change-point modelling that have been used in copy number variant analyses [65,66], coupled with experimental validation, could also be deployed in the future. Secondly, we have not performed clonal-allele analysis of candidate regions to definitively prove that they are bona fide ICRs, with methylation patterns consistent in cis for the parental alleles. Thus, we do not expect all candidates to overlap true ICRs; however, based on the characteristics of known ICRs, we are confident that a majority of human ICRs are captured within the 1,488 candidate ICRs.
Despite these limitations, we have used a combination of bioinformatic and sequencing approaches to provide the first draft of the complete repertoire of human ICRs -the human imprintome. To facilitate its use, we have also developed an online tool, the imprintome browser, which is linked to the UCSC Genome Browser to visualize the data at: https://humanicr.org/. Further refinement is needed to identify bona fide ICRs and their exact boundaries. This will be bioinformatically iterative and will require experimental validation. To our knowledge, however, this is the first study to create a human ICR compendium. As these sequence regions are also present in accessible peripheral blood DNA, we anticipate our data will greatly facilitate the ability to quantify the contribution of exposures in early development to a wide range of adult-onset chronic diseases; pave the way for novel earlydetection tools; and eventually reveal the molecular underpinnings of vulnerability in disease processes, especially in early life.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability
Data will be available from the authors upon request, and at the website https://humanicr.org/.