Epigenetic Memory via Concordant DNA Methylation Persists in Mammalian Toti- and Pluripotent Stem Cells

In storing and transmitting epigenetic information, organisms must balance the need to maintain information about past conditions with the capacity to respond to information in their current and future environments. Some of this information is encoded by DNA methylation, which can be transmitted with variable fidelity from parent to daughter strand. High fidelity confers strong pattern matching between the strands of individual DNA molecules and thus pattern stability over rounds of DNA replication; lower fidelity confers reduced pattern matching, and thus greater flexibility. • We are interested in the strategies that various cell types, organisms, and species use to achieve balance between flexibility and stability. Here, we present a new conceptual framework, Ratio of Concordance Preference (RCP), that uses double-stranded methylation data to quantify the flexibility and stability of the system that gave rise to a given set of patterns. • We confirm previous observations that differentiated cells in mammals operate with high DNA methylation stability. Stem cells in culture and in embryos, in contrast, operate with reduced, albeit significant, methylation stability. We conclude that preference for concordant DNA methylation is a consistent mode of information transfer, and thus provides epigenetic stability across cell divisions, even in stem cells and those undergoing developmental transitions. Broader application of our RCP framework will permit comparison of epigenetic-information systems across cells, developmental stages, and organisms whose methylation machineries differ substantially or are not yet well understood.

