s2fieldsofstudy: Biology Chemistry

text:
Metabolic folate profiling as a function of time during cultivation suggests potential C2-metabolism in Saccharomyces cerevisiae

Yeasts are reported to be rich in folates, a group of vitamers known to be involved in several biosynthetic reactions such as methylation reactions, oxidation and reduction processes, and nucleotide synthesis. Not being able to synthesize folates, humans rely on external folate supply. Here, we show the application of LC/MS-MS methods using SIDA (stable isotope dilution analysis) assays for the quantitative analysis of different folate mono- and polyglutamates during growth of Saccharomyces cerevisiae. Molecular networking (MN) was applied for detailed analysis of further folate metabolites. Highest folate contents of 13,120 μg/100 g were observed after 20 h of cultivation. The main vitamers 5-CH3-H4folate and H4folate decreased during cultivation, while 5-CHO-H4folate increased during cultivation. The hexa- and heptaglutamate of 5-CH3-H4folate accounted for >96% of the total 5-CH3-H4folate content. A shift of the major polyglutamate from hexa- to heptaglutamate was observed after 29 h. MN unraveled two groups of novel folates which could be assigned to a potentially existing C2-metabolism in yeast. In detail, 5,10-ethenyl-tetrahydrofolate and a further CO-substituted 5-CH3-H4folate were identified as hexa- and heptaglutamates. The latter was neither identified as 5-acetyl-tetrahydrofolate nor as EthylFox, the oxidation product of 5-ethyl-tetrahydrofolate. The structure needs to be elucidated in future studies.

Introduction
Folates are a group of more than 150 structurally similar compounds, their entirety giving the vitamer class B 9 . The most stable form of these vitamers is the fully oxidized pteroyl glutamic acid (PteGlu), commonly known as folic acid. The common structural element of all compounds is given by the core structure consisting of a 4-amino-4-oxo-pteridine ring being linked to p-amino benzoic acid and one to fourteen γ-linked glutamic acid moieties (1)(2)(3). Different oxidation states within the pteridine ring, variation within the substitution at positions N 5 and N 10 (typically C 1 -units in the form of methyl, formyl, or methenyl groups), and varying length of the glutamyl tail result in the large number of folates known so far. Generally, folates serve as cofactors in the transfer of C 1 -units in a variety of biochemical reactions within the organism (4-6). Such reactions are for instance purine and thymidylate synthesis and thus DNA replication as well as methylation reactions (1,3). A part of this C 1metabolism is shown in Figure 1 in which the metabolic active 5-methyl-tetrahydrofolate (5-CH 3 -H 4 folate) serves as methyl donor for the reactivation of methionine by methylation of homocysteine (7). A lack of folates has been linked to several chronic diseases such as neurogenerative disease, cardiovascular diseases, or Alzheimer's disease (8-10). Most severe, however, is the incidence of neural tube defects in newborns as result of a maternal malnutrition with folates during pregnancy (11). Yet mammals are not able to synthesize folates de novo and thus rely on a sufficient external supply. Yeast has been proven to be a good folate source due to high folate contents above 4,000 µg/100 g.
To date, folate analysis mostly focused on quantitative analysis of different vitamers after enzymatic deconjugation into the respective monoglutamates (1,3,12). Yet little is known about the folate metabolism especially with regards to polyglutamates which has been shown to account for the majority of the folate pool in organisms (13)(14)(15). Consequently, there is a clear need for deeper research on the metabolism of folates as a whole.
In a previous study we could show that high resolution mass spectrometry (HR-MS) applying quadrupole time of flight MS (Q-ToF-MS) provides a more profound tool for folate metabolic profiling (16). However, folates are only rarely included into one of the numerous spectral databases available (a manual search in a selection of databases was applied). If at all, only spectral information for monoglutamates is accessible which makes a comparison of experimental data rather challenging. Therefore, other tools for data analysis need to be considered. A promising tool for this purpose is molecular networking (MN). MN offers the possibility for undirected analysis of fragmentation information of an experimental setup. Global Natural Products Social Molecular Networking (GNPS) offers a web-based application of MN which is freely available (17). GNPS uses the idea of clustering nearly identical tandem-MS spectra in one consensus spectrum (18,19). Those consensus spectra are matched against each other to find similarities within these spectra, assuming that structurally similar molecules show similar fragmentation behavior in MS/MS analysis. Spectra exceeding the user-defined similarity threshold will subsequently be visualized within one cluster "family." In those clusters, consensus MS/MS spectra are represented by nodes. Significant similarities above the chosen match score are visualized by connecting edges.
In this study we show the application of MN for a deeper analysis of folate metabolites. S. cerevisiae was chosen as model organism due to its reportedly high folate contents of >4,000 µg/100 g (20,21) and the possibility for cultivation according to our needs. The chosen yeast strain was cultivated in minimal medium with varying length of cultivation time to trace the development of the folate composition during yeast growth. Samples were quantitatively analyzed for their mono-and polyglutamate content by means of LC-MS/MS (liquid chromatography tandem mass spectrometry) to show the development of different monoglutamates and 5-CH 3 -H 4 folate polyglutamates during cultivation. Furthermore, samples were analyzed by means of Q-ToF-MS for further analysis of folate polyglutamates and subsequent application of MN. MN was applied to give deeper insight into the folate metabolism in yeasts and to uncover potential further folate vitamers with focus on the polyglutamates.