Organismal development is characterized by a shift from the phenotypic flexibility of 2 embryonic cells to the canalized identities of differentiated cells. To achieve stable 3 gene-regulatory states in terminally differentiated cells, organisms ranging from Archaea to 4 humans use a variety of epigenetic mechanisms, including DNA methylation. Perturbation of the 5 state of DNA methylation at various loci in differentiated cells is associated with several human 6 cancers [1][2][3]. In turn, restoring epigenetic flexibility of some loci has proven challenging in 7 efforts to create induced pluripotent stem (iPS) cells [4]. Together, these findings highlight the 8 importance of shifting ratios of epigenetic flexibility and stability in establishing cellular 9 identity. 10 There exists an extensive literature documenting changes in single-locus and genome-wide 11 methylation frequencies at various stages of development [5,6]. Most genomic regions in 12 primordial germ cells (PGCs), for example, are known to undergo dramatic and rapid shifts in 13 DNA methylation frequency [7]. It is now clear that mammalian stem cells can utilize active 14 demethylation [8], highlighting the potential for both gain and loss of cytosine methylation to 15 impact the overall methylation frequency and, perhaps, stability of a given genomic region 16 during development. 17 High concordance of methylation in differentiated cells, with matching states for parent and 18 daughter DNA strands at individual CpG/CpG dyads, is considered to be a hallmark of 19 conservative epigenetic processes [9][10][11][12][13]. For earlier stages of development, the extent of 20 concordance is far less clear. For example, do methylation patterns in dividing embryonic stem 21 cells arise entirely by random placement of methyl groups, or is concordance favored to some 22 degree? 23 Recent work has begun to address these issues [7,[14][15][16][17][18]. Shipony et al. [16] analyzed 24 single-stranded DNA methylation patterns in populations of cultured cells established from 25 single founder cells. Under this approach, the degree of stability was inferred from the extent of 26 congruence among patterns collected from cultured descendant cells. The observation of 27 substantial pattern diversity among cells separated by many rounds of division led Shipony et 28 al. [16] to conclude that the bulk of methylation in human embryonic stem (ES) and induced 29 pluripotent stem (iPS) cells arises through "dynamic" -that is, non-conservative -DNA 30 methylation processes rather than through the "static" -that is, conservative -processes that 31 were emphasized in earlier studies [10,11,19]. Using data collected using hairpin-bisulfite 32 PCR [13], which yields double-stranded DNA methylation patterns, other studies suggested that 33 dynamic processes contribute substantially to DNA methylation in cultured mouse ES cells, but 34 perhaps not to the exclusion of the conservative processes that dominate at many loci in adult 35 differentiated cells [7,14,15,17,18]. 36 To fully characterize the balance between conservative and non-conservative methylation 37 processes, it is necessary to quantify the extent to which the arrangement of methylation in a 38 given set of patterns deviates from the null assumption of random placement. To assess and 39 visualize such deviations, we here introduce a new metric, Ratio of Concordance Preference 40 (RCP), which utilizes double-stranded methylation data. Here, as previously, we use the term 41 double-stranded DNA methylation pattern to refer to the overall pattern of methylation on both 42 top and bottom strands of an individual double-stranded DNA molecule. Double-stranded 43 patterns provide information on the extent of matching between methylation states on parent and 44 daughter strands, which are separated by exactly one round of DNA replication. RCP requires

73
Ratio of Concordance Preference is Defined for All Possible Configurations 74 of Methylation at Symmetric Nucleotide Motifs 75 We have developed Ratio of Concordance Preference (RCP) to assess the strategy of binary 76 information transfer, with focus on the degree to which exact information is conserved. We apply 77 our RCP framework to DNA methylation in mammalian cells. Our goal is to infer whether and 78 how much the system of processes that established a given set of methylation patterns prefers 79 concordant to discordant methylation states. This general formulation is free of assumptions 80 about the molecular mechanisms whereby methylation is added to and removed from DNA. 81 In our data from double-stranded DNA molecules from human and mouse, methylation 82 occurs principally at the CpG motif. This symmetric motif may be written as CpG/CpG, here 83 termed "CpG dyad". CpG dyads have opportunities for methylation on both strands. The 84 methylation state of a dyad thus takes one of three forms: fully methylated, at frequency M , 85 with methylated cytosines on both strands; hemimethylated, at frequency H, with a methylated 86 cytosine on only one strand; and unmethylated, at frequency U , with neither cytosine 87 methylated. The RCP framework can also be extended to non-CpG methylation at symmetric 88 nucleotide motifs. 89 To infer concordance preference for sets of double-stranded methylation patterns, we use the overall frequency of methylation, m, and the frequency of unmethylated dyads, U , of each data set. Because m is derived from the three dyad frequencies, the pair (m, U ) encompasses the full information available from the implicit dyad frequencies, M and H. We evaluate the extent of deviation from expectations under a random model in which the system has no preference for either concordant or discordant placement of methyl groups, using RCP, defined as: RCP can also be expressed in a form more familiar in biology if dyad frequencies are 90 considered as genotype frequencies for a gene with two alleles. RCP 2 is 4M U/H 2 , which is 91 expected to equal 1 under the Hardy-Weinberg equilibrium [21,22]. Thus, RCP can be 92 considered as a measure of deviation from random expectations.

93
The random expectations, for which RCP = 1, are met both with truly random placement of 94 methyl groups, and with equal contributions from processes operating with strong preference for 95 concordance and processes operating with strong preference for discordance. Under this random 96 model, the frequency of unmethylated dyads is given by U = (1 − m) 2 , leading to dyad 97 frequencies as expected under the binomial distribution (Fig 1a,b dashed curve; Fig 1c). A 98 system in which methyl groups are added de novo without regard to the methylation state of the 99 other strand [23], such as one dominated by the activity of mammalian Dnmt3s, is expected to 100 behave largely in accordance with random expectations. 101 One set of deviations from the random expectation is characterized by preference for 102 concordant placement of methyl groups, such that the two classes of concordant dyads -fully 103 methylated and fully unmethylated -are more frequent than expected under the random model. 104 This situation occurs under conservative systems of methylation where strong contributions from 105 maintenance-like processes, such as the activity of Dnmt1 in mammals [11,13,24], lead to high 106 frequencies of concordant dyads. In the extreme form of this deviation from random, methyl and unmethylated dyads (U ) locate each data set on the continuum from complete preference for concordance to complete preference for discordance. Ratio of Concordance Preference (RCP) is indicated for each contour line. (b) Expanded view. For this schematic, individual double-stranded methylation patterns (c-f) are used to illustrate different methylation configurations that lie along this continuum. Individual patterns, with methylated and unmethylated cytosines indicated in red and blue, respectively, can reflect (c) a random methylation system; (d) a fully conservative system, with complete preference for generating concordant dyads; (e) a fully dispersive system, with complete preference for generating discordant dyads (partial preference for dispersive placement is also possible); or (f) a partially conservative system, with more concordant dyads than expected under random processes, but fewer than expected under fully conservative processes. (g) A representative partial double-stranded DNA methylation pattern collected using hairpin-bisulfite PCR. The experiment-specific batchstamp is shown in green, and can be used to monitor for PCR contamination; the molecule-specific barcode shown in gray, generalized as "DDDDDDD", can be used to identify redundant sequences. The batchstamp and barcode are encoded by the hairpin oligonucleotide used to join the top and bottom strands. Primer-binding sites are underlined at the left end of the molecule. groups are observed only in fully methylated dyads (Fig 1d), such that unmethylated dyads 108 occur at frequency U = 1 − m (upper diagonal line in Fig 1a-b).

109
The other set of possible deviations from random is characterized by preference for 110 discordant placement of methyl groups, leading to an overabundance of hemimethylated dyads. 111 This situation occurs under dispersive systems of methylation such as those that yield transient 112 hemimethylation following DNA replication and prior to daughter-strand methylation, and 113 perhaps in genomic regions undergoing demethylation during periods of epigenetic transition.

114
When methylation is maximally dispersive and methylation frequency m is less than 0.5, all 115 dyads with methylation will be hemimethylated (Fig 1e), such that U = 1 − 2m (lower diagonal 116 line in Fig 1a,b); when m is greater than 0.5, not all methyl groups can be accommodated in 117 hemimethylated dyads, and so a combination of hemimethylated and fully methylated dyads -118 but no unmethylated dyads -is expected (lower horizontal line in Fig 1a).

119
The two extreme deviations from random form the boundaries of the comprehensive space of 120 possible configurations of methylation at symmetric motifs (Fig 1a). Sets of double-stranded 121 methylation patterns are expected to fall on the continuum between the extreme expectations 122 5/33 (Fig 1b), and can be located within this space to characterize the strategy of information transfer 123 employed to give rise to a given data set, ranging from conservative to dispersive.

124
As noted above, a system with an RCP value of 1 has no preference for either concordance 125 or discordance of methylation, and is analogous to the distribution of genotype frequencies at a 126 two-allele locus in a population that is at Hardy-Weinberg equilibrium [21,22]. An RCP value of 127 2 indicates two-fold preference for concordance, while an RCP of 1  activities [24,27]. 143 We apply RCP to investigate further the conclusion of Shipony et al. [16] that methylation in 144 cultured stem cells is dominated by non-conservative processes, with little or no preference for 145 concordance. Using double-stranded methylation patterns collected by our group, by Arand et 146 al. [14,17], and by Zhao et al. [15], we assess and compare methylation concordance in cultured 147 human and murine stem cells, as well as in murine cells undergoing early developmental 148 transitions that give rise to totipotent embryonic cells.

150
Our previous work with human single-copy loci in uncultured, differentiated cells revealed a 151 substantial role for maintenance methylation, a conservative process, with a comparatively 152 minor role for non-conservative de novo processes [24,27]. We therefore anticipated that RCP 153 analysis of double-stranded methylation patterns from such cells would indicate substantial 154 preference for concordant methylation states. Data published previously for G6PD, FMR1, and 155 LEP, in uncultured differentiated cells and new data presented here for FMR1 in cultured, 156 human differentiated cells represent blood, connective, and adipose tissues. We applied RCP 157 analysis to these data sets and found 13.2-to 85.7-fold preferences for concordant methylation. 158 This confirms, as anticipated, that methylation is predominantly conservative in these 159 differentiated cells (Fig 2a; see table accompanying Fig 2 for 95% CIs). We note that there is a 160 good correspondence between RCP and hemi-preference ratio, a statistic we computed for the 161 same data sets in the previous study [24] (further discussion in S2 Text). 162 We also found a substantial role for conservative methylation processes at single-copy loci in 163 both cultured and uncultured murine differentiated cells. Data sets from Arand et al. [14] for 164 Afp, Igf2, Snrpn, and Tex13 from murine embryonic fibroblasts (MEFs), and from Stöger [28] 165 for Lep from somatic tissues, gave RCP point estimates indicating a greater than 18-fold 166 preferences for concordant methylation (Fig 2b).  inferred to have strong contributions from conservative processes, using data sets that span a wide range of methylation frequencies. For each locus or multi-copy family, we inferred the RCP point estimate of the m, U pair and the two-dimensional 95% confidence region, determined by the uncertainty in the two variables (S7 Text). The intensity of coloration at a given point in a confidence region reflects the confidence level at that point. The m, U point estimates for most data sets are indicated with white asterisks, and the corresponding RCP values are given in the associated table. A larger, colored asterisk is used when the confidence interval of a data set is too small to be readily visible. RCP point estimates and 95% bias-corrected bootstrap confidence intervals are shown in figure-associated tables. (a) Three single-copy human loci -G6PD, FMR1, and LEP -and one human repeat family, L1, all from various tissues as indicated. (b) Four single-copy loci and four repeat families -Afp, Igf2, Snrpn, Tex13, B1, IAP, L1, and mSat -from murine embryonic fibroblasts, one single-copy murine locus -Lep -from somatic tissue, and one repeat family -L1 -from murine gametes. Data in (a) were collected using hairpin-bisulfite PCR and dideoxy sequencing, and taken from published [24,[27][28][29] and previously unpublished work (Table S3). Data in (b) were collected by Arand et al. [14] using hairpin-bisulfite PCR and pyrosequencing, with the exception of murine Lep for which data from somatic tissues were collected by Stoger [28], using dideoxy sequencing. Our analyses of these data sets applied bootstrapping approaches and accounted for inappropriate and failed conversion of methylcytosine using methods described in Supporting Information (S4 Text). Dyad counts, conversion-error rates, and inferences for methylation frequencies and RCPs are summarized in Tables S3, S4. Do multi-copy sequence families also have high preference for concordant methylation in 168 differentiated cells? We inferred RCP for four repeat families -B1, IAP, L1, and mSat -using 169 murine data collected by Arand et al. [14]. Three of these families were found to have preference 170 for concordant methylation in the same range inferred for single-copy loci (RCP point estimates 171 between 14.4 and 19.3; Fig 2b). The fourth -mSat -had an RCP estimate of 7.33, lower than 172 other families and single-copy loci examined in MEFs, but still indicative of strong preference 173 7/33 for concordant methylation. For human cells, data from two independent lines of cultured 174 embryonic fibroblasts were available for the repeat family L1. Inferred RCP values were within 175 the range found for single-copy loci in both human and murine differentiated cells (Fig 2a). 176 Overall, we find appreciable preference for concordance across a diverse group of data sets 177 from differentiated cells. These sets span a more than five-fold range in methylation frequency, 178 underscoring the independence of RCP from m, and, more generally, highlighting the capacity 179 of methylation systems to propagate specific epigenetic states, even when methylation is sparse. 180 We conclude that preference for concordant methylation, albeit to variable degrees, is present in 181 differentiated cells across broad classes of genomic elements, cell and tissue types, and culture 182 states.

183
Concordance Preference is Reduced but still Substantial in Cultured Stem 184 Relative to Differentiated Cells 185 We next ask whether substantial preference for concordance, as we infer above for 186 differentiated cells, is also evident in data from cultured stem cells. In doing so, we compare our 187 findings using RCP to the expectation from Shipony et al. [16] that methylation in such stem 188 cells occurs primarily through non-conservative, random processes.

189
The broadest data set available for our analysis comes from the near-genome-wide 190 double-stranded methylation data presented by Zhao et al. [15]. These data give an inferred RCP 191 of 5.22 for "all CpGs" in DNA from undifferentiated, cultured murine ES cells (Fig 3a; 192 Table S5). For other classes of genomic elements in these near-genome-wide data [15], we infer 193 RCP values of 4.31 or greater (Table S5). These RCP values are significantly higher than 1, the 194 value predicted under Shipony et al.'s proposal of dynamic methylation (p < 10 −16 , maximum 195 likelihood comparispon tests (MLCTs)). 196 We next ask whether our inference of appreciable concordance preference in the murine ES substantial preference for concordant methylation (Fig 3a). One single-copy locus, Afp, had a 203 methylation level too high, 0.99, to permit reliable inference of RCP. Murine double-stranded 204 methylation patterns for the four repeat families were available for two more stem cell lines, E14 205 and WT26 [14]. These additional repeat-family data sets, too, had RCP values significantly 206 greater than 1, although one data set, that for mSat in WT26, had an RCP value closer to 1 than 207 did others (p = 0.045, one-tailed BT). Data were available for a single-copy locus, Lep, for a 208 fourth murine ES line, CGR8. Here, too, RCP was significantly greater than 1 (p < 10 −16 , 209 one-tailed BT).

210
Human stem cell lines also followed the pattern of preference for concordant methylation.  Table S3). Together, these values reveal concordance 215 preference that is reduced relative to differentiated cells, but still greatly exceeds expectations 216 under random placement of methyl groups (p < 10 −16 , one-tailed BT). 217 We now consider the possibility that spontaneous differentiation had produced  stem cells were consistently inferred to have substantial contributions from conservative processes, with concordance greater than the random expectation. RCP point estimates and 95% biased-corrected bootstrap confidence intervals are shown for individual loci and "All CpGs". (a) Three single-copy loci -Igf2, Snrpn, and Tex13 -as well as four multi-copy loci -B1, IAP, L1, and mSat -from murine ES J1 cells were assayed by Arand et al. [14]. "All CpGs" data, collected by Zhao et al. [15]), reflect methylation at 17.3% of CpG dyads in the murine genome (Table S5). Data from both Arand et al. [14] and Zhao et al. [15] were collected using hairpin-bisulfite PCR and pyrosequencing. (b) Human L1 and LEP data were collected using hairpin-bisulfite PCR and dideoxy sequencing (data from published [28] and previously unpublished work (Table S3)). Our analyses of these data sets applied bootstrapping approaches and accounted for inappropriate and failed conversion of methylcytosine using methods described in Supporting Information (S4 Text). Dyad counts, conversion-error rates, and inferences for methylation frequencies are given in Tables S3, S4, and S5.
subpopulations of cultured stem cells that might account for the inference of RCP values 219 substantially greater than 1 at the seven different loci and genomic elements examined. Our 220 calculations revealed that a possible subpopulation of differentiated cells operating at much 221 higher RCP than that of undifferentiated cells would need to comprise more than 50% of the 222 population to account for our finding (S9 Text). Morphological inspection of the cultured human 223 stem cells under non-differentiating conditions did not suggest the presence of a substantial 224 subpopulation of differentiated cells in any of these lines. 225 We conclude that RCP values significantly greater than 1 are a consistent feature of cultured 226 embryonic stem cells, and hold across a broad set of stem cell lines, genomic locations and 227 element categories. 228