Chemicals and materials
Unless otherwise stated, solvents and chemicals were of analytical grade. Methanol (MeOH), acetonitrile (ACN), and formic acid were of LC-MS grade (Sigma-Aldrich, Germany). Water was derived from a Merck Millipore Integral 3 water purification system (Billerica, MA, USA). Folate standards and internal standards were purchased form Schircks Laboratories (Jona, Switzerland). Only MeFox was purchased from Merck & Cie (Schaffhausen, Switzerland). Chicken pancreas and rat serum were obtained from Biozol (Eching, Germany) and Difco (Sparks, MD, USA). Strata strong anion exchange (SAX) cartridges (quaternary amine, 500 mg, 3 mL, 1 g, 6 mL, and 2 g, 12 mL) were purchased from Phenomenex (Aschaffenburg, Germany); Chromabond hydrophilic lipophilic balanced (HLB) cartridges (60 µm, 500 mg, 3 mL) were obtained from Macherey-Nagel (Düren, Germany).  Munich. Yeast cultivation and the sampling procedure are described in detail in the Supplementary material. Baker's yeast for matrix matched calibration functions was bought in a local supermarket (Freising, Germany), freeze-dried, homogenized using a pestle and mortar, and stored at −20 • C until usage.
Sample extraction, analysis, and data preprocessing for folate profiling by UHPLC-Q-ToF-MS Sample preparation was performed under subdued light. Sample extraction of 450 mg of lyophilized yeast sample was as described in a previous publication for high resolution measurements by ultra-high pressure liquid chromatography (UHPLC) Q-ToF-MS (16). Samples were spiked with the internal standard (ISTD) [ 13 C 5 ]-5-CH 3 -H 4 folate right before injection into the UHPLC-Q-ToF-MS system. Parameters used for analysis can be found in the Supplementary material. Data preprocessing of UHPLC-Q-ToF-MS data was as in the above mentioned publication. Within the corresponding retention time window of folates, MS-features were run against an inhouse database to unravel putative folate vitamers. Identification was achieved by comparison of the respective tandem-MS data.

Quantification of the total folate content and of -CH -H PteGlu −
Detailed information about the folate stable isotope dilution assays applied for the quantification of the total folate content can be found elsewhere (16,22). The content of 5-CH 3 -H 4 PteGlu 1−7 was quantified by a modified method also using the internal standard (ISTD) [ 13 C 5 ]-5-CH 3 -H 4 folate as outlined in the Supplementary material.