9/33
Preference for Concordance is Minimal or Absent in ES Cells Deficient in 229 Maintenance Methylation Activity 230 Our finding of substantial preference for methylation concordance in data from cultured, 231 undifferentiated stem cells contrasts with the inference of Shipony et al. [16] that DNA 232 methylation in such cells is dominated by non-conservative, random processes. This disparity 233 led us to ask whether our approach here for data acquisition and analysis is indeed capable of 234 identifying sets of methylation patterns established under exclusively random processes, which 235 are expected to yield RCP values of 1 (see "Defining the Ratio ...", above).

236
To assess this capacity, we consider methylation patterns from two murine embryonic stem 237 cell lines that have impaired maintenance methylation: a Dnmt1 knockout (KO) line and an  RCPs inferred for cultured murine ES cells were close to 1, the random expectation, at some loci. Data from Lep, collected using hairpin-bisulfite PCR and dideoxy sequencing, are from Al-Alzzawi et al. [30]. Data from seven additional loci are from Arand et al. [14], who used hairpin-bisulfite PCR and pyrosequencing. Our analysis of these published data revealed two of the eight loci analyzed -Afp and B1 -to have very low RCP values not significantly different from 1. (b) In the absence of Np95, a protein critical for recruiting Dnmt1 to hemimethylated regions in newly replicated DNA, one of the four loci analyzed -B1 -was inferred to have a very low RCP value, not significantly different from 1. Our analyses of these data sets from Arand et al. [30] applied bootstrapping approaches and accounted for inappropriate and failed conversion of methylcytosines using methods described in S4 Text. Point estimates and 95% confidence intervals on RCP are given above, and also along with conversion-error-rate estimates in Tables S3 and S4.

10/33
Np95 KO line. The Dnmt1 enzyme is principally responsible for addition of methyl groups to 239 daughter-strand CpGs complementary to CpGs methylated on the parent strand [11,24,31]; 240 Np95 facilitates interaction of Dnmt1 with these hemimethylated sites [32]. Absence of either 241 protein is therefore predicted to markedly diminish maintenance activity. If our approach is able 242 to detect essentially random placement of methyl groups, RCP values in these knockout lines 243 should be ∼1 for loci for which Dnmt1, its actions facilitated by Np95, is principally responsible 244 for conservative methylation.

245
Significant reductions in RCP were inferred for all single-copy loci and repeat families 246 examined in Dnmt1 and Np95 KO lines [14,30], compared to the parent stem-cell lines. Some 247 reductions were sufficient to bring RCP values in the knockout lines to that expected for random 248 placement of methyl groups: one single-copy locus -Afp -in the Dnmt1 KO line and one repeat 249 family -B1 -in both knockout lines had RCP values not significantly different from 1 (Afp in 250 Dnmt1 KO: 1.16, p = 0.10; B1 in Dnmt1 KO: 1.14, p = 0.17; B1 in Np95 KO: 1.02, p = 0.38; 251 one-tailed BTs; Fig 4). These findings in the two mutant cell lines thus reveal that RCP analysis 252 is, indeed, able to detect methylation established with random placement of methyl groups, and 253 thus with little or no preference for concordance or discordance.

254
The ability of RCP to detect random methylation has important implications for our work.

255
First, we can conclude that our inference of persistent preference for concordant methylation in 256 cultured stem cells reflects a bona fide property of those cells, rather than an artifact of our 257 approach. Second, we can infer from our finding of RCP > 1 for nine of the twelve data sets 258 examined in the Dnmt1 and Np95 KO lines that methyltransferases other than Dnmt1 can 259 contribute to conservative methylation. This inference is consistent with earlier conclusions that 260 contributions of Dnmt3s can include low levels of maintenance activity [19,24,33]. region inferred for single-copy loci in differentiated somatic cells (Fig 2 and Fig 5a). Using 273 near-genome-wide data for cultured murine cells [15], we inferred significant RCP increases 274 upon cell differentiation for most genomic elements (p < 10 −16 , MLCTs), with the exception of 275 low-complexity and satellite DNAs (Table S5). These RCP increases were greatest at promoters, 276 CG islands, and CG shores, and are more modest at other regions. We conclude that the onset of 277 differentiation in cultured human and murine cells is associated with a shift towards a greater 278 role for conservative processes. 279 Does the dedifferentiation that occurs in culture upon production of an iPS line from a 280 differentiated cell have an opposite effect on concordance preference? To address this question, 281 we compare methylation at L1 elements in three iPS lines to that in the two cultured human 282 fibroblast lines from which they were derived. As predicted, RCP values for all three iPS lines 283 were much reduced compared with values observed for the parent fibroblast lines (p < 10 −16 , 284  human fibroblast lines had substantial preference for concordance. Upon dedifferentiation to iPS cells, RCP values were reduced, indicating diminished preference for concordance. We present data for two different iPS lines established independently from the fibroblast line, FX. These iPS lines differed in their methylation frequency, m, at L1 elements, but had similar RCP values. Some data were previously shown in Fig 2 and Fig 3, and are included again here to illustrate more directly the relationships of L1 methylation in the differentiated and undifferentiated cultured cells. We applied bootstrapping approaches and accounted for inappropriate and failed conversion of methylcytosine using methods described in S4 Text. Point estimates of RCP, with 95% confidence intervals, are shown above, and also with dyad counts, conversion-error rates, and inferences for methylation frequencies in Tables S3  and S4. two-tailed PTs; Fig 5b). Dedifferentiation in tissue culture is thus associated with a shift in DNA 285 methylation toward a greater role for non-conservative processes. It will be useful to investigate 286 whether changes in methylation systems as measured by RCP drive or merely reflect the cellular 287 differentiation process. 288