Molecular networking
GNPS molecular networking parameters were set to a minimum requirement of five ions match, a minimum cluster size of 2 and a cosine score of >0.6 (https://gnps.ucsd.edu). Precursor ion and fragment ion mass tolerance were set to 0.02 Da. To further putatively identify folate metabolites in the yeast samples analyzed, mass spectral molecular networking was applied with a preliminary dataset analyzed by UHPLC-Q-ToF-MS. We then extracted the m/z of compounds clustering within the folate molecular cluster into measurements of the complete dataset by adapting the preference list for fragmentation accordingly (Supplementary Table 2). We then performed a second iteration of MN with the whole dataset. Different cultivation times were processed individually to receive independent networks for each time point. The obtained folate clusters were exported and further analyzed using RStudio software. Precursor masses of the folate clusters obtained by GNPS were compared with an inhouse database containing several folate vitamers. Annotation was verified with an annotation error of <10 ppm. Unknown features within the networks were manually reanalyzed in Compass DataAnalysis 4.3 (Bruker Daltonics, Bremen, Germany) and common features were searched.

Software and statistics
OriginPro 2020 (OriginLab Corporation, MA, USA) and Cytoscape 3.7.2 (cytoscape.org) were used for the generation of plots. Compass DataAnalysis 4.3 (Bruker Daltonics, Bremen, Germany) and RStudio software version 3.6.1 (Boston, MA, USA) were used for data processing. OriginPro 2020 was used for statistical evaluation.

Total folate contents
Total folate analysis by LC-MS/MS revealed decreasing folate contents in cultured yeast with increasing cultivation time as shown in Figure 2. Folate was not released into the culture medium as no significant amount of different folate vitamers could be observed in the culture media (data not shown). Pairwise t-tests were performed to evaluate significant differences between folate contents in dependence on the length of cultivation. Total folate contents were significantly different (p < 0.05) between the different timepoints tested except for 20 h of cultivation compared to 24 h and 24 h compared to 26 h. The 32 h timepoint did not show any significant differences compared to all the other timepoints tested expect for a cultivation time of 20 h. Highest folate contents could be observed after harvest at 20 h with a total sum of 13,119 ± 444 µg/100 g, whereas lowest contents were reached after harvest at 80 h with a total sum of 5,995 ± 253 µg/100 g. Main vitamers observable were 5-methyltetrahydrofolate (5-CH 3 -H 4 folate), 5-formyl-tetrahydrofolate (5-CHO-H 4 folate) representing the sum of all formyl vitamers, and tetrahydrofolate (H 4 folate) as depicted in Figure 2A. While 5-CH 3 -H 4 folate contents decreased from the first harvest timepoint on (79% of the total content at 20 h compared to 52% of the total content after 80 h), 5-CHO-H 4 folate showed clearly increasing tendency with 5% of the total sum after 20 h and 31% after 80 h. H 4 folate, however, increased from 13% (20 h) to 20% (26 h), before a decrease to 10% after 80 h could be determined. The oxidation products folic acid (PteGlu), MeFox, and 10-formyl-folic acid (10-CHO-PteGlu) (shown in Figure 2B) generally accounted for <2% of the total content highlighting the relative stability of the folate vitamers in yeast during cultivation. Solely 10-CHO-PteGlu reached > 4% after 80 h. The degradation product pABG (para-aminobenzoic acid) could not be quantified with the applied method. However, profiling measurements by UHPLC-Q-ToF-MS did not show any spectra with MS 1 features relating to pABG -neither as .
/fnut. . monoglutamate nor as polyglutamate with a varying length of the polyglutamate tail. Trends of folate contents of the oxidized forms were in accordance with tendencies observable for the corresponding reduced forms (PteGlu replicated the tendency of H 4 folate and 10-CHO-PteGlu replicated the tendency of 5-CHO-H 4 folate). Only MeFox showed a slightly different course during cultivation compared to 5-CH 3 -H 4 folate.

Profiling of folate metabolites during yeast cultivation by UHPLC-Q-ToF-MS
To investigate folate metabolites in cultured yeast in more depth, we analyzed the samples by means of ultra-high pressure liquid chromatography (UHPLC) Q-ToF-MS. Generated MS-features were filtered for the respective chromatographic behavior known by the injection of folate standards. Only features eluting within the retention time window of folates (±0.2 min) were kept as described in an earlier publication (16). The profiling of folate metabolites on different polyglutamate levels revealed the occurrence of glutamate moieties in the glutamate tail ranging from n = 2 to n = 8 as shown in Figure 4. In total, 45 different polyglutamates could be determined by their m/z (hereafter stated as MS 1 level, yellow nodes) of which 29 could further be identified by tandem MS/MS spectra (stated as MS 2 hereafter, black edging) for 24 h of cultivation. Twenty one folate polyglutamates were identified on MS 2 level for 80 h of cultivation. Decreasing intensities could be observed with increase of cultivation time (represented by the node size in Figure 4A). Furthermore, a shift in the main polyglutamate form from the hexaglutamates (after 24 h of cultivation) to the heptaglutamates (after 80 h of cultivation) was observed for all vitamers regardless their substitution (methylated, formylated and unsubstituted polyglutamates shown in Figure 4B). Supplementary Figure 11 illustrates comparable transformation rates throughout the cultivation process from the hexa-to the heptaglutamates for each of the vitamer groups. The oxidative degradation products (shown in Supplementary Figure 11B), however, were characterized by slightly higher hexaglutamate .
/fnut. . proportions in comparison to their corresponding metabolic active vitamers (shown in Supplementary Figure 11A).

Molecular networking and putative identification of unknown folates
To further putatively identify folate metabolites in the yeast samples analyzed, mass spectral molecular networking was applied (Supplementary Figure 12) for each individual cultivation time point. Compounds not eluting within the respective retention time of folates (known of the profiling part) ± 0.2 min were excluded from the networks. On average, more than 60% of the clustered features referred to already known folates for each of the cultivation time points. For identification of the remaining features, we focused on common features throughout different cultivation times. Of those, we selected the nodes being observable in more than two thirds of the obtained networks. Furthermore, we checked for metabolic activity of these nodes. Therefore, we plotted the corresponding intensity of the selected nodes as a function of time during yeast cultivation. Applying those filtering functions, we aimed at only selecting those features showing metabolic significance in the context of the already known folate metabolism. Doing so, we identified four promising features which could be grouped into two subgroups based on their fragmentation behavior (shown in Supplementary Figure 13A). Within both subgroups, the corresponding MS 1 -spectra revealed a mass difference accounting for an additional incorporation of one glutamate moiety between the features of the subgroup visualized in Figure 5A. Based on the number of neutral losses observable in each of the corresponding MS/MS-spectra, the features were putatively annotated as hexa-and heptaglutamates. The exemplary visualization of the stepwise neutral loss of glutamate moieties for one of the putative hexaglutamates is shown in Supplementary Figure 13B.
For means of identification, we searched for further mass differences between the unknown compounds and already known folate vitamers which are shown in Figure 5B. The hexaglutamate of the first subgroup (m/z 558.1984) was linked with 5,10-CH + -H 4 folate (m/z 551.1895) by a neutral loss of 7.0089 Da referring to CH 2 when accounting for the double charged state of the ions. Furthermore, it was linked with the hexaglutamate of the second subgroup (m/z 567.2024) by the addition of H 2 O (9.0050 Da). The latter revealed further mass differences accounting for the addition of CO (13.9972 Da) and CH 2 (7.0072 Da) to 5-CH 3 -H 4 folate (m/z 553.2060) and MeFox-Glu 6 (m/z 560.1952), respectively. As 5-CHO-H 4 folate, 10-CHO-H 4 folate and MeFox are isomeric folate forms, m/z 560.1952 could also represent the formyl vitamers.
The first group of metabolites was putatively identified as a C 2 variant of the C 1 -donor 5,10-CH + -H 4 folate, namely 5,10-ethenyl-tetrahydrofolate (5,10-CH + -CH 3 -H 4 folate) due to the additional incorporation of CH 2 into the 5,10-CH + -H 4 folate core structure. The existence of 5,10-CH + -CH 3 -H 4 folate recently has been described in Lactobacillus reuteri (23) and thus, its prevalence in other organisms such as yeast would be highly plausible. Compliance in the retention behavior as well as the fragmentation behavior of the C 1 -and the C 2 -metabolites observed in our studies, supported this identification further (shown in Supplementary Figure 14). Furthermore, the trend  of the two polyglutamates during yeast cultivation constantly increased (visualized in Supplementary Figure 15A). Thus, this gave further hint for their identification as formyl variant as the latter had been shown to gain in content towards the end of the yeast cultivation. However, due to the huge structural similarity of folate vitamers and generally low prevalence, a further purification could not be conducted. Consequently, the identification could not further be proven.
Several more mass differences to already known folates had been observed for the second group of metabolites enabling its annotation as C 21 H 25 N 7 O 7 for the monoglutamate form. Hence, different theories could be postulated for the structure .
/fnut. .  Figures  15A,B), theory 1 could be rejected. If they would represent different ions but refer to the same compound, they should show the same trends. Thus, according to theory (2) and (3), identification as a further C 2 -metabolite was more plausible.
In accordance with the potential existence of a C 2 -metabolism shown in Figure 6, the second group of metabolites could represent the precursor of 5,10-CH + -CH 3 -H 4 folate within the postulated C 2 -cylce, the acetyl vitamer form 10-acetyltetrahydrofolate (10-CO-CH 3 -H 4 folate). The latter would be in compliance with the additional incorporation of either a formyl group into the 5-CH 3 -H 4 folate core structure [theory (2)] or an additional methyl group into the 10-CHO-H 4 folate core structure [theory (3)]. As the core structure of theory (3) could also be the isomeric form MeFox, the second group of metabolites could also be the C 2 analog of the latter, which we name "EthylFox" throughout the rest of the manuscript. This compound would represent the oxidation product of a further C 2 -metabolite, namely 5-CH 2 -CH 3 -H 4 folate. To verify the identity of the second group of metabolites, two C 2 -metabolites, 5-CO-CH 3 -H 4 folate and EthylFox, were synthesized as monoglutamates. Information about the synthesis and verification are provided in the Supplementary material. However, comparative measurements of the synthesized standards by UHPLC-Q-ToF-MS unraveled different retention behavior of both standards. Thus, the group of metabolites could neither be identified as 5-CO-CH 3 -H 4 folate nor as the oxidation product of 5-CH 2 -CH 3 -H 4 folate, namely EthylFox. It has to be mentioned that the 5-CO-CH 3 -H 4 folate was synthesized instead of the metabolic active 10-CO-CH 3 -H 4 folate due to higher reactivity of the position N 5 . However, the C 1 -metabolites 5-and 10-CHO-H 4 folate are hardly separable by liquid chromatography and thus also for the C 2 -analog a similar elution was assumed. Admittedly, the possibility of identification of the second metabolite group 2 as 10-CO-CH 3 -H 4 folate remains even though this is highly implausible.

Discussion
Values for the total folate content decreased during cultivation of brewer's yeast from 13,200 µg/100 g (after 20 h) to 6,000 µg/100 g after 80 h of cultivation ( Figure 2). Those values are in agreement with previous findings reporting highest total folate content in cultured yeast of 12,500 µg/100 g (24). However, these authors observed highest folate contents after 10 h of cultivation, whereas we cultured for 20 h to reach this maximum. This could be due to different yeast strains used for experimental setup. Both studies, however, revealed 5-CH 3 -H 4 folate as main vitamer with about 80% of the total folate content at maximal folate contents.
Comparing the development of 5-CH 3 -H 4 folate and H 4 folate (Figure 2A) during yeast growth in both studies, .
/fnut. .  Of the oxidation products which we could determine with our methods, PteGlu and 10-CHO-PteGlu ( Figure 2B) showed the same tendencies compared to their metabolic active equivalents H 4 folate and 5-CHO-H 4 folate, respectively. Solely MeFox appears to rather follow the development of 5-CH 3 -H 4 PteGlu 7 ( Figure 3A) than the development of total 5-CH 3 -H 4 folate. Thus, it can be suspected that the heptaglutamate form accounts for the preferred folate polyglutamate form in yeasts and thus is more prone to oxidation reactions. Furthermore, the shift of the main polyglutamate from the hexa-to the heptaglutamate occurred during 29 and 32 h of cultivation for each of the polyglutamates analyzed by UHPLC-Q-ToF-MS (Supplementary Figure 3A). This time point represents the beginning of the stationary growth phase which is also characterized by decreasing need for folates in general. It is known that longer polyglutamates show higher affinity towards folate binding proteins (FBP) but reduce the activity of the latter (27). Consequently, it can be suspected that hexaglutamates reduce the activity of FBPs tremendously. Only a reduced number of folate monoglutamates competing with these hexaglutamates to be bound and metabolized by FBPs, enables a further elongation of the hexaglutamate tail. The reasons for the further elongation, however, are still unknown.
At the end of the exponential growth phase, oxidative degradation products were affected similarly by this shift in the polyglutamate tail length. However, ratios between hexa-and heptaglutamate form of these degradation products (shown in Supplementary Figure 3B) did not shift towards the heptaglutamate form to the same extent compared to the metabolic active vitamer forms. Three reasons for this observation are plausible: (1) longer metabolic active polyglutamates tend to be more stable (28) and a relatively higher proportion of hexaglutamates is being oxidized. (2) with elongation of the polyglutamate tail, metabolic inactive polyglutamates tend to be less stable (28) once having been oxidized and are more susceptible towards further oxidation reactions. (3) the observed results are the results of generally low abundancies of oxidation products compared to the corresponding active folate vitamers in combination with reduced response of longer polyglutamates during mass spectrometric analysis. Theory (1) is likely to be rejected as total MeFox contents were mostly deriving from 5-CH 3 -H 4 PteGlu 7 contents which seemed to be less stable compared to the hexaglutamate form. Theories (2) and (3) are somewhat difficult to prove. However, it is more likely that both theories are contributing to the observed experimental results.
To identify similarities between obtained MS/MS spectra acquired by UHPLC Q-ToF-MS, molecular networking was performed (17). MN revealed several features mostly referring to the main occuring hexa-and heptaglutamate derivatives. Clustering features consisted up to >60% of known folate vitamers. Further known degradation products than the above mentioned, quantified oxidation products were not detected within the MS 1 features. Detailed analysis of unknown but promising features throughout the whole sample set unraveled typical mass differences of CO, CH 2 , and H 2 O between already known and putatively as folate vitamer annotated species. The appearance of mass differences corresponding to the addition of a further methyl-or formyl group suggested the hypothesis of folate vitamers not only serving as C 1 -but also as C 2donors. Within the postulated C 2 -metabolism, only 5,10-CH + -CH 3 -H 4 folate could be found yet not finally identified as hexaand heptaglutamate. A further group of compounds consisting of hexa-and heptaglutamates could neither be identified as 5-CO-CH 3 -H 4 folate nor as the oxidation product of 5-CH 2 -CH 3 -H 4 folate, namely EthylFox. Solely, an annotation as C 21 H 25 N 7 O 7 could be made for the corresponding monoglutamate. Thus, this class of compounds requires further approaches to identify its structure. Any further C 2− reaction partners such as 5,10-CH 2 -CH 3 -H 4 folate or 5-CH 2 -CH 3 -H 4 folate could also not be found within the MS 1 features remaining the C 2 -metabolism questionable.
The detected 5,10-CH + -CH 3 -H 4 folate could have further origins other than the already discussed acetate. In accordance with the folate metabolism in yeasts, 5,10-CH + -H 4 folate is produced via three different pathways: (1) FTHFS (formyltetrahydrofolate synthetase) transfers formiate to the coenzyme H 4 folate forming 10-acetyl-tetrahydrofolate which will subsequently be metabolized into 5,10-CH + -H 4 folate, (2) SHMT (serine hydroxymethyl transferase) converts serine and H 4 folate into 5,10-CH 2 -H 4 folate which will subsequently be transformed into 5,10-CH + -H 4 folate, (3) the GCS (glycine cleavage systems) metabolizes H 4 folate and glycine yielding 5,10-CH 2 -H 4 folate which is further oxidized to 5,10-CH + -H 4 folate. In analogy, 5,10-CH + -CH 3 -H 4 folate could be traced back to acetate, homoserine, and alanine. Acetate occurs in yeasts ubiquitously as it is a by-product of alcoholic fermentation (29). Homoserine represents an intermediate of the synthesis of methionine which itself is involved in methylation reactions within the organism as well as the regeneration of H 4 folate (30, 31). Another key player in the folate metabolism, namely glutamate, is the product of the transamination of alanine and transfer of its amine function onto α-ketoglutaric acid yielding pyruvate (32). Consequently, all three precursors are involved in the metabolism of yeast with homoserine and alanine in the broad sense even being involved in the folate metabolism. Which one of the precursors does indeed account for the incorporation of the C 2 -moiety into the folate metabolism of yeasts needs to be addressed in future studies. Furthermore, more proof for the existence of the C 2 -metabolism needs to be found as well as its importance and function within the metabolism of yeasts. Last, the structure of the compound class annotated as C 21 H 25 N 7 O 7 in the respective monoglutamate form needs to be elucidated. Therefore, different foodstuffs need to be screened to identify . /fnut. . organisms with high concentrations of this compound which could then be used for further purification and identification.