12/33
Murine Embryos Mirror the Transitions in Concordance Preference Ob-289 served in Cultured Cells

290
The ∼3-fold-or-greater preference for concordant methylation we infer in many different 291 cultured ES and iPS cell lines (Fig 3a; Table S3) far exceeds the concordance expected under the 292 null hypothesis that methyl groups are placed at random. Here we ask whether this appreciable 293 preference for concordance is an artifact of growing stem cells in culture, or whether it is shared 294 by uncultured stem cells taken directly from an embryo. 295 We first consider whether totipotent cells from an embryo have evidence of conservative 296 processes. We applied RCP to double-stranded methylation patterns collected by Arand et 297 al. [17] for three multi-copy loci in mouse embryos: L1, mSat, and IAP. Our analyses revealed 298 that these totipotent embryonic cells, from post-replicative zygote to morula stage (through 3 299 days post conception, dpc), also exhibit moderate preference for concordance. Each of the 300 eighteen data sets we considered yielded RCP point estimates greater than 1, and confidence 301 intervals that exclude 1 (p < 0.005; Fig 6a; Table S4).

302
Pluripotent stem cells from mouse embryos (gastrula, 3.5 dpc) also exhibit moderate  Table S4). Thus, we conclude that moderate preference for concordance is an epigenetic 305 feature of uncultured embryonic stem cells of disparate developmental potential, and is not an 306 artifact of the establishment of embryonic stem cells in culture.

307
Though RCP values for stem cells from embryos clearly indicate some preference for 308 concordance, the extent of this preference is lower than for differentiated cells at most of the loci 309 examined (Fig 2 and Fig 5). Does this lower preference for concordance originate in gametes, or 310 does it instead arise post-fertilization? To address this question, we expanded our inference of 311 RCP values for sperm and oocytes from just L1 (Fig 3) to all three loci analyzed by Arand et 312 al. [17]. At the three multi-copy loci, RCP values for gametes are within the range observed for 313 other differentiated cells, with average point estimates ranging from 13.4 to 45.0 (Fig 6b; 314 Table S4). This high preference for concordance in gametes implies that the lower RCP values 315 characteristic of zygotes and stem cells must arise post-fertilization rather than in gametes.

316
To pinpoint the timing of this transition to lower RCP values observed in stem cells, we 317 consider data from post-fertilization nuclei and cells [17]. Data available for pronuclear stages 1, 318 2, and 3 revealed high RCP values, similar to those observed in gametes. In pronuclear stages  Table S4).

321
Is this transition dependent on the DNA replication event that occurs from pronuclear stage 3 322 to stages 4-5? To address this question, we assess data from Arand et al. [17], in which 323 aphidicolin was used to block DNA replication in the fertilized egg. Our RCP analysis reveals 324 that methylation patterns at L1 and mSat in these treated cells, while having somewhat lower 325 RCP values relative to pronuclei at earlier stages, had not undergone the major reduction in RCP 326 that we infer for unmanipulated pronuclei at stages 4-5 (p < 10 −16 , two-tailed PTs; Table S4). 327 Thus, the shift to lower RCP in stem cells following fertilization appears to require either DNA 328 replication or some later event that is itself replication-dependent. This conclusion is consistent 329 with the inference of Arand et al., using a different metric [17], that DNA replication in the      (4)(5) were accompanied by a sharp decrease in RCP at L1 elements, to a level similar to that observed in cultured stem cells (Fig 3). Further embryonic development was accompanied by minor increases and subsequent decreases in methylation and RCP values. (b) In PGCs at the earliest stage for which data are available, 9.5 days post conception (dpc), L1 elements had RCP values that were unusually low but still significantly greater than 1 (p < 10 −16 ). RCP values increased during PGC maturation to stage 13.5 dpc even as methylation frequencies decreased by more than 50%. RCP values for eggs and sperm, shown previously in Fig 3, are included here to highlight the transition from early primordial germ cells to terminally differentiated gametes. (c) Tracking RCP at repeat families during development. We inferred RCP for the embryonic and differentiated stages for which data were published by Arand et al. [17]. The topology of our longitudinal RCP plot highlights the transitions also evident in Arand et al.'s plot of the percentage of hemimethylated CpG dyads relative to all methylated CpG dyads (Figure 6b in [17]). Their metric captures relative shifts in concordance, but, in contrast to RCP, does not include a null model enabling quantitative comparison of inferred concordance across data sets with disparate methylation frequencies [17]. Point estimates and 95% confidence intervals for all but one data set were estimated by bootstrapping and applying the BCa correction (S7 Text). Dyad counts, methylation frequencies, and conversion error rates are in Table S3. In cases where data were available for multiple replicates, dyad frequencies were pooled. For one data set, 3 dpc at mSat, double-stranded sequences and methylation patterns were not available to us for bootstrap analysis; we therefore used a likelihood-based method for the estimation of confidence intervals (S8 Text).
14/33 divisions? Murine primordial germ cells (PGCs), in their maturation to differentiated gametes, 335 offer an opportunity to approach this question [34,35]. Using data from Arand et al. [17], we 336 infer that RCP values for PGCs increased by factors ranging from 2.5 to 5 during their 337 progression from 9.5 dpc to 13.5 dpc. This increase was not sudden, but occurred over the 338 four-day period, and so spanned an interval of substantial cell proliferation [34] (Fig 6b). The 339 murine embryo data thus provide examples of both abrupt and gradual transitions in RCP 340 through development (Fig 6c).

341
The RCP values for PGCs at 9.5 dpc, the earliest stage for which data were reported by 342 Arand et al. [17], were strikingly low: 1.40 for L1, 1.57 for mSat, and 2.38 for IAP. Nonetheless, 343 confidence intervals for all three loci indicated RCP values greater than 1 (p < 2 × 10 −5 ), the 344 value expected expected under wholly random placement of methyl groups, indicating persistent, 345 low-level preference for concordance. This low, residual preference for concordance in maturing 346 PGCs perhaps reflects both the epigenetic memory needed to maintain the poised state of stem 347 cells, and the epigenetic flexibility required for the production of differentiated gametes.