Conclusion
We could show that 5-CH 3 -H 4 folate and H 4 folate monoglutamates increased until a cultivation time of 20 h and 26 h, respectively. 5-CHO-H 4 folate, however, increased throughout the whole cultivation process analyzed. Interestingly, the sum of 5-CH 3 -H 4 folate and 5-CHO-H 4 folate in relation to the total folate content remained constant from 20 h to 80 h, both vitamers were quantified up to the heptaglutamate. 5,10-CH + -CH 3 -H 4 folate could be found as hexa-and heptaglutamate but not finally identified as such. A further folate like group of metabolites was found as hexa-and heptaglutamate which according to the sum formulae obtained could also belong to a potential C 2 -metabolism. The latter group of folates could neither be identified as 5-acetyl-tetrahydrofolate nor as the oxidation product of 5-ethyl-tetrahydrofolate, here named as EthylFox. To identify the unknown group of folate metabolites, further studies need to be conducted.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions
Conceptualization: LSc, PS-K, and MR. Methodology: LSc, DW, NW, FS, and LSt. Data curation: LSc and DW. Writing: LSc. Review and editing: PS-K and MR. Funding: MR. All authors contributed to the article and approved the submitted version.

Funding
We thank the TUM Open Access Publication Fund and the TUM School of Life Sciences for supporting the green open access publication of this study.


Synonymous Codon Pattern of Cowpea Mild Mottle Virus Sheds Light on Its Host Adaptation and Genome Evolution

Cowpea mild mottle virus (CpMMV) is an economically significant virus that causes severe disease on several legume crops. Aside from recombination, other factors driving its rapid evolution are elusive. In this study, the synonymous codon pattern of CpMMV and factors shaping it were analyzed. Phylogeny and nucleotide composition analyses showed that isolates of different geography or hosts had very similar nucleotide compositions. Relative synonymous codon usage (RSCU) and neutrality analyses suggest that CpMMV prefers A/U-ending codons and natural selection is the dominative factor that affects its codon bias. Dinucleotide composition and codon adaptation analyses indicate that the codon pattern of CpMMV is mainly shaped by the requirement of escaping of host dinucleotide-associated antiviral responses and translational efficiency.

Introduction
Cowpea mild mottle virus (CpMMV), a whitefly (Bemisia tabaci)-transmitted virus of the genus Carlavirus within the family Betaflexviridae, has filamentous particles of about 650 nm in length and a diameter of about 13 nm [1]. The genome of CpMMV is composed of a positive-sense, single-stranded RNA (+ssRNA) that is capped at the 5 -end and polyadenylated at the 3 -end. The genome of CpMMV encodes a total number of six open reading frames (ORFs). The largest ORF on the 5 -termini encodes the viral RNA-dependent RNA polymerase (RdRp). ORF2-4 are partially overlapped and are composed of the so-called triple-gene-block (TGB) module, which encodes TGBp1, TGBp2, and TGBp3, respectively. ORF5 encodes the viral coat protein (CP), while ORF6 encodes a small cysteine-rich protein that may function as the pathogenicity factor and RNA silencing suppressor [2][3][4]. ORF1 is translated from genomic RNA, while ORF2-4 and ORF5-6 are translated from two subgenomic RNAs [5].

Phylogenetic Analysis
A phylogenetic tree was constructed using the MEGA 11 software with the neighborjoining (NJ) algorithm [23]. The nucleotide substitution model, mutation rate, and pattern were determined using the Model Selection function in MEGA 11. The bootstrapped confidence interval was based on 1000 replicates.

Calculation of the Relative Synonymous Codon Usage (RSCU)
RSCU is defined as the ratio of observed to expected codon frequency under equal codon usage without being affected by the amino acid compositions or the CDS sizes of different gene samples [24]. Synonymous codons with RSCU values lower than 1, equal to 1, and higher than 1 represent negative, no, and positive codon usage bias, respectively. Synonymous codons with RSCU values > 1.5 and <0.6 can be recognized as significantly overrepresented and underrepresented codons, respectively [25]. The RSCU of CpMMV was calculated with a homebuilt BioPython script. To calculate the RSCU values of host, open reading frames (ORFs) within mRNA sequences were extracted using a homemade BioPython script, which also removed any incomplete and erroneous sequences. The RSCU values of soybean, cowpea, and papaya were calculated by a BioPython script from a total number of 89,784, 41,014, and 25,430 full-length coding sequences, respectively. RSCU values of soybean leaf-specific and seed-specific genes were derived from [26], which were calculated from 58 and 21 highly expressed leaf-specific and seed-specific genes, respectively. All BioPython scripts are available upon request.

Nucleotide Composition Analysis
The ORFs in each CpMMV genome were retrieved manually and then were aligned by Clustal W 2.0 with default parameters [27]. Overall nucleotide composition (%A, %C, %U, and %G), nucleotide composition at the third codon position (%A3s, %C3s, %U3s, and %G3s), and mean frequencies of G + C nucleotide at the first, second, and third codon position (GC1s, GC2s, and GC3s) were calculated using the CodonW software version 1.4.2 with default parameters, which is available at http://codonw.sourceforge.net (accessed on 1 March 2022).

Analysis of Effective Number of Codons (ENC)
ENC is a simple and efficient method to quantify codon bias independent of gene length and amino acid composition [28]. The value of an ENC analysis ranges from 20 (when only one synonymous codon is used for the corresponding amino acid) to 61 (when all synonymous codons are used equally). The codon usage patterns were investigated by the NEC plot, a plot of NEC vs. GC3s at synonymous sites. In an ENC plot, the observed and expected ENC values are compared to determine the influence of structuring synonymous codon usage bias. An ENC plot is commonly used to determine the effect of G + C compositional constraints on codon usage bias. When the corresponding points fall near the expected curve, the mutation is the main force shaping codon usage. Alternatively, natural selection is the main force shaping codon usage when the corresponding points fall considerably below the expected curve.

Neutrality Plot Analysis
A neutrality plot was used to identify the effects of natural selection and mutation pressure on the codon pattern [29]. In a neutrality plot, the average GC content at GC12s was plotted against GC3s. A regression line was plotted between the GC3s-variable and the GC12s-variable. Unlike GC3s, GC1 and GC2 are subject to functional constraints, because a mutation at these positions usually leads to an amino acid change. As a result, the evolutionary speed of the mutation and natural selection pressure is expressed as the slope of the regression line, while the regression coefficient implies the mutation-selection equilibrium coefficient [30]. Theoretically, the mutation is recognized as the main force shaping the codon pattern when the slope of the regression line is close to 1, while natural selection is accepted as the main force when the slope of the regression line is close to 0.

Dinucleotide Odds Ratio
The dinucleotide odds ratio was calculated as described using a homebuilt BioPython script that is available upon request [31]. Based on the statistical theory, the dinucleotide can be recognized as significantly underrepresented if the ratio is equal or less than 0.78 and significantly overrepresented if the ratio is equal or higher than 1.23 [31].

Phylogenetic Analysis
To gain insight into the phylogenetic relationship of CpMMV isolates, a neighborjoining (NJ) phylogenetic tree was constructed using thirty-three full-length genomes of CpMMV (Table 1). Consistent with the previous study [32], the phylogenetic tree was apparently clustered into two major clades, namely Clade I and Clade II ( Figure 1A). Clade I could be further separated into three subclades, namely subclades I-1, I-2, and I-3. Subclade I-1 contained thirteen isolates from China; subclade I-2 contained eleven isolates from America and a papaya isolate from Kenya (MK984605); subclade I-3 contained two isolates from India and Sudan, respectively. Clade II could be further divided into two subclades: subclade II-1 contained two isolates from Ghana; subclade II-2 contained three isolates from Pakistan. We further analyzed the phylogenetic relationship of the thirty-three CpMMV isolates based on hosts that were isolated. We found that isolates from soybean and cowpea were located in Clades I and II; isolates from other legumes and non-legume hosts were clustered in Clade I ( Figure 1A). The GenBank accession number, origin, and host of each sequence are indicated. Isolates from soybean, other legumes, and non-legumes are in blue, purple, and black, respectively. Clades I and II are also highlighted in light blue and light green backgrounds, respectively.
Since recombination is prevalent in CpMMV isolates [17][18][19], which may mask the real phylogenetic relationship of these isolates; therefore, we constructed a phylogenetic tree based on the CP gene, which has a low recombination ratio. A total number of 161 sequences of the CP gene were retrieved from the GenBank database. For clarity, sequences with higher than 99% nucleotide sequence similarities were removed. A phylogenetic tree was then constructed using the rest 51 sequences of the CP gene ( Figure 1B). Consistent with the NJ tree based on full-length genomes, the 51 sequences of the CP gene were also clustered into two major Clades; however, the two isolates of subclade II-1 in the NJ tree constructed based on the full-length genome was comprised of the fourth subclade in Clade I. Geographic distribution analysis showed that isolates from India display a high degree of diversity since they were present in Clade I and Clade II. In this CP-based NJ tree, isolates from soybean were distributed in Clades I and II, isolates from cowpea were located in Clade I, and isolates from other legumes and non-legume hosts were found in Clades I and II ( Figure 1B).

Nucleotide Composition Analysis
To gain insight into the genetic forces driving the evolution of CpMMV, we first analyzed the nucleotide composition of the full-length genome of the thirty-three CpMMV isolates (excluding the poly-A tail). Results showed that the overall nucleotide composition of CpMMV was T and U-rich; U and C were the most abundant and infrequent bases, respectively (Table 2). We further calculated the nucleotide composition in the proteincoding regions of the thirty-three CpMMV isolates. The GC1 position had the highest G + C content in the CpMMV coding sequence (0.41), while the GC2 position had the lowest G + C content (0.40). Further analysis of the third wobble codon position revealed a significantly higher abundance of U3 (0.41) and A3 (0.35) than G3 (0.28) and C3 (0.24) ( Table 2). These results suggest that CpMMV displays a codon pattern that prefers A/U-ending codons; and, the codon pattern of CpMMV may be affected by compositional constraints. CpMMV had an overall effective number of codons (ENC) value of 52.23, suggesting a moderate but highly conserved codon pattern bias of CpMMV. The variations of nucleotide composition and ENC among isolates of Clade I and Clade II or of hosts were ignorable ( Table 2), indicating that different isolates have similar mutation bias and natural selection.

Natural Selection Is the Major Force Influencing the Codon Usage of CpMMV
To evaluate the forces shaping the codon pattern patterns of CpMMV, an ENC plot and neutrality analysis were carried out. In the ENC plot, all CpMMV isolates were clustered together and located below the expected curve, suggesting that the codon pattern of Cp-MMV is subjected to GC compositional constraints (Figure 2A). A detailed survey showed that isolates of Clade I and Clade II of different hosts were distributed in the same area ( Figure 2B,C), further confirming that different CpMMV isolates are subjected to similar GC compositional constraints. In the neutrality analysis ( Figure 2D), average values of GC3s were negatively correlated with the average values of GC12s (r = 0.715; p < 0.01); however, the slope of the regression line was 0.3863, which suggests that mutation pressure accounts for approximately 38.63% of variations and natural selection accounts for 61.37% of variations. Neutrality analysis was also performed according to their phylogenetic relationships and hosts. Results showed that the average values of GC3s were also negatively correlated with the average values of GC12s ( Figure 2E,F) and slopes values were closer to 0 than to 1. Thus, natural selection instead of mutation pressure plays a dominant role in shaping the codon pattern of CpMMV.

RSCU Patterns of CpMMV
To gain insight into the codon pattern of CpMMV, we calculated the overall RSCU values of the thirty-three CpMMV isolates (Table 3). Among the 59 synonymous codons, 9 and 12 were identified as significantly preferred codons (RSCU value > 1.5) and unpreferred codons (RSCU value < 0.6), respectively. An analysis of the base composition at the third codon showed that 7 out of 9 significant codons and 24 out of 27 preferred codons (RSCU value > 1.0) were A or U-ending, while 9 out the 12 significantly unpreferred codons were G or C-ending, suggesting that CpMMV prefers A/U-ending codons.
however, the slope of the regression line was 0.3863, which suggests that mutation pressure accounts for approximately 38.63% of variations and natural selection accounts for 61.37% of variations. Neutrality analysis was also performed according to their phylogenetic relationships and hosts. Results showed that the average values of GC3s were also negatively correlated with the average values of GC12s ( Figure 2E,F) and slopes values were closer to 0 than to 1. Thus, natural selection instead of mutation pressure plays a dominant role in shaping the codon pattern of CpMMV.