348
The first few rounds of PGC division involve developmental reprogramming and 349 commitment [35], and establishment of lineage-specific gene expression patterns. The 3.5-fold 350 average increase of RCP across these divisions (Fig 6;  across early rounds of division [36]. Seisenberger et al. [7], Arand et al. [17], and von Meyenn 357 et al. [18] concluded that this reduction in methylation frequency is driven by partial impairment 358 of maintenance methylation.

359
Our RCP framework permits a closer look at the likely extent of this proposed maintenance 360 impairment. If maintenance methylation were completely impaired and no other methylation 361 processes were active, the resulting passive demethylation would be fully dispersive: methyl 362 groups would remain only in hemimethylated dyads. A fully dispersive process such as this 363 should halve the methylation frequency at each DNA replication event and yield an RCP value of 364 0. We infer an extremely low value of RCP only for cells treated with S-Adenosylmethionine-ase 365 (0.20; Table S4). As this treatment either impairs or eliminates a cell's ability to methylate DNA, 366 it reveals RCP trajectories that would be observed with complete or nearly complete suspension 367 of all methylation processes. Our finding of RCP values greater than 0, as well as 1, in maturing 368 PGCs suggests not only that conservative processes remain, but also that their contributions are 369 sufficient to outweigh any dispersive effects of passive demethylation.

371
Because RCP makes no explicit enzymatic or mechanistic assumptions about the methylation 372 machinery, it permits quantification and comparison of strategies for symmetric methylation 373 across cell types, developmental periods, and organisms, despite potential and likely differences 374 in exact mechanisms. Application of RCP to double-stranded DNA methylation patterns reveals 375 that preference for concordance in DNA methylation is a persistent though quantitatively 376 variable feature of mammalian cells of disparate developmental potential. Specifically, we find 377 that: (i) in cultured human and murine ES and iPS cells, preference for concordance is lower 378 than in differentiated cells, but not absent; (ii) for cultured human stem cells, cellular 379 differentiation is characterized by increasing preference for concordance, whereas, for cultured 380 15/33 differentiated cells, dedifferentiation is characterized by declining preference for concordance; 381 and (iii) during early murine development, transitions in RCP mirror those found in cultured 382 cells, with pluripotent and totipotent stem cells showing appreciable concordance preference 383 throughout. We also observe that substantial changes in RCP can be either abrupt, requiring only 384 one DNA replication event, or gradual, occurring over multiple rounds of replication.

385
Although preference for concordance is substantial throughout early murine development, 386 there is an instance of concordance preference near the expectation under entirely random 387 processes. We infer RCP values close to, albeit significantly different from, 1 in the early 388 primordial germ cell stage at the three repetitive element families examined by Arand et al. [17] 389 (Fig 6). The instance of low, yet present, concordance preference may reflect both the epigenetic 390 stability required to maintain the poised state of the stem cells and the epigenetic flexibility 391 needed en route to production of functional gametes. Flexibility, indicated by RCP values near 1, 392 may result from near-random processes or instead from a balance of conservative and dispersive 393 methylation. Existing data and conclusions of Seisenberger et al. [7], Arand et al. [17] and von 394 Meyenn et al. [18] are more consistent with the latter interpretation.  al. [16], such as for stem cells dividing in culture, preference for concordant methylation may be 410 less important than other mechanisms of epigenetic memory. For example, regulation of 411 promoter activity by DNA methylation can occur via an ensemble effect rather than by 412 methylation of specific CG dyads within a promoter [37]. In such cases, propagation of exact 413 methylation patterns may be less important than the density of methylation that influences gain 414 or loss of methylation and states of transcriptional activity over many cell divisions [38]. 415 Epigenetic mechanisms other than DNA methylation also contribute to epigenetic memory at 416 various timescales. RCP analysis in combination with histone-modification data from 417 ENCODE [39] and Roadmap [40] will provide unprecedented opportunities to infer interactions 418 between DNA-methylation machinery and histone modification, the developmental timing of 419 epigenetic stability, and its variation across the genome.

420
The value of RCP analysis will be enhanced and broadened by emerging DNA sequencing 421 technologies that yield longer, more informative double-stranded methylation patterns. Longer 422 sequence reads will enable inference of RCP for single cells, permitting study of cell-cell 423 epimosaicism, such as arises in cancer [2] and other syndromes characterized by epigenetic 424 heterogeneity and change [29,41]. Some of these methods also reduce data corruption arising 425 through errors in bisulfite conversion and amplification, and can distinguish between methyl-426 and hydroxy-methyl cytosine [42,43]. High-resolution RCP estimates available through these 427 advances will provide new insight into the flexibility and potential sensitivity of individual loci 428 and cell types to environmental conditions encountered during embryogenesis and beyond. 429

430
Mathematical Foundations of the RCP framework 431 Overview of the mathematical foundation is given in Results and Discussion, and developed 432 further in S1 Text.

434
The six human ES and iPS cell lines for which we collected methylation patterns were 435 derived from either embryos or fibroblasts described as normal [44] or from fibroblasts of 436 individuals with disorders not known to affect the basic biochemistry of DNA methylation.

437
Many of the sets of human DNA methylation patterns analyzed here were presented in 438 previous publications, which include information on University of Washington Human Subjects 439 approval for collection and use. These data include G6PD and FMR1 from leukocytes of normal 440 individuals [24,27]; FMR1 from males with fragile X syndrome [29]; LEP from male leukocytes 441 and from female lymphocytes and adipose tissue [24,28]. Methylation patterns from murine ES cells, and the origin and culturing of these cells, have 464 previously been described [14,15,17,30]. 465

467
The DNA methylation patterns collected in our lab and analyzed here, both those published 468 previously and those presented here for the first time, were collected using the hairpin-bisulfite 469 PCR approach [13], with barcodes and batchstamps to authenticate each sequence [47]. Details 470 for collection of each data set are given in Table S2. 471 The data presented by Arand et al. [14,17], and Zhao et al., [15] and analyzed here, were 472 collected in the absence of molecular batch-stamps and barcodes, raising the possibilty that the 473 reliably of those data sets is undermined by PCR clonality. However, both groups used alternate 474 strategies that revealed that PCR clonality was not rampant. Zhao et al. found that essential 475 features of data sets did not differ appreciably between the "real" data set collected 476 conventionally, using PCR, and a test, PCR-free data set that excluded opportunities for clonality 477 by including only one read from each locus, providing no evidence of impacts from PCR     Table S3. RCP values and associated confidence intervals inferred for the 24 data sets from our labs. We collected double-stranded DNA methylation patterns from two speciesmouse and human -and several loci using bisulfite conversion under either low-molarity/temperature ("LowMT") or high-molarity/temperature ("HighMT") conditions [48]. For each data set, we counted methylated (M ) , hemimethylated (H), and unmethylated (U ) dyads, and used these values to infer methylation frequency, m, unmethylated dyad frequency, U , and the ratio of concordance preference, RCP. Inferences for (m, U ) and for RCP incorporated correction for conversion error, as described in S3 Text. We estimated 95% confidence intervals by bootstrapping and applied the BCa correction to obtain bias-corrected confidence intervals and point estimates, as described in S7 Text. Both uncorrected and BCa-corrected intervals and point estimates are listed here. Table S4. RCP values and associated confidence intervals inferred for data reported in Arand et al. [14,17]. Dr. Julia Arand kindly shared raw double-stranded sequences for samples described in these publications. We accounted for failed and inappropriate conversion, as described in S3 Text, in our point estimates of m, U , and RCP. The inappropriate conversion rate of 0.031 was inferred because the bisulfite conversion conditions used by Arand et al. resembled lowMT conditions [48]. We estimated 95% confidence intervals by bootstrapping and applied the BCa correction to obtain bias-corrected confidence intervals and point estimates, as described in S7 Text. Both uncorrected and BCa-corrected intervals and point estimates are listed here. For one sample, Arand et al. (2015) mSat 3dpc, only dyad counts, but no raw sequences, were available. We used the dyad counts to estimate the point estimate and the confidence interval for this sample while assuming independent sampling of dyads (S8 Text). Table S5. RCP values and associated confidence intervals inferred for data reported in Zhao et al. (2014) [15]. Dr. Hehuang Xie kindly shared raw data on M , H, and U values for samples described in Figure S6 of their paper, and also provided information on error rates from bisulfite conversion: for the "Day 0" sample (cultured stem cells grown under non-differentiating conditions), the failed conversion rate was 0.011, and the inappropriate conversion rate was 0.0109; for the "Day 6" sample (cultured cells grown for 6 days under differentiating conditions), the corresponding error rates were 0.012 and 0.0099. Dr. Xie also commented, "'All' refers to the total CG dyads and 'DNA' refers to the CG dyads within 'DNA repeat elements (DNA)' annotated in the UCSC genome database." For these data sets, whose individual sequence reads only contain 2 to 3 dyads on average, we assumed independent sampling of dyads to obtain the confidence intervals (S8 Text). Table S6. Pairwise comparisons of RCP values inferred for replicate samples of cultured human ES and iPS cells, and for mouse embryonic cells sampled at various developmental stages. We applied our test for heterogeneity, as described in S5 Text, to assess evidence of significant differences between RCP values inferred from various sample replicates presented by Arand et al. (2015) [17]. A color gradient encodes approximate levels of significance, with dark green indicating non-significant differences, and red indicating highly significant differences.

S1 Text: Deriving the Ratio of Concordance Preference, RCP
The goal of RCP is to quantify the extent to which the DNA methylation machineries that gave rise to each data set deviate from random expectations under the binomial distribution, as indicated by an over-or under-abundance of concordant dyads. Here we approach this problem by modeling the formation of concordant and discordant dyads as transitions between unmethylated and hemimethylated dyads and between hemimethylated and fully methylated dyads, without regard to the molecular processes that facilitate these transitions. We then take a mathematical approach to derive an expression for RCP.
We seek the equilibrium frequencies of fully methylated (M ), hemimethylated (H), and unmethylated (U ), dyads. Consider a continuous-time Markov chain operating on the probability distribution of the dyads M, H, U : where η's and γ's represent the rates of methylation addition and removal, respectively, as shown in Fig S1.

23/33
In the meantime, we define RCP as the ratio between β c , the rate of dyad conversions yielding concordant dyads, and β d , the rate of dyad conversions yielding discordant dyads. We thus define β c = √ η 1 γ 0 , the geometric mean of the methyl-addition and methyl-removal rates yielding concordant dyads. Likewise, we define β d = √ η 0 γ 1 .
We can then solve for the steady state distribution for the Markov chain in Equation (2) to arrive at the ratio. We can also express in in terms of m and U , as shown below.
This formulation of RCP does not require the assumption that methylation frequency is constant over time. Here, "steady state" refers to the dyad frequencies expected under a given system of methylation processes, regardless of whether the methylation frequency is constant over rounds of cell division.
It is notable that RCP 2 is 4M U/H 2 , which is expected to equal 1 under the Hardy-Weinberg equilibrium [21], if dyad frequencies are considered as genotype frequencies of a gene with two alleles. Following this, RCP can be considered as a measure of deviation from the null equilibrium.
RCP is therefore a metric for the degree to which the system of methylation processes prefers concordance (RCP > 1), discordance (RCP < 1), or, possibly, exhibits no preference in either direction (RCP = 1). If we set RCP = 1 and solve this expression for U , we find that U = (1 − m) 2 . This is consistent with the expectation under the binomial distribution that RCP will be 1 when there is no preference for either concordance or discordance. As RCP approaches ∞, U approaches 1 − m. Setting RCP = 0 results in two solutions: U = 1 − 2m and U = 0. These solutions are congruent with the boundaries that define the space of (m, U ) as given in "Defining the Ratio of Concordance Preference for All Possible Methylation Configurations" of the main text, and in Fig 1. S2 Text: Comparing RCP Values and Hemi-Preference Ratios from HMM "Hemi-preference ratio", a parameter inferred under our earlier analysis with a hidden-Markov model (HMM) [24], evaluates the preference of a given DNA methyltransferase for acting at hemimethylated as compared to unmethylated dyads [11], and thereby measures its preference for creating concordant dyads. Fu et al. [24,27] calculated this ratio for DNMT1, a mammalian maintenance methyltransferase. Because RCP measures the concordance preference of the entire ensemble of enzymes that give rise to methylation patterns, the hemi-preference ratio of a given enzyme and the RCP value inferred from the same data set are expected to have good agreement if that enzyme is the primary actor. The congruence between these two metrics is expected to decline with increasing contributions from other enzymes.
Three of the four data sets that had been analyzed under HMM showed very good agreement between the RCP values we infer here and the hemi-preference ratios previously inferred for the maintenance methyltransferase DNMT1: 58.0 vs. 58 for FMR1, 13.3 vs. 15 for G6PD, and 89.1 vs. 94 for LEP [24] (Table S1). The close correspondence between these values indicates that for these loci in leukocytes, methylation dynamics are driven primarily by conservative, maintenance-type processes such as accomplished by DNMT1, and that neither active demethylation nor de novo processes have a substantial role. Furthermore, when there is such a 24/33 correspondence, RCP strongly suggests that the mechanistic assumptions made for the enzymatic model hold for the given data set.
A large discrepancy, on the other hand, may suggest shortcomings of the mechanistic model. The fourth data set that had been analyzed under HMM, LEP in human adipose tissue, had an inferred RCP value of 34, well within the range of RCPs inferred for other data sets from single-copy loci including LEP in leukocytes (Fig 1a; Table S1). This RCP value was, however, more than eighteen-fold lower than the DNMT1 hemi-preference ratio estimate of 628 that we obtained under the earlier HMM approach (Table S1). What might account for this discrepancy? A hemi-preference ratio of 628 is unrealistically high, even for a maintenance enzyme, compared to the hemi-preference ratios inferred in other data sets, including those from methylation patterns established by DNMT1 in vitro [11]. This could reflect the inability of the HMM to yield a reasonable estimate for a data set impacted by demethylation, a process that was not considered in the HMM design. This lack of correspondence between HMM and RCP estimates is consistent with the possibility that active removal of methylation has a heightened role at loci with temporally variable transcription levels, and may reflect the role of the LEP locus as a sentinel of adipose but not blood [28].

S3 Text: Assessing the Potential Impact of Bisulfite-Conversion Errors
Two classes of error that occur during bisulfite-conversion can lead to misinterpretation of cytosine methylation states. Failure to convert unmethylated cytosines to uracil occurs at rate b, and can result in overestimation of methylation frequency. Inappropriate conversion of methylated cytosines to thymine, first noted by Shiraishi and Hayatsu, [49] occurs at rate c, and can result in underestimation of methylation frequency. We have reported previously that the rates of failed and inappropriate conversion depend strongly on the chemical and thermal conditions of bisulfite conversion. [48] In particular, we found that bisulfite treatment prolonged beyond that required to attain complete or nearly complete conversion of unmethylated cytosines can yield high rates of inappropriate conversion. Historically, conversion protocols have been designed to minimize the failed-conversion rate, with little or no attention to the rate of inappropriate conversion events. Thus, while both classes of error can alter parameters used to infer RCP -the methylation frequency, m, and the unmethylated dyad frequency, U -errors arising through inappropriate conversion are of particular concern.
How severely can conversion error impact RCP? And to what extent does its potential impact depend on the true methylation frequency of the target sequences? Densely methylated sequences contain a large number of fully methylated dyads, such that the most likely conversion error is inappropriate conversion yielding apparent hemimethylated dyads. Such events elevate the apparent level of discordance in a given data set and artifactually depress the inferred RCP values. To quantify this potential impact, we calculated how RCP would be altered by a single artifactual hemimethylated dyad introduced by inappropriate conversion, assuming that in its true state the relevant data set had one unmethylated dyad, one hemimethylated dyad, and 98 methylated dyads. The inappropriate conversion of one cytosine among 98 dyads whose true state is fully methylated is close to the number expected for an inappropriate conversion rate of about 0.5%. For this hypothetical data set, the introduction of a lone artefactual hemimethylated dyad by an inappropriate conversion event reduces RCP from about 20 to about 10, i.e., by a factor of 2. Thus, in the absence of mathematical corrections of the sort implemented here, even very low levels of inappropriate conversion can severely impact inference of RCP.
Sparsely methylated sequences contain a large number of unmethylated cytosines that are potential targets for failed conversion. We calculated that for such sequences the impact of 25/33 conversion error depends on whether most unmethylated cytosines are in hemimethylated or in unmethylated dyads. For example, when most unmethylated cytosines in a sparsely methylated sequence occur in hemimethylated dyads, the level of discordance is already high, such that the addition of an artifactual hemimethylated dyad by failed conversion only slightly decreases RCP. By contrast, when most unmethylated cytotosines are in unmethylated dyads, production of an artifactual hemimethyated dyad by failed conversion reduces RCP by a factor of two. Mathematical correction for conversion error is therefore critical, not only because error has potentially large impacts on absolute RCP values, but also because the magnitude of these impacts differs among data sets, with the potential either to magnify or to dampen variation among them.

S4 Text: Mathematical Correction for Bisulfite-Conversion Error
To account for conversion errors occurring at known rates, we express the observed dyad frequencies, M f obs , H f obs , and U f obs , as functions of the true dyad frequencies M t , H t , and U t , that would have been observed had conversion error not occurred: In cases where the mathematical correction yielded negative counts for one or two dyads, the negative counts were redistributed to the remaining dyads in proportion to their original counts, such that no dyad counts were negative after correction.
For each of the data sets reported here, DNAs were converted under one of two conditions: low molarity-temperature (LowMT), with failed-conversion rate 0.0030 and inappropriate conversion rate 0.031, and high molarity-temperature (HighMT), with failed-conversion rate 0.0086 and inappropriate-conversion rate 0.017 [48]. Error rates used in analysis of published data from other groups are given in Table S4 (Arand et al. [14,17]) and Table S5 (Zhao et al. [15]). Using these rates and observed dyad frequencies, true dyad frequencies can be inferred by maximum likelihood.

S5 Text: Assessment of Heterogeneity Among Replicates and Quasi-Replicates
To ask about possible heterogeneity among RCP values inferred for replicates -for example, for individual developmental stages as investigated by Arand et al. [17] -and for quasi-replicates -for example, for multiple cell lines at similar cell-culture and differentiation conditions -we utilized pairwise, two-tailed permutation tests.
These calculations revealed an intriguing pattern. There was no evidence of heterogeneity among samples of cultured human ES and iPS cells (p > 0.29, pairwise two-tailed PTs for the six cultured stem cell lines at the L1 locus; p-values summarized in Table S6; Fig 3). By contrast, several pairwise comparisons of mouse embryonic cells sampled independently for a given developmental stage yielded evidence of significant differences in RCP (p-values summarized in Table S6). Notably, however, most of the significant heterogeneity observed was for cells in pronuclear stages, at which data collection at identical stages is made difficult by rapid developmental transitions. On the other hand, there was evidence of only limited heterogeneity among totipotent and pluripotent cells between the 2-cell stage and the 3.5-dpc stage. Similarly, when experimental interventions were taken to prevent DNA replication or methylation (i.e. aphidicolin treatment at PN4/5 stages and SAMase treatment just after fertilization and before the first round of DNA replication), no evidence of heterogeneity was observed. Cell culture techniques are focused on minimizing opportunities for developmental differences to arise among cells, and could account for the observation of no heterogeneity among the cultured human ES and iPS cells.
Given the low levels of heterogeneity observed and the fact that some level of heterogeneity is inevitable in rapidly developing cells, we pooled several sets of replicates for analysis, as described in S6 Text (Fig 6).

S6 Text: Pooling Data Sets Across Replicates
In some instances, multiple independent replicates were collected and used to assess RCP for cells under a given biological condition -for example, multiple embryos at a given pronuclear stage (Fig 6). In these cases, the sequence-level dyad counts of individual replicate data sets were corrected for conversion errors independently before the replicate data sets were pooled.
When using likelihood methods to quantify and to assess pooled data sets -thus assuming independence of dyads (S8 Text) -we explicitly allow m to vary across replicates while estimating a single RCP value in the likelihood-maximization process. When using bootstrap methods (S7 Text), as we do whenever individual-molecule data are available, we do not explicitly account for the possibility that m may vary across the replicate data sets. This may result in a slight enlargement of confidence intervals in the bootstrap process.

S7 Text: Inferring and Comparing RCP Without Assuming Independent
Sampling of Dyads CpG dyads typically are not sampled individually, but instead as members of sequence reads that can contain from a few to many neighboring dyads. Often, there is correlation among the methylation states of these neighboring dyads. The processivity of the DNA methyltransferases, especially Dnmt1, is a substantial contributor to these correlations. As the mean number of dyads per read increases, so too does the potential for dyad-dyad correlation to undermine the accuracy of confidence intervals inferred under the assumption of independent sampling of dyads.
Of the three groups of data we analyzed here, one -Zhao et al. (2014) [15] -consists of reads that contain, on average, only two to three dyads, such that correlation among the dyads ascertained on a given molecule is anticipated to have a fairly minor impact on sampling. Reads for our data have as many as 22 dyads; mean dyad counts for data from Arand et al. [14,17] are intermediate between our data and those of Zhao et al. (2014) [15].
For our own data and those of Arand et al. [14,17], we used an approach that, at slightly higher computational cost, models and seeks to account for the potential impacts of dyad-dyad correlations. Our data yielded moderately larger confidence intervals under the bootstrapping approach as compared to under the likelihood approach with the assumption of independent sampling of dyads. By contrast, the data from Arand et al. [14,17] yielded almost identical confidence intervals without and with the assumption of independence. In view of the even smaller mean number of dyads per read in the data of Zhao et al. (2014) [15], we chose to make the not unreasonable assumption that dyads were sampled independently in their data. Although 27/33 it is possible that two or more reads originated from nearby regions of a single molecule and thus have dyad-dyad dependence, we assumed that the effect of such occurrences, if at all present, is very small, given the large amount of starting material.
The first of the two methods is described in this section, and the second in the next section, S8 Text.
Inferring RCP point estimates and confidence intervals. For our own data and those of Arand et al., we used a bootstrapping approach to model the uncertainty in RCP values introduced by possible within-sequence correlations in methylation states, and to make point estimate and confidence-interval inferences that account for this uncertainty.
For each data set of n sequences, we sampled n sequences with replacement, B = 2,000,000 times, and calculated the RCP for each of the B bootstrapped data sets. We inferred the true distribution of RCP for a given biological condition using the bootstrap distribution. It was clear from the resulting distributions that RCP is a biased estimator, as many of these distributions had longer right tails than left.
The simplest, and potentially misleading, approach for inference of point estimates and construction of confidence intervals from bootstrap distributions is to assume normality, and then to exclude right and left tails at the intended level of confidence.   [50] commented that, for biased estimators, this approach can lead to inference of inappropriately exclusive limits at the long-tailed end of the distribution, and inappropriately inclusive limits at the right-tailed side.
To address this problem, we applied Efron and Diccicio's "bias-corrected and accelerated" (BCa) method. Under the BCa method, the cumulative-density function observed for the distribution of bootstrap replicates is compared to that expected under normality. Bias-corrected point estimates are then inferred as the 50th-percentile values in the BCa corrected distributions. Similarly, critical points for intervals of a given confidence level are inferred from the values at the relevant percentiles of the distribution of bootstrapped values.
Because methylation states are predicted to be correlated across dyads within a molecule, but not across molecules, we resample at the level of molecules. Furthermore, as noted above, our molecular-barcoding procedures enable exclusion of redundant reads, such that each methylation pattern in our resulting data set is known to derive from a unique molecule in the original sample.
Moreover, it appears that molecules are sampled without bias due to methylation: in all eight data sets from murine DNA methyltransferase knockout lines (Fig 4) and all six data sets from wildtype human cells (Fig 3) there was no evidence of correlation between methylation frequency and RCP. Thus, it is reasonable to consider sampled molecules as independent and identically distributed draws from a population.
For each sampled molecule we derive a vector of values, (M i , H i , U i , n i ), where n i is the number of dyads. These are the vectors we resample in our procedure. We see that RCP can be written, using this vector, as whereM = niMi ni ,H = niHi ni andŪ = niUi ni . Now considering the vector (n i M i , n i H i , n i U i , n i ) and noting that the ratio of sums can be written as a ratio of means, we see that this falls directly under the "smooth function of means" framework introduced in [51] for the standard bootstrap and applied in [52] to the BCa bootstrap. From results in Section 6 of [52], it can be seen that our bootstrap procedure will give confidence intervals with 28/33 asympotically correct coverage.
Assessing whether a data set has RCP value greater than 1. If methyl groups are placed completely at random -that is, with preference for neither concordance or discordance -RCP is expected to be 1. In a previous report, Shipony et al. [16] interpreted their data to indicate that methyl groups are, indeed, placed essentially at random in undifferentiated cells. As there is very little evidence in any of our data sets to indicate possible preference for discordance, we opted to perform a one-tailed test, asking for the probability that our data sets do not have RCP values greater than 1.
To do so, we first calculated RCP point estimates for 200,000 bootstrap replicates and performed the BCa correction, using the method described above, and then calculated the p-value as the fraction of those point estimates that were less than or equal to RCP of 1. For example, the finding that only 20 of the 200,000 bootstrap replicates yielded an RCP point estimate less than 1, we would conclude that, for this data set, RCP is significantly greater than 1 (p = 0.0001).
Assessing whether RCP values differ significantly between data sets. A key goal of our study is to assess possible evidence for RCP differences between data sets. For example, we ask whether RCP for a given cell type differs between samples collected under differentiating as compared to non-differentiating conditions. When the two compared data sets can share a null distribution, as is the case for most comparisons we make here, we performed a permutation test (PT). For comparisons in which a shared null distribution cannot be established (such as would be the case if they were from different loci and thus had different numbers and locations of dyads), we used a bootstrap test (BT). We describe these two methods below.
To compute the significance of observed differences between RCP values for Data Sample A and Sample B, which can share a null distribution, we used permutations to compute a null distribution of ordered differences (for example, RCP(Sample A) − RCP(Sample B)) expected under the null assumption that the sequences in the two sets were drawn from a single population. To do this, we pooled sequences from Sample A and Sample B and then repeatedly drew from that pool, without replacement, to generate, at random, versions of Sample A and Sample B with sequence counts equal to those of the observed data. For each pair of randomly generated sets, we calculated the difference between their RCP point estimates (θ * ) , and obtained the distribution of RCP-difference values under the null hypothesis. For one-tailed tests, for example, to ask whether one data set has RCP value greater than the RCP value for another, we took the proportion of the distribution greater than the difference of the point estimates (θ) as the p value. For two-tailed/equal-tailed tests, for example, to ask whether RCP values for two data sets differ significantly, we first calculated the proportion of permuted differences smaller than (θ). We then calculated the proportion of permuted differences greater than (θ). We took the twice the lesser of these two values as the p-value for observing a difference this great in the event that the sequences in the two data sets were, in reality, drawn from the same distribution.
To compute the significance of observed differences between RCP values for two samples with different generating distributions, we used a bootstrap-based comparison test. Instead of pooling the data sets to form a single population of molecules, we bootstrapped a single RCP value from each of the two separate populations of molecules and computed the difference. We repeatedly sampled this difference to draw a bootstrap distribution of the difference in RCP values. For one-tailed tests, with which we examine directional differences, we determined the p value as the proportion of bootstrap-difference samples to the left of 0. For two-tailed tests, with which we can detect differences in any direction, we determined the p value as twice the smaller 29/33 proportion of the bootstrap difference samples on either side of 0.
Defining the 95% confidence region for a data set in the (m, U ) space, without assuming independent sampling of dyads. Using the methods described above, we generated 2,000,000 bootstrap samples of the data set. Instead of estimating the RCP value for each of the bootstrap samples, we calculated m and U . We constructed the two-dimensional confidence region for the two parameters for plotting using ci2d function in the gplots R package.

S8 Text: Inferring and Comparing RCP With Assuming Independent
Sampling of Dyads Calculating the likelihood of proposed true dyad frequencies, M t , H t , and U t , given the observed dyad counts, M c obs , H c obs , and U c obs . Failed-and inappropriate-conversion events create observed dyad frequencies that differ from true dyad frequencies. If the rates of these two types of error are known, the likelihood of a set of proposed true frequencies -M t , H t , and U t -can be calculated as follows, given the observed dyad counts -M c obs , H c obs , and U c obs -and Equation (4). M f obs , H f obs , and U f obs , which indicate the observed frequencies of the dyads, can be easily calculated from the dyad counts.
Inferring RCP point estimates and confidence intervals. The RCP point estimate of a data set is calculated directly from the conversion-error-corrected observed dyad frequencies.
We determine the 95% confidence interval for RCP as the interval that includes all values of RCP for which the natural log likelihood lies within χ 2 0.95; df =1 = 1.92 units of the maximum natural log likelihood point estimate. [53] Although bias is just as much a concern with the assumption of independent sampling of dyads, we did not perform a bias correction for the data from Zhao et al. [15], because without bootstrapping, we lacked a simple method for bias estimation. Nonetheless, our observations from our analyses of other data sets suggest that most of these samples are likely not severely affected by bias. From bootstrapping of smaller data sets collected by our lab and by Arand et al. [14,17], we have observed that the asymmetry in the distribution is small in the lower ranges of RCP (1∼10), but greater as RCP increases. Therefore, bias is likely to be small for most samples presented by Zhao et al., for which RCP point estimates rarely exceed 10. For data sets presented by Zhao et al. that yielded RCP estimates greater than 10, bias are unlikely to affect our analyses and the general conclusion that their methylation behavior is characterized by strong preference for concordance.
Assessing whether RCP values differ significantly between data sets. To compare the RCP values between two data sets while assuming independence among all dyads, we can use a likelihood approach, comparing a model in which two true RCP values are required to described the two data sets to an alternate model in which both data sets can be explained by a single RCP value. We implemented that test as follows: Under the assumption of large sample of dyads, D is approximately χ 2 distributed with 1 degree of freedom.
Assessing whether a data set has RCP value greater than 1. We again take a likelihood-based approach as we did in the section above. Here, the null model states that the system operates under the specified RCP value. The alternate model states that the system operates under another RCP value.
Under the assumption of a large sample of dyads, D is approximately χ 2 distributed with 1 degree of freedom.
Defining the 95% confidence region for a data set in the (m, U ) space. We determine the 95% confidence region in the space of two parameters -here m and U -as the region that includes all proposed pairs of parameter values for which the natural log likelihood lies within χ 2 0.95; df =2 = 3.00 units of the maximum natural log likelihood point estimate [53].

S9 Text: Could Spontaneous Differentiation of a Subset of ES and iPS
Cells Substantially Influence the Inference of RCP?
One possible explanation for the inference of conservative processes operating in cultured, undifferentiated cells is that these cells may in reality be a mixture of differentiated and undifferentiated cells. We calculated how large the differentiated subpopulation would need to be under this scenario, given that the subpopulation operated with the RCP inferred for the corresponding differentiated cells. The remainder of the population was assumed to operate at RCP of 1, per the alternate hypothesis. We allowed specification of m for each of the two populations.
Let p 1 and p 2 represent the proportions of the two putative subpopulations that we wish to estimate, such that p 1 + p 2 = 1. We start with RCP and m for each of the two subpopulations, which we denote by RCP 1 , RCP 2 , m 1 , and m 2 . We can find U 1 and U 2 using the RCP and m values. We have the following set of equations: p 1 m 1 + p 2 m 2 = m overall p 1 U 1 + p 2 U 2 = U overall Using these equations and the observed overall RCP value, and the expression for RCP given m and U (Equation (3)), we can find p 1 and p 2 .
Assuming that m of the undifferentiated subpopulation is at least 0.2, we found that over half of the cells would need to be differentiated if methylation in the remaining undifferentiated cells operates at RCP of 1. Morphological inspection of the cultured human stem cells did not suggest the presence of a substantial subpopulation of differentiated cells. Additionally, this calculation is not critical for the work presented here, because RCP values for totipotent and pluripotent stem cells from murine embryos are very similar to those inferred for all cultured ES and iPS cell lines that we examined. Such calculations may, however, be useful for future evaluation of cultured cells that are early in the process of differentiation. 32/33 S10 Text: Comparing RCPs of Dnmt3-Knockout Lines with Those of Wildtype Lines The relative contributions of conservative processes are expected to be more substantial when the fraction of methylation events achieved through maintenance-type activity is elevated through loss of one or both of the de novo enzymes. Indeed, for the Lep locus, we inferred RCP values of 8.11 (Dmnt3a KO) and 6.97 (Dnmt3b KO), both higher than the RCP value for the wildtype murine ES cell line (p = 0.008, Dnmt3a KO; p = 0.04, Dnmt3b KO) (Fig S2a;  Table S3).
Indeed, at the Lep locus, we inferred RCP values of 8.11 for Dmnt3a KO cells, a value significantly higher than that for wildtype cells (p = 0.046, two-tailed PT; Fig S2a; Table S3). Although the point estimate was higher for Dmnt3b KO cells as well, the difference was not significant (p = 0.11, two-tailed PT).
Examination of data from Arand et al. [14] for a broader set of loci in cell lines with knockouts at one or both of the Dnmt3s, however, revealed a more complex role for these enzymes in shaping methylation concordance. Overall methylation levels were diminished in the absence of the de novo methyltransferases 3a and 3b, as expected (Fig S2b- Table S3); RCP values increased as expected for several, though not all loci. Results for cells that lacked only one of the two DNA methyltransferases were even more variable across loci ( Fig S2, Table S3). In some cases, loss of a single Dnmt3 enzyme had the predicted impact of increasing RCP (L1 in 3a KO, p < 10 −16 , Fig S2h; L1 in 3b KO, p = 0.002, Fig S2h; IAP in 3a KO, p = 0.002, Fig S2g; IAP in 3b KO, p = 0.029, Fig S2g; two-tailed PTs), as was observed for Lep in our data. For one locus, B1, knockout of both of the Dnmt3s resulted in significantly increased RCP (p < 10 −16 , two-tailed PT), though neither of the single knockouts did (Fig S2f). Most surprisingly, for two loci, loss of either of the Dnmt3s decreased RCP (Igf2, p < 0.016, Fig S2c; Igf2, p < 10 −16 , Fig S2c; two-tailed PTs). These counterintuitive results likely reflect variation across loci in the roles of the individual DNA methlytransferases -and possibly the demethylation machinery -in shaping overall methylation levels across loci and categories of genic elements [54].