RSCU Patterns of CpMMV
To gain insight into the codon pattern of CpMMV, we calculated the overall RSCU values of the thirty-three CpMMV isolates (Table 3). Among the 59 synonymous codons, 9 and 12 were identified as significantly preferred codons (RSCU value > 1.5) and unpreferred codons (RSCU value < 0.6), respectively. An analysis of the base composition at the third codon showed that 7 out of 9 significant codons and 24 out of 27 preferred codons  We further analyzed the RSCU of CpMMV based on their phylogenetic relationship. The average RSCU values of Clades I and II were very close to that of overall RSCU values (Table S1). Especially, the significantly preferred and unpreferred codons of Clades I and II were identical to that of the overall RSCU values ( Figure 3A; Table S1). The RSCU values of the thirty-three CpMMV isolates according to their hosts were also calculated (Table S1). Noticeably, RSCU values of several synonymous codons, e.g., CCU, CCU, UGU, UGC, UUA, and UUG, showed high variations ( Figure 3B). Specifically, CCA had the highest range (R = 0.734), followed by CCU (R = 0.689), UGU/UGC (R = 0.612), and UUA (R = 0.528) ( Figure 3C). The extremum RSCU values were observed between soybean, cowpea, mungbean, and papaya-infecting CpMMV isolates, indicating that the codon pattern of CpMMV may be strongly affected by these hosts.

Dinucleotide Frequency Affects the Codon Pattern of CpMMV
Since the codon pattern of CpMMV is subject to GC compositional constraint and previous studies showed that several pairs of dinucleotides in RNA viruses were strongly subjected to host selection [31,33], the dinucleotide frequency of CpMMV was analyzed to gain insight into the constraint of dinucleotide bias on the codon pattern. Results showed that CpG and UpA were significantly underrepresented (observed/expected ratio ≤ 0.78) in both the coding region and untranslated regions (UTRs), while no dinucleotide pairs were found to be significantly overrepresented (observed/expected ratio ≥ 1.23) (Figure 4), suggesting that the codon pattern of CpMMV may be affected by CpG and UpA dinucleotide bias. A direct comparison of the RSCU values of codons with an underrepresented dinucleotide showed that the codon pattern is highly affected by the dinucleotide frequency. For instance, four of the six synonymous Arg containing CpG on the first and second codon position (CGN; N indicates A, U, G, or C) were underrepresented ( Table 2). Codons containing CpG on the second and third codon position (ACG (coding Thr), UCG (coding Ser), GCG (coding Ala), and CCG (coding Pro)) were also significantly unpreferred (Table 3). Although UpA is contained in two stop codons, we still could find the RSCU values of three out of the four codons containing UpA on the second and third codon position (AUA (coding Ile), GUA (coding Val), and CUA (coding Leu)) were lower than expected ( Table 2).

Codon Adaptation Analysis
To further gain insight into the impact of the host on CpMMV codon pattern, the overall codon pattern of soybean, cowpea, and papaya were calculated using all available nuclear genes in the GenBank nucleotide database (mungbean was not included due to data availability; Table S2). Corresponding analyses were then performed between RSCU values of CpMMV isolates from soybean (CpMMV-soybean), CpMMV isolates from cowpea (CpMMV-cowpea), or CpMMV isolates from papaya (CpMMV-papaya) and that of soybean, cowpea, and papaya. As expected, the codon pattern of CpMMV-cowpea was most adapted to cowpea according to the correlation coefficient ( Figure 5A-C). Interestingly, the codon patterns of CpMMV-papaya and CpMMV-soybean also displayed the highest adaptation to that of cowpea but not papaya or soybean ( Figure 5D-I), indicating that cowpea might be the original host of CpMMV. Compared to soybean, the codon pattern of CpMMV-papaya was more similar to that of papaya ( Figure 5D,E), suggesting that the infection of papaya causes the readaptation of the codon pattern from cowpea to papaya.

Dinucleotide Frequency Affects the Codon Pattern of CpMMV
Since the codon pattern of CpMMV is subject to GC compositional constraint and previous studies showed that several pairs of dinucleotides in RNA viruses were strongly subjected to host selection [31,33], the dinucleotide frequency of CpMMV was analyzed to gain insight into the constraint of dinucleotide bias on the codon pattern. Results showed that CpG and UpA were significantly underrepresented (observed/expected ratio ≤ 0.78) in both the coding region and untranslated regions (UTRs), while no dinucleotide pairs were found to be significantly overrepresented (observed/expected ratio ≥ 1.23) (Figure 4), suggesting that the codon pattern of CpMMV may be affected by CpG and UpA second codon position (CGN; N indicates A, U, G, or C) were underrepresented ( Table 2). Codons containing CpG on the second and third codon position (ACG (coding Thr), UCG (coding Ser), GCG (coding Ala), and CCG (coding Pro)) were also significantly unpreferred (Table 3). Although UpA is contained in two stop codons, we still could find the RSCU values of three out of the four codons containing UpA on the second and third codon position (AUA (coding Ile), GUA (coding Val), and CUA (coding Leu)) were lower than expected (Table 2).

Codon Adaptation Analysis
To further gain insight into the impact of the host on CpMMV codon pattern, the overall codon pattern of soybean, cowpea, and papaya were calculated using all available nuclear genes in the GenBank nucleotide database (mungbean was not included due to data availability; Table S2). Corresponding analyses were then performed between RSCU values of CpMMV isolates from soybean (CpMMV-soybean), CpMMV isolates from cowpea (CpMMV-cowpea), or CpMMV isolates from papaya (CpMMV-papaya) and that of soybean, cowpea, and papaya. As expected, the codon pattern of CpMMV-cowpea was most adapted to cowpea according to the correlation coefficient ( Figure 5A-C). Interestingly, the codon patterns of CpMMV-papaya and CpMMV-soybean also displayed the highest adaptation to that of cowpea but not papaya or soybean ( Figure 5D-I), indicating that cowpea might be the original host of CpMMV. Compared to soybean, the codon pattern of CpMMV-papaya was more similar to that of papaya ( Figure 5D,E), suggesting that the infection of papaya causes the readaptation of the codon pattern from cowpea to papaya. The codon pattern of CpMMV-soybean displayed a similar correlation coefficient to the codon pattern of soybean and papaya ( Figure 5G-I), implying that papaya and soybean have similar codon usage or CpMMV-soybean possibly originated from papaya.
To further gain insight into the codon adaptation, we compared the codon pattern of CpMMV with that of soybean leaf and seed. It has been reported that the codon usage of genes that were highly expressed in leaf tissues dramatically differed from that of genes expressed in seed [26]; therefore, it is expected that the codon usage of CpMMV is closer to that of leaf-specific genes than seed-specific genes since CpMMV is a phloem-limited virus. As expected, we found that the codon usage of CpMMV is highly correlated to that of leaf-specific genes (r = 0.883; p < 0.01), but not seed-specific genes (r = 0.402; p > 0.05) ( Figure 5); thus, the codon pattern of CpMMV is highly adapted to the host it infects.

Discussion
In this study, possible forces other than recombination that drive the evolution of CpMMV were analyzed using phylogeny, nucleotide composition, codon pattern, neutrality and ENC plot, dinucleotide bias, and host adaptation. Phylogenetic analyses showed that isolates from Africa and Southwest Asia were present in Clades I and II; isolates from East Asia and America were distributed only in Clade I (Figure 1), suggesting that CpMMV possibly originated from Africa and/or Southwest Asia and was introduced to East Asia and America through the international trade of agricultural products or along with the spread of its vector: silverleaf whitefly (Bemisia tabaci MEAM1). The availability To further gain insight into the codon adaptation, we compared the codon pattern of CpMMV with that of soybean leaf and seed. It has been reported that the codon usage of genes that were highly expressed in leaf tissues dramatically differed from that of genes expressed in seed [26]; therefore, it is expected that the codon usage of CpMMV is closer to that of leaf-specific genes than seed-specific genes since CpMMV is a phloem-limited virus. As expected, we found that the codon usage of CpMMV is highly correlated to that of leaf-specific genes (r = 0.883; p < 0.01), but not seed-specific genes (r = 0.402; p > 0.05) ( Figure 5); thus, the codon pattern of CpMMV is highly adapted to the host it infects.

Discussion
In this study, possible forces other than recombination that drive the evolution of CpMMV were analyzed using phylogeny, nucleotide composition, codon pattern, neutrality and ENC plot, dinucleotide bias, and host adaptation. Phylogenetic analyses showed that isolates from Africa and Southwest Asia were present in Clades I and II; isolates from East Asia and America were distributed only in Clade I (Figure 1), suggesting that CpMMV possibly originated from Africa and/or Southwest Asia and was introduced to East Asia and America through the international trade of agricultural products or along with the spread of its vector: silverleaf whitefly (Bemisia tabaci MEAM1). The availability of more sequences from other regions of the world could confirm or rule out this hypothesis. Codon adaptation analysis supported the fact that cowpea is the original host of CpMMV, while papaya may be an important intermediate host of CpMMV. Nevertheless, further studies are needed to confirm these possibilities.
Nucleotide composition, codon pattern, and ENC plot neutrality analysis showed that CpMMV had a codon pattern that prefers A/U-ending codons. Neutrality analysis showed that the RSCU of CpMMV is subjected to GC constraint and that natural selection instead of mutation pressure plays a dominant role in shaping the codon pattern of CpMMV. Dinucleotide frequency analysis showed that CpG and UpA were significantly underrepresented in both the untranslational regions and coding regions of CpMMV ( Figure 4). Interestingly, codons with underrepresented dinucleotides were always significantly unpreferred ( Table 2), suggesting that the RSCU of CpMMV is heavily affected by dinucleotide bias, especially CpG and UpA. The dinucleotide pattern is a stable genomic signature, which is related to genome compartmentalization, gene evolution, and gene expression regulation [34]. Studies showed that UpA is underrepresented in almost all organisms tested, whereas CpG is differentially represented in the genomes of eukaryotic organisms [35]. It was proposed that the TpA depletion may be caused by its presence in two out of three stop codons and in transcriptional regulatory motifs (e.g., the TATA box sequence), while CpG repression has been directly linked to cytosine methylation in DNA [36]. The dinucleotides CpG and UpA are also significantly underrepresented in the genome of retroviruses and riboviruses [31,33,37]. Importantly, increasing or decreasing the CpG and UpA dinucleotide frequencies in the genome could significantly attenuate or increase viral proliferation, respectively [38][39][40]. Recently, it was found that the attenuation of viruses with elevated frequencies of CpG and UpA dinucleotides in animals was associated with the ZAP and OAS3/RNAseL-mediated antiviral pathways, respectively [41,42]. How CpG and UpA affect the codon pattern of plant RNA viruses such as CpMMV is an interesting question for investigation in the future since neither ZAP nor OAS3/RNAseL counterparts have been identified in plants.
It has been proposed that the translation of mRNA with rare codons is less efficient than those with optimal codons since the latter have more cognate tRNA species [43]. As a result, a higher translation efficiency certainly will benefit the proliferation of viruses. Indeed, studies showed that the codon pattern of many viruses is highly adapted to their hosts [21,44,45]. Our results also showed that the codon pattern of CpMMV was overall highly adapted to its hosts, e.g., cowpea, papaya, and soybean. Moreover, our results showed that the codon pattern of soybean-isolated CpMMV and papaya-isolated CpMMV was closer to that of soybean and papaya, respectively ( Figure 5), which further confirmed the impact of the host on the codon pattern of CpMMV. Moreover, the codon usage of CpMMV resembles the codon usage of leaf-specific genes but not that of seed-specific genes ( Figure 6). Thus, mimicking the host codon pattern for maximum translation efficiency is also an important force shaping the codon pattern of CpMMV. highly adapted to its hosts, e.g., cowpea, papaya, and soybean. Moreover, our results showed that the codon pattern of soybean-isolated CpMMV and papaya-isolated CpMMV was closer to that of soybean and papaya, respectively ( Figure 5), which further confirmed the impact of the host on the codon pattern of CpMMV. Moreover, the codon usage of CpMMV resembles the codon usage of leaf-specific genes but not that of seedspecific genes ( Figure 6). Thus, mimicking the host codon pattern for maximum translation efficiency is also an important force shaping the codon pattern of CpMMV.

Conclusions
The synonymous codon pattern of CpMMV is shaped by the host it infects, possibly through pressures from dinucleotide bias and translational machinery.


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/.