Comparing DNA replication programs reveals large timing shifts at centromeres of endocycling cells in maize roots

Summary

Plant cells undergo two types of cell cycles–the mitotic cycle in which DNA replication is coupled to mitosis, and the endocycle in which DNA replication occurs in the absence of cell division. To investigate DNA replication programs in these two types of cell cycles, we pulse labeled intact root tips of maize (Zea mays) with 5-ethynyl-2’-deoxyuridine (EdU) and used flow sorting of nuclei to examine DNA replication timing (RT) during the transition from a mitotic cycle to an endocycle. Comparison of the sequence-based RT profiles showed that most regions of the maize genome replicate at the same time during S phase in mitotic and endocycling cells, despite the need to replicate twice as much DNA in the endocycle and the fact that endocycling is typically associated with cell differentiation. However, regions collectively corresponding to 2% of the genome displayed significant changes in timing between the two types of cell cycles. The majority of these regions are small with a median size of 135 kb, shift to a later RT in the endocycle, and are enriched for genes expressed in the root tip. We found larger regions that shifted RT in centromeres of seven of the ten maize chromosomes. These regions covered the majority of the previously defined functional centromere, which ranged between 1 and 2 Mb in size in the reference genome. They replicate mainly during mid S phase in mitotic cells but primarily in late S phase of the endocycle. In contrast, the immediately adjacent pericentromere sequences are primarily late replicating in both cell cycles. Analysis of CENH3 enrichment levels in 8C vs 2C nuclei suggested that there is only a partial replacement of CENH3 nucleosomes after endocycle replication is complete. The shift to later replication of centromeres and possible reduction in CENH3 enrichment after endocycle replication is consistent with a hypothesis that centromeres are inactivated when their function is no longer needed.

Keywords:

Gene expression – Plant genomics – Maize – DNA replication – Cell cycle and cell division – Centromeres – Synthesis phase – Genomic signal processing

Introduction

Developmentally programmed DNA replication without nuclear breakdown, chromosome condensation or cell division, a phenomenon known as endoreduplication or endocycling, occurs in a wide variety of plants and animals [1–3]. In plants, endoreduplication is a systemic feature [4] that is often an important step in the development of tissues and organs such as fruit, endosperm, leaf epidermal cells and trichomes [5]. Initiation of endocycling is frequently associated with a transition from cell proliferation to cell differentiation and expansion [6]. In plant roots, cells at the tip divide actively by normal mitosis, while endocycling cells become frequent further from the tip, in a zone associated with differentiation and increased cell size [7, 8]. In the maize (Zea mays) root, the first 1-mm primarily contains actively dividing mitotic cells [8, 9]. Mitotic activity continues in the region between 1 to 3-mm from the tip, but in this region about 30% of the cells undergo a single endocycle. Instead of undergoing mitosis and returning to a 2C nuclear DNA content, these cells replicate their DNA without undergoing mitosis and transition directly from 4C to a final DNA content of 8C [8]. Beyond 3 mm, cell differentiation predominates, and DNA replication activity becomes rare [8, 10]. Thus, while not all differentiating cells undergo an endocycle, in those cells that do, the endocycle precedes or accompanies cellular differentiation.

In animal systems, different cell types often exhibit differences in the temporal order of DNA replication along the chromosomes, also called the “replication timing (RT) program”. These differences in RT programs are likely related to changes in chromatin and gene expression during the differentiation process (e.g. [11–14]). However, as yet there is no information concerning possible changes in RT programs associated with endoreduplication or differentiation in plant systems. Other than in a few model systems (e.g. [15, 16]), it is difficult to separate individual cell types from plant organs, and DNA replication occurs mainly in meristematic regions, which are small and often difficult to access. In addition, plant cells do not maintain their differentiated state when grown in suspension culture [17, 18]. Hence, there are few opportunities to compare RT programs in individual cell types. However, the occurrence of endocycling and differentiating cells near the apical meristem of maize root tips, in which we previously characterized the mitotic RT program [19], offers a unique opportunity to compare the two modes of replication in an intact plant organ.

We developed a system to analyze DNA replication in maize roots [8, 20], using similar approaches to those being applied in our work with Arabidopsis cell suspensions [21]. In this system, newly replicated DNA is labeled in vivo with the thymidine analog, 5-ethynyl-2’-deoxyuridine (EdU), and labeled nuclei are separated by flow cytometry into populations representing different stages of S phase. One key advantage of using EdU as the nucleotide analog instead of 5-bromo-2’deoxyuridine (BrdU), which has traditionally been used in replication timing by sequencing (Repli-seq) protocols [12], is the ability to sort nuclei on both DNA content as well as EdU content, followed by immunoprecipitation of EdU labeled DNA. This modification allows for clean separation during flow sorting of non-replicating, unlabeled G1 and G2 nuclei from the labeled S-phase nuclei.

Cytological analysis of sorted, EdU-labeled nuclei showed that replication activities in early and mid S are more closely interspersed in the maize nucleus than in animal cells [22, 23]. We characterized the RT program in mitotic cells of the apical 1-mm root segment [19], using the modified Repli-seq protocol [24]. In mitotic cells, we found evidence for a gradient of early replicating, open chromatin that transitions gradually into less open and less transcriptionally active chromatin replicating in mid S phase. We also confirmed cytological observations showing that heavily compacted classical heterochromatin, including knobs and pericentromeres, replicate primarily in late S phase [22, 25]. While these relationships between RT and chromatin packaging are generally similar to those found in other systems, we did not find evidence for megabase-scale replication domains like those that have been characterized in mammalian cells (reviewed in [26] and references therein).

Although replication in the first 1-mm of the root is mostly mitotic, with DNA contents of labeled nuclei ranging from 2C to 4C, flow cytometry profiles of nuclei derived from root tissue between 1 and 3-mm from the tip also included a substantial population of nuclei with DNA contents beyond 4C. Appearance of endopolyploid (8C) nuclei in this zone would be expected, as some 4C nuclei are known to enter the endocycle rather than undergo division [7, 8]. Cytological analysis of replicative labeling showed that the spatiotemporal patterns of replication during the 4C to 8C transition in these endocycling nuclei are very similar to those in mitotic nuclei [22]. However, it remained to be determined whether the entire genome is uniformly replicated during the endocycle, and whether the temporal program is altered in differentiating cells when replication occurs without an intervening mitosis.

Both under-replication and over-replication (amplification) have been observed in multiple animal systems during developmentally programmed endocycles, notably including Drosophila (reviewed in [27]). In addition to the well-known amplification of chorion genes and under-replication of heterochromatin, under-replication also occurs in a number of euchromatic regions, with a degree of tissue specificity suggesting a possible role in differentiation [28–30].

Even though endopolyploidy is common in plants, there are very few reports dealing with over- or under-replication of specific sequences. Some orchids exhibit a phenomenon in which only a fraction of the genome is endoreplicated [31, 32], but in most cases, endopolyploid cells have DNA contents that are multiples of the 2C value. Both highly repetitive heterochromatic regions and highly expressed genes are extensively endoreduplicated in maize endosperm nuclei, as would be expected for uniform replication of the entire genome [33]. More definitively, whole genome sequencing in Arabidopsis showed that leaf nuclear DNA is evenly endoreduplicated in wild-type plants, although the same series of experiments clearly demonstrated selective over-replication in atxr5 and atxr6 mutants [34].

To address the question of over- or under-replication and whether there are developmentally associated changes in RT programs in the maize root tip system, we performed a detailed comparison of RT dynamics in mitotic and endocycling cells. We found very little evidence for over- or under-replication in endocycling cells, consistent with the few previous reports on this topic from plant systems. We also found that the RT programs for the vast majority of the genome are very similar. However, we found significant changes in timing for a number of loci that together correspond to 2% of the genome. Most notably, we found major changes in the RT of centromeres, which replicate mainly during mid S phase in mitotic cells but primarily in late S phase of the endocycle.

Results

Separating endocycling from mitotic nuclei

As reported previously and described in Methods, we used a 20-min pulse of the thymidine analog, EdU, to label newly replicated DNA in intact maize roots. This was followed by formaldehyde fixation and isolation of nuclei from defined segments of root tips (Fig 1A). Incorporated EdU was conjugated with Alexa Fluor 488 (AF-488) by “click” chemistry [35]. The nuclei were then stained with DAPI and fractionated by two-color fluorescence activated flow sorting to generate populations at different stages of the mitotic cell cycle or the endocycle [8, 20]. Fig 1B and 1C show flow cytometry profiles obtained for root segments 0–1 mm and 1–3 mm from the tip, respectively. Fluorescent signals from nuclei that incorporated EdU during S phase of a normal mitosis form an “arc” between 2C and 4C DNA contents (Fig 1B). Previous EdU pulse-chase time course experiments in the 0–1 mm region showed that most labeled 4C nuclei return to 2C as would be expected in a mitotic cycle (see Fig 1B in [36]). In the 1–3 mm zone, a substantial fraction of labeled nuclei is undergoing endoreplication, forming a similar arc with DNA contents ranging between 4C and 8C (Fig 1C). As seen in Fig 1C, the endocycle arc is more prominent in nuclei preparations from 1–3 mm root segments. To analyze the endocycle RT program, which is described in detail below, labeled nuclei representing early, mid, and late S-phase fractions were separated using the sorting gates shown in Fig 1C, adjusting the endocycle early gate to avoid contamination with mitotic nuclei in late S phase. Reanalysis of the sorted nuclei confirmed that there was good separation between the nuclei populations from the adjusted early sorting gate and the mid sorting gate (S1 Fig). To complete the Repli-seq protocol, described in more detail below, the DNA was extracted from the nuclei in each gate and sheared. The labeled, newly replicated DNA from each S-phase fraction was then immunoprecipitated and sequenced.

The flexibility of the EdU labeling and flow sorting system also allowed us to collect unlabeled nuclei, representing non S-phase cells with 2C, 4C and 8C DNA contents. These nuclei were used to characterize selected histone marks following mitotic or endocycle replication and to investigate the copy number of individual loci across the genome.

Evidence for complete genome replication during the endocycle

Given the well documented examples of over- and under-replication during the endocycle in animal systems, we investigated whether there are local copy number differences in the maize genome after endocycle replication. To do this, we used the non S-phase 2C, 4C, and 8C nuclei populations described above, and carried out whole genome paired-end sequencing. To gain a better representation of the copy number of repeat regions in the genome, reads that could not be uniquely mapped to a single location were included, but we retained only the primary alignment location for each read pair. This approach was used exclusively for the copy number analysis, while all subsequent analyses included only uniquely mapping reads (see Methods and S1 Text). The data for each ploidy level were examined for regions in which normalized read frequencies in 5-kb windows differed between 8C and 4C or 4C and 2C nuclei, using procedures described by Yarosh et al. ([37]; S1 Text). We found about 5% of the 5-kb windows had ratio values that fell outside of two standard deviations of the mean ratio for 4C and 2C or 8C and 4C (1.0 ± 0.2 S. D. for both; S2A and S2B Fig). However, these windows all either occurred as singleton 5-kb windows scattered around the genome (S2C Fig) or coincided with regions that had very low read mapping in the 2C sample, indicating they are likely the spurious result of making a ratio between windows with very few reads in both samples. As such, there is very little evidence of meaningful over- or under-replication of genomic regions in nuclei with different ploidy levels.

To further investigate whether there is complete replication of high-copy repeats that are not well represented in the genome assembly, we used BLAST software to query all reads, not just those that can be mapped to the genome, to determine the percentage of reads corresponding to each of several consensus sequences for high-copy repeats (S1 Text). Analyzed sequences included the knob repeats knob180 and TR-1 [38, 39], 5S and 45S rDNA repeats [40], and centromere-associated CentC satellite repeats [41]. We also queried consensus sequences for centromere retrotransposons of maize (CRM) families 1–4 [42–45]. In all cases, we found the percentages to be similar in the 2C, 4C and 8C samples (S2D and S2E Fig), further suggesting that there is little or no over- or under-replication.

Replication timing analysis

As described above, we sorted endocycling nuclei from the S-phase populations in Fig 1C, and extracted and sheared the DNA in each fraction. EdU-containing DNA fragments were immunoprecipitated (IP) with an antibody to AF-488, resulting in sequence populations representing DNA replicating during early, middle, or late S phase of the endocycle. Given that EdU/AF-488-labeled DNA is immunoprecipitated from a population of only EdU-labeled nuclei (see Fig 1B and 1C), the level of unlabeled DNA background in each IP is substantially lower than would be present if nuclei were sorted simply by DNA content. We also prepared DNA from the unlabeled 2C nuclei pool to provide a reference dataset representing pre-replicative nuclei. DNA from three biological replicates of each sample was sequenced to generate paired-end reads.

To compare the RT programs in endocycling and mitotic nuclei, we mapped our previous Repli-seq data for mitotic nuclei [19] and our new data for endocycling nuclei to the maize B73 RefGen_v4 genome, which includes improved assemblies of centromeres and more complete annotations of transposable elements (TEs) [46, 47]. Including only uniquely mapped reads resulted in a read depth that varied between 65.7 million and 261.2 million reads per S-phase fraction (including reads from 3 biological replicates per S-phase fraction). The data from all S-phase samples for the mitotic cycle and endocycle were then randomly downsampled to the lowest read depth (65.7 million reads) to ensure comparable results (see Methods and S1 Spreadsheet).

We used the Repliscan analysis pipeline [24] to generate continuous normalized data profiles representing the intensity of replication activity across the genome (RT profiles) in early, mid and late fractions of each S phase. These RT profiles were generated by aggregating the Repli-seq read densities for each S-phase sample in 3-kb static windows, scaling the reads to 1× genome coverage, and then dividing by the scaled read counts from the unlabeled 2C reference data and smoothing by Haar wavelet transform (see Methods and [24]). Normalizing with the pre-replicative 2C reference provided a uniform 2C copy number and corrected for differences in sequence mappability and collapsed repeats that caused “spikes” in the data (illustrated for late replication in the endocycle in S3 Fig) while preserving replication signal [24, 48]. The 3-kb windows identified as having no or extremely low read coverage in the 2C reference sample (see Methods) were excluded from all analyses. These windows include both unmappable and multi-mapping regions (“blacklist” windows, indicated by black tick marks in Fig 1E). After these RT profile normalization steps, the result is an estimate of the intensity of replication activity in each 3-kb window, which we refer to as “replication signal”.

Fig 1D shows that the global RT profiles are remarkably similar in endocycling and mitotic nuclei, and overlays of the corresponding profiles show mostly minor differences (Fig 1E). Pearson’s correlation coefficient values between corresponding S-phase fractions from the mitotic and endocycle data are very high (r values of 0.91, 0.89 and 0.96 for early, mid and late, respectively). These values are similar to those found between individual biological replicates within each sample (S4 Fig).

Identifying regions of altered timing

Despite the global similarity of the RT programs of mitotic and endocycling cells, there are regions scattered around the maize genome that show a shift in RT. To identify regions with differences in RT (DRT), we first calculated the difference in normalized replication signal between the mitotic and endocycle data at each genomic location for the early, mid and late profiles separately (S1 Table; S5 Fig). We then constrained our analysis by focusing only on regions where there was an equal and opposite DRT in at least one other S-phase fraction (for example, regions in which a decrease in early replication signal in endocycling cells was associated with a corresponding increase in mid and/or late replication signal at the same location). We allowed a gap distance of 6 kb when searching for regions with DRT to account for small blacklist regions that break up larger regions of change. We found that 11% of the genome showed a difference in replication signal of at least 10% of the total difference range for a given S-phase fraction (absolute difference in replication signal ≥ 0.4; S1 Table), with an opposite difference in replication signal at the same threshold criterion at the identical location in another S-phase fraction. Many of these regions are small, with the lower 50% of regions ranging in size from 3 kb to the median size of 33 kb (S2 Table). Since the units of replication initiation, elongation and termination (“replicons”) in monocot plants are estimated to be on average 47 ± 13 kb in size [49], the biological relevance of RT differences across regions much smaller than this is not clear. However, instead of implementing an arbitrary size cutoff, we chose to focus on regions that were associated with at least one core region with a larger (≥ 25%) DRT, as described in the next paragraph.

To identify these more robust RT differences, designated Regions of Altered Timing (RATs), we identified regions in which the DRT was ≥ 25% of the total difference range for a given S-phase fraction (absolute difference in replication signal ≥ 1.0; S1 Table), and which also met the criterion of having an opposite DRT in at least one other S-phase fraction. To highlight larger and contiguous regions of change, we included ≥ 10% regions that were adjacent to the original ≥ 25% regions. However, RATs had to have at least one core region where the DRT was at least 25% (S2 Table) to be included in our analysis. Representative ≥ 25% and ≥ 10% regions are indicated by various shades of red and blue bars in Fig 2 (additional examples are in S6 Fig). Finally, we examined the RT profiles for the RATs in individual biological replicates to verify there was good agreement between the replicates (Figs 2B and S6). By selecting only the most robust RATs we excluded other regions where RT changes are less dramatic–for example those indicated by dashed boxes in Fig 2. In such regions, the DRT did not meet our criteria of a ≥ 25% difference in replication signal (box 2 in Fig 2A) and/or there is not an equal and opposite (“compensated”) DRT (box 3 in Fig 2A).

Robust RATs fall into two categories, those in which the strongest replication signal occurs later in the mitotic cycle than it does in the endocycle (“Later-to-Earlier” shift), and those in which the strongest replication signal occurs earlier in the mitotic cycle than in the endocycle (“Earlier-to-Later” shift). In addition, we separately characterized a subset of the Earlier-to-Later RATs that are located in functional centromeres (“Earlier-to-Later-CEN”) using centromere (CEN) coordinates from [46]. Our stringent criteria identified RATs comprising only about 2% of the maize genome (Table 1), with 233 of the 274 total regions (representing 1.6% of the genome) in the Earlier-to-Later category. Non-CEN Later-to-Earlier and Earlier-to-Later RATs have similar size distributions, with median sizes of 141 and 135 kb, respectively (Fig 2C and Table 1). All of the CEN RATs fall into the Earlier-to-Later category and have a median size of 132 kb, similar to the non-CEN RATs. It is important to note, however, that the sizes of CEN RATs are underestimated because of numerous blacklist regions within the centromeres that break what are likely continuous RATs into several smaller parts in our analysis. Even though maize centromeres are remarkably well sequenced [46], they still contain some gaps and regions where reads cannot be uniquely mapped in the current B73 RefGen_v4 genome assembly. To account for this, we calculated an upper estimate of continuous RAT size by incorporating the coverage from the multiple RATs called within each centromere as well as coverage from RATs that extend past the previously determined CEN boundaries (“presumed CEN RATs”), and including the interspersed blacklist regions. This upper estimate of CEN RAT size, which is addressed in more detail below, ranges from about 0.9–2 Mb in each CEN (Table 2).

Non-centromeric RATs

We analyzed the non-CEN RATs for the content of genes and TEs, as well as the presence of histone modifications and functional annotations related to the genes within RATs. To assess whether the percentage of RATs containing genes differed from random expectation, we randomly shuffled coordinates corresponding to the non-CEN Later-to-Earlier and Earlier-to-Later RATs around the genome 1000 times and calculated the percentage of randomly shuffled regions that overlap genes in each of the 1000 sets. These “expected” random distributions were compared to the observed percent overlap values found for Later-to-Earlier and Earlier-to-Later RATs and a permutation P value was calculated (see Methods). We found that 93% and 96% of Later-to-Earlier and Earlier-to-Later RATs, respectively, contain at least one annotated gene and usually contain a small cluster of, on average, 2–3 genes (Fig 3 and S3 Table). The observed 96% percent overlap of Earlier-to-Later RATs with genes is significantly greater than expected by chance (permutation P value = 0.001; Fig 3). To assess which of these genes are expressed, we used root-tip RNA-seq data that are not specific to mitotic or endocycle cells, and found that although only 50% of the 682 genes found in non-CEN RATs are expressed at a meaningful level (FPKM ≥ 1), 83% and 91% of Later-to-Earlier and Earlier-to-Later RATs, respectively, still contain at least one expressed gene (S3 Table). The observed 91% overlap of Earlier-to-Later RATs with expressed genes is also significantly greater than expected by chance (permutation P value = 0.001; Fig 3). In contrast, we found only 51% and 68% of Later-to-Earlier and Earlier-to-Later RATs, respectively, contain genes not expressed in the root (Fig 3). These values for non-expressed genes are not significantly different from random expectation, indicating that the enrichment of genes in Earlier-to-Later RATs is mainly driven by the expressed genes.

The significant enrichment of genes expressed in the root in Earlier-to-Later non-CEN RATs suggests the possibility that these regions may be related to shifts in gene expression. However, we were unable to directly compare expression of genes in RATs in mitotic and endocycling cells because we could not obtain RNA of sufficient quality to sequence from fixed, sorted nuclei. Instead, we assessed a selection of gene-associated histone post-translational modifications in sorted non S-phase 2C, 4C and 8C nuclei. In our previous work in maize root mitotic cells, we showed that trimethylation of H3 lysine 4 (H3K4me3) and acetylation of H3 lysine 56 (H3K56ac) modifications tend to colocalize on active genes and are associated with earlier replicating regions, while trimethylation of H3 lysine 27 (H3K27me3) tends to be on repressed genes regardless of their RT [19]. For each ploidy level, we quantified the percentage of genes within RATs that have each mark, as well as the fold enrichment relative to input for called peaks within genes. There are very few differences between ploidy levels in the number of genes bearing these marks (S7D Fig), but there are some minor shifts in the peak enrichment in 8C nuclei compared to 2C (S7A–S7C Fig). The clearest shift is a decrease in H3K4me3 enrichment found on expressed genes in Earlier-to-Later RATs (S7B Fig), which suggests these genes may have decreased expression in endocycling cells.

We also performed a gene ontology (GO) analysis for the genes found in non-CEN RATs to ask if there are functional annotations enriched in genes that shift RT. For this analysis, we focused on the genes that we identified as expressed in the root tip (S2 Spreadsheet). We found 44 significantly enriched GO terms for genes within Earlier-to-Later RATs, including biological process and molecular function terms related to gene expression, DNA/RNA metabolism, and the cell cycle (S8 Fig). A wide variety of significant cellular component GO terms were also found, which may relate to various differentiation processes occurring in endocycling cells. There are no significant GO terms for genes within Later-to-Earlier RATs, though the presence of only 52 expressed genes in this RAT category made it difficult to fully assess significance. Taken together, these analyses of transcription-related histone modifications and functional annotations suggest a role for gene expression changes in the Earlier-to-Later RATs. Given that these regions are shifting to a later RT in the endocycle, a decrease in gene expression would be expected [19]. However, more work will be needed to confirm this hypothesis.

The general organization of the maize genome is genes clustered in “islands” interspersed with blocks of transposable elements [50–52]. We used a permutation strategy similar to that described above to estimate the significance of any differences in percent coverage of TEs and individual TE superfamilies in non-CEN RATs. Observed coverage values were compared to random expectation, estimated from 1000 randomly shuffled sets. The TE annotations were from the recent B73 RefGen_v4 TEv2 disjoined annotation, where every bp is assigned to a single TE [47]. We found no TE superfamilies with percent coverage values in non-CEN RATs that are significantly different from random expectation at the permutation P value threshold of 0.001, although one TE superfamily, RLG/Gypsy, is very close to the threshold for being called significantly depleted in Earlier-to-Later RATs (permutation P value = 0.002; S9 Fig). RLG/Gypsy elements make up 32% of the coverage of Earlier-to-Later RATs, and from 32–40% of the coverage of the randomly shuffled sets. We also found that the percent AT content in RATs is similar to that of the genome as a whole, with median values of 55% and 56% for Later-to-Earlier and Earlier-to-Later RATs, respectively, compared to a median value of 55% for the whole genome (S10 Fig). Taken together, there is no evidence of a difference in AT content or a major enrichment or depletion in the coverage of specific TE superfamilies in non-CEN RATs, with the possible exception of a minor depletion of RLG/Gypsy element coverage in Earlier-to-Later RATs. Given that RLG/Gypsy elements have by far the most abundant coverage across the B73 reference genome [46], this minor depletion may be related to the increased presence of genes in Earlier-to-Later RATs (Fig 3).

Centromeric RATs

Functional centromeres are defined by their content of nucleosomes containing the centromere-specific histone variant known as CENH3 in plants and CENP-A in animals. CENH3/CENP-A makes up only a small percentage of the total H3 population in centromeres, but plays an important role in recruiting kinetochore proteins [53–55]. Maize is unusual among higher eukaryotes in that a majority of centromeric reads can be uniquely mapped [56]. In our Repli-seq data, for example, we found that on average 45% of all reads that map to centromeres could be uniquely mapped to a single location (S11 Fig). Only these uniquely mapping reads were used for further analysis. In addition, most of the maize centromere assemblies are relatively intact, and functional centromeres have been located by mapping ChIP-seq reads for CENH3 [46]. When combined with our replication timing data, these features of the maize system create a unique opportunity to assess RT programs for centromeres.

Our analysis found large, robust RATs across seven of the ten centromeres (Figs 4C, 4D and S12). These seven centromeres (CEN 2, 3, 4, 5, 8, 9 and 10) were previously classified as “complex” because they contain a mixture of retrotransposons with some centromere satellite repeat arrays (CentC; [56, 57]). In the RefGen_v4 genome assembly, CEN 9 has two called CENH3-binding regions [46], which we refer to as CEN 9a and 9b (Fig 4C; grey bars). Interestingly, we only found robust RATs in the larger CEN 9a, with the smaller CEN 9b showing very little RT shift (Fig 4C and Table 2). The cumulative coverage of RATs in each complex CEN ranges from 0.45–1.5 Mb (Table 2). However, because each centromere includes blacklist regions that vary in size and number, automated analysis did not identify the true sizes of the RATs. Therefore, we also calculated an upper estimate of continuous RAT size in each centromere, including the multiple RATs called by automated analysis as well as the interspersed blacklist regions. We also included coverage from CEN RATs that extend somewhat past the previously determined CEN boundaries (“presumed CEN RAT”). This upper estimate ranges from about 0.9–2 Mb in different centromeres, and covers between 65–95% of the previously defined functional centromere for the seven complex CENs (Table 2).

Upon inspecting the RT profiles for each of the complex CENs we found the strongest replication signals occurring mainly in mid S in mitotic cells, but changing to primarily late S in endocycling cells. It is also noteworthy that though replication occurs mainly in mid S in mitotic cells, there are some distinct peaks of early replication inside or directly adjacent to the called centromere (indicated by black arrowheads in Figs 4 and S12) in all but one (CEN 5) of the complex centromeres. These early peaks remain in the endocycle, though usually there is some reduction in early signal with a concomitant increase in mid signal at the same location.

The remaining three centromeres (CEN 1, 6, and 7) were previously characterized as “simple” because they mainly contain large arrays of the CentC repeat [56, 57]. In our analysis, the simple centromeres showed, at most, small RT shifts that did not meet our criteria for a robust RAT (Table 2, Figs 4A, 4B and S12). However, CentC repeats are not well represented in the reference genome assembly, so it is not possible to analyze RT profiles for the complete simple centromeres. Portions of CEN 7 that are present in the assembly replicate mainly in mid S phase in both mitotic and endocycling cells (S12 Fig), while sequences in the assemblies for CEN 1 and CEN 6 are mostly late replicating in both types of cells, with some minor RT changes across small regions (Fig 4A and 4B and Table 2).

Because of the issues with computationally identifying continuous RATs we chose to focus the following set of analyses on the entire CENH3-binding region of each chromosome (excluding blacklist regions). We calculated the difference in early, mid and late replication signal (endocycle minus mitotic) from RT profiles by averaging across 100-kb static windows. For comparison, we also calculated the replication signal differences in pericentromeres, which were arbitrarily defined as the ± 1 Mb flanking the CENH3 region. We inspected all RT differences in the centromeres and pericentromeres by not requiring that the DRT be compensated by an opposite shift in the other S-phase fractions. Early and mid replication signals across the complex centromeres decrease and late replication signals increase in endocycling cells, reflecting a large shift toward late replication. The DRT values for the late profile in centromeres and pericentromeres are shown in Fig 4E and 4F, respectively, while the DRT values for early and mid profiles are shown in S13 Fig. Interestingly, the DRT tapers off towards the edges of the functional centromere (see profiles in Figs 4C, 4D and S12), and there is striking congruity in the replication signals for mitotic and endocycling cells in the immediately adjacent pericentromere regions (Fig 4A–4D). The few RT shifts in pericentromeric regions are smaller in size and much less dramatic than those in the centromere proper (Fig 4F). Moreover, very few (8%) of pericentromeric windows with DRT are compensated by an equal and opposite shift in the other S-phase profiles (S4 Table), suggesting many of these uncompensated differences may result from technical variation rather than from meaningful biological differences. In contrast, nearly all (85%) of the centromeric windows have compensated RT shifts.

Genomic elements and features in centromeres

Maize centromeres contain varying amounts of tandemly arrayed CentC repeats (single repeats of 156 bp in length; [41]) as well as several CRM retrotransposon families interspersed with elements from a few other retrotransposon families [44, 52, 58, 59]. CentC repeats and CRM elements are also present in the adjacent pericentromeres where there is no CENH3 binding [52, 58]. In RefGen_v4, there are also fifty annotated genes within centromeres. We asked if all of these sequence elements in centromeres behave similarly in the mitotic to endocycle transition, or if certain elements show larger RT shifts than others. We also asked if all three types of sequence elements show similar RT changes in centromeres versus pericentromeres. Given that the replication signal values were aggregated in 3-kb windows, we only included elements that covered at least half a window (1.5 kb) in our analysis. Fig 5 summarizes data on these questions for the complex centromeres, while data for the simple centromeres are shown in S14 Fig. Similar results were found when all elements were included (S14 Fig).

The results for the two dominant CRM families, CRM1 and CRM2, are similar (S14 Fig), so these families were grouped together in Fig 5C. When present in centromeres, all three major classes of elements–genes, CRM1/2, and CentC repeats–clearly replicate later during the endocycle than in the mitotic cycle (Fig 5). In contrast, genes and CRM elements in the pericentromere show little or no timing shifts. A full analysis of the replication times of CentC repeats in pericentromeres is hampered by the limited representation of this repeat class in the genome assembly (Figs 5D and S14E).

Chromatin features in centromeres

We also examined activating (H3K56ac and H3K4me3) and repressive (H3K27me3) histone H3 post-translational modifications to look for epigenetic changes in centromeres after endocycle replication. It was previously reported that some H3K4me3 and H3K27me3 peaks of enrichment occur in the centromere, mainly associated with genes [60]. We asked whether genes that have these modifications continue to have them after mitotic and endocycle replication, and found very few changes in the number of genes with these modifications at each ploidy level (S15 Fig). There was also very little change in the fold enrichment of these histone marks in centromere genes when comparing 2C, 4C and 8C nuclei.

We also investigated the levels of dimethylation of histone H3 lysine 9 (H3K9me2) enrichment in each centromere after mitotic and endocycle replication. Previous work indicated there is a depletion of H3K9me2 in centromeres relative to adjacent pericentromeres [61, 62], which we observed as well (S16 Fig). Traditional peak calling tools are not effective for H3K9me2 because of its even distribution across the maize genome. Instead, we estimated the fold enrichment relative to the corresponding DNA input control by calculating the percent of total H3K9me2 ChIP reads in a given centromere region (using coordinates from [46]) and dividing by the percent of total input reads corresponding to that centromere in three biological replicates. We found a similar H3K9me2 average depletion for all centromeres and for 2C, 4C and 8C nuclei relative to their corresponding input (averages for individual CENs 0.75–0.88, 0.83–0.91 and 0.85–0.93 for 2C, 4C and 8C nuclei, respectively), although values for 4C and 8C nuclei were consistently slightly higher than those for 2C nuclei (S16A Fig). The maize CENH3 N-terminal tail lacks some of the conserved motifs found in canonical histone H3 (see S3 Table in [63]), so H3K9me2 enrichment is likely to occur in the interspersed H3 nucleosomes. Altogether, we found very little change in centromeres between non S-phase 2C, 4C and 8C nuclei in the histone H3 post-translational modifications we assessed by ChIP-seq.

Centromeric histone H3 in mitotic and endocycling centromeres

Unlike the canonical histone H3, CENH3 is not replaced in a replication dependent manner in higher eukaryotes, resulting in a dilution of CENH3 relative to centromeric DNA during S phase [64, 65]. New CENH3 is incorporated into nucleosomes after the completion of S phase, but the timing of its integration into centromeric chromatin differs for plants, flies and humans (reviewed in [66]). In the plants tested thus far, deposition of CENH3 has been reported to occur between late G2 and metaphase [67–70].

Because mitosis does not occur in the endocycle and centromere function is presumably not required, we speculated that CENH3 might remain at low levels following DNA replication in endocycling cells. This hypothesis is supported by cytological studies of Arabidopsis endopolyploid nuclei showing the CENH3 signal does not increase in parallel with the total DNA content or the signal for 180-bp centromeric repeats [68, 69]. To test this hypothesis with maize centromeres, we used a maize anti-CENH3 antibody [58] for ChIP-seq analysis of CENH3 binding in sorted non S-phase 2C, 4C and 8C populations of nuclei. It is important to note that the 4C nuclei come from a mixture of cells, some of which will return to the mitotic cycle, while others will remain at 4C, and still others will continue on to the endocycle (at least 15% of nuclei in the combined 0–3 mm region). We asked whether the location or level of CENH3 enrichment changed after DNA replication in the mitotic cycle or the endocycle. For visualization of CENH3 localization, ChIP-seq uniquely mapping read counts from three biological replicates for each ploidy level were aggregated in 3-kb windows and normalized to the level of a uniform 1× genome coverage, so that corresponding windows in the different ploidy level profiles were comparable. The normalized read count in each 3-kb window was then divided by the normalized read count for input DNA of the corresponding ploidy to calculate a fold enrichment value for CENH3 binding sequences in that window. The spatial distribution of CENH3 enrichment across the centromeres remained the same in 2C, 4C, and 8C cells. This is illustrated for CEN 9 and CEN 10 in Fig 6A and 6B, and data for the rest of the centromeres are shown in S17 Fig. There are also a few small spikes of CENH3 enrichment outside the called centromere (e.g. seen in Fig 6 and S17 Fig, but also occasionally further out on the arms). These spikes also remain in the same location between 2C, 4C and 8C cells. Some of them could be related to misassembly of the reference genome. However, if real, these ectopic CENH3 peaks are less numerous and more persistent in G2 (4C) than those recently observed in HeLa cells [71].

To compare total CENH3 content of entire centromeres at different ploidy levels, we calculated the percent of total CENH3 reads found in a given centromere and made a ratio to the percent of total reads from the corresponding input DNA in that centromere separately for each biological replicate, as described above for H3K9me2. The CENH3 average fold enrichment relative to total DNA content is similar for 2C and 4C nuclei in each of the complex centromeres (Fig 6C), with an average 4C/2C enrichment ratio of 1.1 (S5 Table). However, CENH3 enrichment decreases with the increase in ploidy from 4C to 8C (Fig 6C). As noted above, the 4C nuclei come from a mixed population of cells, only a fraction (ca. 15% in the combined 0–3 mm root region) of which will enter the endocycle. Because of this ambiguity, we chose to focus on the CENH3 enrichment in 8C nuclei. In these nuclei, we found, an average 8C to 2C enrichment ratio of only 0.7 (S5 Table). CENH3 enrichment values for simple centromeres were lower and slightly more variable, likely because of assembly issues. In both cases, however, the ratio of CENH3 enrichment in 8C cells to that in 2C cells is higher than 0.5, the ratio that would be expected if there was no incorporation of new CENH3 after endocycle replication, but smaller than the 1.0 ratio expected if there was full replacement (S5 Table). It is worth noting that these data refer to post-replication 8C nuclei, which exited S phase prior to the time of analysis, and that post-replication 4C nuclei show no dilution of CENH3 relative to DNA content. Thus, our data are consistent with a hypothesis in which the average CENH3 to DNA ratio in the 8C population is only partially restored after completion of S phase.

Discussion

The maize root tip includes a naturally occurring developmental gradient, with cells in the meristem region (ca. 0–1 mm) primarily undergoing mitotic cell cycles, while a subpopulation of cells in the transition zone (ca. 1–3 mm) enters a developmentally programmed endocycle prior to further differentiation [8, 9]. Even though endocycling is very common in plants and plays essential roles in differentiation and the development of specialized tissues, cell size increases, and stress responses [2, 5, 72, 73], replication timing (RT) programs have not yet been characterized for alternative cell cycles, such as the endocycle.

We generated whole genome Repli-seq data for root cell nuclei undergoing DNA replication in either the mitotic cycle or the endocycle, making use of in vivo EdU labeling of intact root tips and two-color fluorescence activated nuclei sorting. By doing so, we avoided potential artefacts caused by cell synchronization [74] and chromosome aberrations often found in plant and animal cell cultures (e.g. [75–77]). We present replication timing activity profiles (RT profiles) for early, mid and late replication separately, instead of collapsing the data into an early:late ratio as many studies do. The rationale for this approach is that, for roughly one third of the maize genome, we previously found heterogeneity in mitotic RT–e.g. regions of the genome in which root tip cells exhibit significant replication activity in both early and mid S, or both mid and late S [19]. An additional advantage to presenting the RT profiles separately is the ability to assess whether there are concomitant or “compensated” changes in a region at multiple stages of S phase. This compensation criterion helped us separate RT shifts that could be subject to technical error, such as alterations in flow sorting gates, from shifts that are more likely to represent meaningful changes in the population preference to replicate a replicon or cluster of replicons at a particular time in S phase.

The current study sought to investigate whether the mitotic RT program is maintained as endocycling cells transition from 4C to 8C, despite the need to replicate twice as much DNA and the initiation of various root cell differentiation pathways. Extending our previous cytological observation that spatiotemporal patterns of replication are similar in mitotic and endocycling cells [22], we found that RT programs at the sequence level are also very similar. Pearson’s correlation coefficient values comparing data from the two types of cell cycles were nearly identical to those for biological replicates within each type. The high level of similarity is particularly noteworthy in the case of the early RT profiles, given that the flow sorting gate for early replicating nuclei in the endocycle had to be adjusted to minimize contamination from late replicating mitotic nuclei (Fig 1C). This overall conservation of the RT program in the endocycle is consistent with a recent study in Drosophila follicle cells which found no regions of differential RT between mitotic and endocycling cells [78]. Additionally, the global maintenance of RT in the two types of cell cycles in maize roots suggests that the process of re-establishing the RT program must be similar in both. In animal systems, re-establishment of the RT program has been shown to occur in G1 of each cell cycle at a “timing decision point” [79], however the details of this process have not been studied in plants.

Most plants fully replicate their genome during endocycles [80], although there are a few exceptions (e.g. various orchid species; [31, 32]). We found very little evidence for over- or under-replication occurring in endocycling maize root cells, unlike the distinctive over- and under-replication found in Drosophila endocycles (reviewed in [27] and references therein). Our result is consistent with earlier cytological reports that whole chromosomes, as well as repetitive knobs and centromeres, are completely replicated in the highly endopolyploid maize endosperm [33].

In contrast to the global maintenance of the RT program, we observed a small fraction of the maize genome that exhibits some difference in RT (DRT) between the two types of cell cycles. Approximately 11% of the genome showed compensated DRT at a stringency level of ≥ 10% difference in replication signal (see Methods). However, with the notable exception of centromeric regions, which are discussed in more detail below, we chose to characterize only the most robust Regions of Altered Timing (RATs), defined by the criteria of containing a core region with compensated DRT at a stringency level of ≥ 25% difference in replication signal. These robust non-centromeric RATs comprise only 1.6% of the genome (centromeric RATs comprise an additional 0.4% of the genome), and the size range of individual non-centromeric RATs (39–387 kb, median 138 kb) is consistent with our previous observation that regions of coordinate replication in maize are ~50–300 kb in size [19]. This may include from one to a few replicons, based on previous estimates of replicon size in monocot plants [49].

The first 1-mm of the maize root contains the meristem and precursors for at least ten different cell types. Only some of these cell types enter the endocycle prior to cell elongation [9]. If there are differences in the RT programs of different cell types, some or all of the non-centromeric RATs may be associated with shifts in the relative contribution of different cell types to the two samples of nuclei, rather than to endocycling per se. Research in metazoans has revealed ~8–20% of their genomes can shift RT between cell types [11–13, 81–83]. In mammals, these RT shifts generally involve large regions or “domains” in the megabase size range (reviewed in [26]). These RT domains are much larger than the non-centromeric RATs in maize, even though the maize genome is similar in size to the human and mouse genomes. However, in the much smaller Drosophila genome, regions that show RT shifts between cell types are more similar in size to the maize non-centromeric RATs [81, 83].

The vast majority of the non-centromeric RATs involved RT shifts from Earlier-to-Later, with a significant enrichment for genes expressed in the root tip, but not for non-expressed genes (FPKM < 1). This result suggests the possibility that RT shifts may be related to shifts in gene expression. Unfortunately, we have been unable to follow transcriptional changes in endocycling nuclei directly, as we have so far not been able to isolate RNA of sufficient quality to characterize transcripts from fixed, sorted nuclei. However, our analysis of activating and repressive histone modifications uncovered only minor changes in the enrichment and location of these marks within RAT genes after endocycle replication. The lack of notable changes in the proportion of RAT genes bearing H3K56ac and H3K4me3 modifications after the endocycle suggests that these histone marks are permissive to changes in RT. Nonetheless, the direction of the change in H3K4me3 enrichment on genes in Earlier-to-Later RATs after endocycle replication (S7B Fig) is consistent with the hypothesis that a shift to later RT may accompany a decrease in gene expression. Many studies have identified a correlation between RT and transcriptional activity (reviewed in [26]), but there are also multiple examples of these processes being uncoupled (e.g. [14, 84]).

In the case of centromeres, it is easy to imagine that the large shifts to later replication are related specifically to endocycling, because endocycling cells presumably no longer require functional centromeres. Though often broken by unmappable and multi-mapping (“blacklist”) regions in the genome assembly, centromeric RATs when combined across blacklist regions are much larger in size (0.9–2 Mb) than the non-centromeric RATs and cover the majority of each of the seven complex centromeres (Table 2). These seven centromeres, which are well assembled in the maize B73 RefGen_v4 genome, contain satellite repeats interspersed with retrotransposons [46, 56], enabling almost 50% of our sequencing reads that map to these centromeres to be uniquely positioned. In most species, in which centromeres contain large numbers of tandemly arrayed satellite repeats, it is difficult to map centromeric sequence reads to unique positions and, thus, to fully assess centromeric RT patterns [85]. Though yeast centromeres replicate in early S phase [86–89], most higher eukaryotes replicate centromeres asynchronously through mid to late S phase [64, 90–95]. Many of the reports in higher eukaryotes are based on cytological observations, membrane hybridization, or PCR data with limited resolution. Even a recent genomic analysis of centromeric RT in human cell lines was significantly limited by the quality of the human centromere assemblies, and could only uniquely map ~15% of centromeric reads [85]. Centromere replication in plant species, assessed mostly by cytological methods, has variously been reported to occur in early, mid or late S [96–99], though it is often unclear if the analysis was of sufficient resolution to distinguish the RT of centromeres from that of adjacent pericentromeres. In contrast, we have provided a high-resolution analysis of the distribution of replication times across maize centromeres, and compared RT of centromeres to adjacent pericentromeres.

These analyses revealed several features shared by the RT programs of the seven complex maize centromeres. For example, in mitotic cells there are a few distinct peaks of early replication (e.g. arrowheads in Fig 4 and S12 Fig), flanked by mainly mid replication activity that is, in turn, flanked by regions of late replication at the edges of the functional centromere. Except for these few early regions, which show an increase in mid replication activity in the endocycle, entire centromeres and the genes, retroelements and CentC repeats within them–undergo a shift to late replication in the endocycle. As a result, the RT of the complex centromeres in the endocycle becomes much more similar to that of the immediately adjacent pericentromeric regions, which replicate primarily in late S phase in both mitotic and endocycling cells. Late replication of pericentromeres is expected based on our previous cytological observations in mitotic and endocycling nuclei [22] and the typical replication time of highly compacted heterochromatin in many systems.

The presence of distinct peaks of early replication in or adjacent to functional centromeres (arrowheads in Fig 4 and S12 Fig) is noteworthy because they signify a population preference for replication initiation in early S phase at these loci. This observation is of particular interest because in the yeast Candida albicans, centromeres contain a replication origin that is the first to initiate on its respective chromosome and plays a role in centromere specification [89]. Work in other yeasts has shown evidence that the phenomenon of early firing origins in and near centromeres is present across a range of yeast phylogeny and that the presence and function of the centromere itself influences RT [100–102]. In maize, there is no evidence that these early regions in and adjacent to centromeres are the first to replicate on the entire chromosome, but they are earlier replicating than their surroundings. Origin mapping experiments (e.g. [103–105]) would be required to distinguish whether these early regions contain single or small clusters of origins, and the location of any other origins in centromeres that may initiate in mid or late S phase.

Unlike complex centromeres, the three simple centromeres of maize show less drastic RT changes that occur over smaller regions. These simple centromeres are not as well assembled as the complex centromeres [56, 57], and we cannot assess RT for the possibly large portions of these centromeres not present in the genome assembly. One potential interpretation of our results is that the simple centromeres have distinct RT programs that show less timing shift in the endocycle, possibly related to their different sequence composition. Alternatively, the missing portions of the simple centromere assemblies could be replicating more like the complex centromeres. Because simple centromeres are known to primarily contain large CentC arrays [56, 57], the second hypothesis is supported by our analysis of mapped CentC satellite repeats in all centromeres, which showed that, as a group, these repeats consistently shift RT from mid to late. Another piece of evidence comes from our analysis of complex centromeres, which showed that the magnitude of the RT change tapers off toward the outer edges of the functional centromere. One can speculate that the simple centromere assemblies are comprised mostly of the sequences at the edges of the actual centromere, which would still be anchored to nonrepetitive regions in the genome assembly. As in complex centromeres, these edge sequences might have a smaller RT shift than internal sequences. Future cytological experiments, using a combination of flow sorted, EdU-labeled nuclei and techniques for identifying maize chromosomes [106, 107] could help address questions related to the RT of simple centromeres.

The centromere-specific histone variant, CENH3 (also called CENP-A in animal systems) plays an important role in recruiting kinetochore proteins [53–55]. In metazoans, it has been shown that CENP-A is distributed among sister centromeres during replication, but the full complement of new molecules is not redeposited until later [65, 108]. However, there are differences in the timing of deposition of CENH3/CENP-A among eukaryotes. Deposition occurs from S phase to G2 in yeasts, while in plants and protozoans it occurs from late G2 to metaphase, and in metazoans it occurs mostly during G1 (with the exception of some Drosophila cell types in metaphase to G1; reviewed in [55, 66, 70]). These interesting differences between phylogenetic groups in the timing of CENH3/CENP-A deposition suggest there may also be differences in the mechanisms and regulation of deposition that need to be explored further [69]. In our analysis of CENH3 enrichment relative to DNA content in maize root cells, the population of 4C nuclei appear to have a full complement of CENH3, which would be consistent with the previous results for plant species. This result supports a hypothesis that the sub-population of 4C cells entering the endocycle also carry a full complement of CENH3. If this hypothesis is correct, our data for 8C nuclei imply that CENH3 is only partially replaced after DNA replication in the endocycle. Because the population of 8C nuclei we analyzed likely represents a mixture of cells that recently exited the endocycle S phase and others that exited some time ago we cannot determine whether CENH3 is fully restored in some cells and not others, or if it might be fully restored at a later time. However, our data suggest that the ratio of CENH3 to DNA is not immediately restored in the 8C population, and that the lower ratio is widely distributed across all ten centromeres.

It is unlikely that endocycling cells will ever re-enter the mitotic cycle [1, 109, 110], and it is not clear why endocycling cells would maintain or redeposit CENH3 nucleosomes at all unless CENH3 has roles outside of mitotic cell division. A recent study in Drosophila midgut cells found that CENP-A is required in post-mitotic and differentiated cells, and proposed that the loading of CENP-A in endocycling cells is essential for maintaining chromosome cohesion [111]. This possibility has not yet been tested.

Centromeres are considered to be epigenetically specified, as there are no unique sequences in the functional centromere that are not also found in the adjacent pericentromere (e.g. reviewed in [53, 112]). With this in mind, we tested whether changes in enrichment levels of CENH3 nucleosomes, or several modifications to canonical H3 nucleosomes, could explain the large shift to later replication of centromeres in endocycling cells. These studies only uncovered very small changes in activating and repressive histone H3 modifications in centromeres after endocycle replication. The magnitude of the change in CENH3, while somewhat larger, was not on the scale of the change in RT. It is possible that more significant changes might be found in epigenetic factors that we did not investigate, for example changes in DNA methylation patterns or other histone post-translational modifications. A variety of modifications to CENP-A nucleosomes have been identified, (reviewed in [113]), but very little is known about CENH3 modifications in plants [114, 115], highlighting an area for future research. Experiments in human cells identified cell cycle related interchanges of acetylation, monomethylation and ubiquitination at the lysine 124 residue of CENP-A [116, 117]. Mutations of this residue led to replication defects and alterations to centromeric RT [117].

An interesting question for future investigation is whether changes in chromatin conformation or 3D positioning in the nucleus are associated with the large shift in centromeric RT. In mammals, RT is considered a functional readout of large-scale chromatin structure [14, 26, 82], and regions that shift RT have been shown to also change 3D localization [118]. Likewise, a study in mouse showed that when late replicating pericentric heterochromatin was experimentally repositioned to the nuclear periphery, a location where mid replicating chromatin is usually found in that system, the RT of those regions was advanced [119].

Additionally, we speculate that centromere transcription could play a role in the shift to later replication of centromeres in the endocycle. Non-coding transcription has been reported to occur in centromeres across a wide range of eukaryotes, and has also been found to be essential for centromere function ([112] and references therein). Furthermore, in human cells, when centromere RNAs transcribed from alpha satellite repeats are specifically degraded, cells arrest before mitosis and have reduced CENP-A levels [120]. In maize, low level transcription has been detected from CentC satellite repeats as well as CRM retroelements, and there is evidence that centromere transcripts are involved in the function of the kinetochore complex and centromere chromatin organization [121–123]. As noted above, technical limitations prevented us from assessing transcripts from endocycling cells. However, it seems possible that a reduction in centromere transcription during the endocycle could contribute to the observed reduction in CENH3 and potential dismantling of the kinetochore function. Outside the centromere region, reduced transcription is often associated with late replication [19, 26]. Future work addressing the technical challenge of sorting maize nuclei for transcript profiling will be needed to test this hypothesis.

One further interesting possibility arises from the observation that small regions in or adjacent to the 7 complex centromeres exhibit early replication during the mitotic cycle which is reduced in the endocycle. In mitotic cells, regions flanking these early peaks often show strong mid-S replication with relatively little replication in late S. In endocycling cells there is still very little late replication at the actual locus of the early peak, but late replication is enhanced in the flanking centromere regions. These small, early-replicating regions are reminiscent of the early replicating origins described in and near centromeres in several types of yeast [89, 101] as well as the “initiation regions” identified in very early S phase recently described in Arabidopsis by Wheeler et al. [105]. Strikingly, some of the Arabidopsis very early initiation regions occur in centromeric/pericentromeric heterochromatin of Arabidopsis chromosomes. It is possible to imagine that these early replicating regions at maize centromeres contain early-firing origins that drive replication of bulk centromeric chromatin during mid-S in mitotic cells, but that in endocycling cells, initiation activity in these regions is reduced. Alterations in chromatin conformation, subnuclear positioning, or other factors could mediate such a reduction in early initiation activity, which in turn could cause, or contribute to, the reduction in overall centromere transcription hypothesized in the previous paragraph.

The three speculative hypotheses presented above are not mutually exclusive, and in fact can easily be imagined to work synergistically together. We believe all three deserve further investigation as appropriate techniques can be developed. Investigating the interplay of chromatin environment, subnuclear organization, transcription and DNA replication in plant systems has proven difficult in the past. Numerous reasons for these difficulties exist, for example, plants have cell walls and are rich in nucleases, actively dividing cells are sequestered in tiny meristematic regions, and many genomes have a high content of retrotransposons and other repeats. As a result, understanding of such critical areas has lagged behind that in yeast and animal systems. However, with recent progress in assembling genomic resources and anticipated advances in the ability to isolate individual cell types [124], perform sophisticated analyses of genome conformation [125, 126] and follow individual chromosome regions using elegant cytological paints [107], the maize root tip system is poised to contribute to rapid progress in these and many other important areas of plant genome biology.

Methods

Plant material

Seeds of Zea mays inbred line B73 (GRIN NPGS PI 550473) were germinated on damp paper towels and grown for three days. Seedling roots were labeled by immersion in sterile water containing 25 μM EdU (Life Technologies) for 20 min, using growth and experimental conditions described previously [8, 19, 20]. Biological replicate material was grown independently and harvested on different days. For the endocycle Repli-seq experiment, after rinsing roots well with sterile water, the 1–3 mm segments (Fig 1A) were excised from primary and seminal roots. The root segments were fixed, washed and snap-frozen as described previously [20].

Flow cytometry and sorting of root nuclei

Details of the flow sorting for Repli-seq analysis were described previously [19, 20]. Briefly, nuclei were isolated from the fixed root segments, and the incorporated EdU was conjugated to AF-488 using a Click-iT EdU Alexa Fluor 488 Imaging Kit (Life Technologies). The nuclei were then resuspended in cell lysis buffer (CLB) [20] containing 2 μg/mL DAPI and 40 μg/mL Ribonuclease A and filtered through a CellTrics 20-μm nylon mesh filter (Partec) just before flow sorting on an InFlux flow cytometer (BD Biosciences) equipped with UV (355 nm) and blue (488 nm) lasers. Nuclei prepared from the 1–3 mm root segments were sorted to collect populations of EdU/AF-488-labeled nuclei with DNA contents in three defined sub-stage gates between 4C and 8C, corresponding to early, mid and late S phase of the endocycle. The early endocycle gate was shifted slightly to the right to exclude mitotic nuclei in late S phase (Fig 1C). For each biological replicate, between 50,000 and 200,000 nuclei were sorted from each fraction of the endocycle S phase. A small sample of nuclei from each gate was sorted into CLB buffer containing DAPI and reanalyzed to determine the sort purity (S1 Fig). Sorting and reanalysis details for the mitotic nuclei are described in [19].

For ChIP-seq experiments, roots were labeled with EdU, and nuclei were isolated from 0–3 mm (H3K27me3 and H3K4me3) or 0–5 mm (H3K56ac) root segments and conjugated to AF-488 as described above. The 2C, 4C and 8C unlabeled, non S-phase populations of nuclei were sorted into 2× extraction buffer 2 (EB2) [127] using the same sorting conditions as in Wear et al. [19]. After sorting, the 2× EB2 was diluted to 1× with 1× STE. All flow cytometry data were analyzed using FlowJo v10.0.6 (TreeStar, Inc.) as described in Wear et al. [19].

DNA and chromatin immunoprecipitations

For endocycle Repli-seq samples, reversal of formaldehyde cross links, nuclear DNA purification and isolation, DNA shearing, EdU/AF-488 DNA immunoprecipitation with an anti-Alexa Fluor 488 antibody (Molecular Probes, #A-11094, lot 895897), and DNA fragment purification were performed as described in Wear et al. [19].

ChIP procedures were performed as in Wear et al. [19] except the chromatin was sheared using a Covaris S220 ultrasonicator to an average fragment size of 200 bp using a peak incident power of 140 W, 10% duty cycle, and 200 cycles per burst for 6 min. Three percent of the chromatin volume was set aside to use as the input control for each of the 2C, 4C and 8C samples and frozen at -70°C until the formaldehyde cross link reversal step. The antibodies used for ChIP were as follows: Zea mays anti-CENH3 antibody at a 1:250 dilution (gift from R.K. Dawe) [58], anti-H3K9me2 antibody at a 1:25 dilution (Cell Signaling Technologies; 9753, lot 4), anti-H3K56ac antibody at a 1:200 dilution (Millipore; 07–677, lot DAM1462569), anti-H3K4me3 antibody at a 1:300 dilution (Millipore; 07–473, lot DAM1779237) and anti-H3K27me3 antibody at a 1:300 dilution (Millipore; 07–449, lot 2,275,589). See S18 Fig for antibody validation experiments for anti-H3K9me2 and anti-CENH3.

Library construction and sequencing

For Repli-seq and ChIP-seq samples, the final purified DNA was used to construct paired-end libraries as described [19]. After adapter ligation, all samples underwent 17 cycles of PCR. For each Repli-seq or ChIP-seq experiment, individual samples from three biological replicates collected on different days were barcoded, pooled and sequenced on either the Illumina HiSeq 2000 or NextSeq platforms. However, in the case of the Repli-seq mitotic late-S samples and CENH3 ChIP 4C samples, one biological replicate failed during library generation or sequencing, resulting in data from only two biological replicates. Repli-seq and ChIP-seq read mapping statistics are shown in S1 Spreadsheet.

Replication timing data analysis

Trimming and quality control of 100-bp paired-end Repli-seq reads were carried out as described previously [19], and reads were aligned to the maize B73 RefGen_v4 reference genome [46] (Ensembl Plants release 33; ftp://ftp.ensemblgenomes.org/pub/plants/release-33/gff3/zea_mays/) using BWA-MEM v0.7.12 with default parameters [128]. Redundant reads resulting from PCR amplification were removed from each of the alignment files using Picard (http://broadinstitute.github.io/picard/) and SAMtools [129]. Properly paired, uniquely mapping reads (MAPQ score > 10) were retained with SAMtools [129] for downstream analysis. The resulting mitotic Repli-seq data were more than 3× the sequencing coverage of the endocycle Repli-seq data (S1 Spreadsheet). Repli-seq results are robust at various sequencing depths [24], but to ensure that the mitotic and endocycle data were comparable, the reads were downsampled by a uniform random process using a custom python script incorporating the BEDTools suite [130] to a total of 65.7 million reads per S-phase fraction (S1 Spreadsheet). We preferred this to normalization so that any possible sampling bias due to sequencing depth would be similar in all samples.

Repli-seq data were analyzed using Repliscan [24]. Individual biological replicates of Repli-seq data were independently analyzed, and after finding good correlation between replicates (Pearson correlation coefficients from 0.80–0.99; S4 Fig) the replicates were aggregated by sum and normalized to 1× genome coverage using the reads per genomic content (RPGC) method. The following changes from the Repliscan default parameters described in [19] were used. Read densities were aggregated in 3-kb windows across the genome (parameter -w 3000). Additionally, we customized the cutoff for reducing type one errors which excluded genomic windows with extremely low coverage in the 2C reference sample. To identify these low read mapping windows, which we labeled “blacklist”, Repliscan natural log-transformed the read counts from the pre-replicative 2C reference sample and windows with read counts in the lower 2.5% tail of a fitted normal distribution were excluded from all samples (parameter—pcut 2.5–100). The upper 2.5% tail containing extremely high coverage windows or “spikes” was not removed at this step, because we found that these data spikes were adequately normalized in the subsequent step of dividing each 3-kb window in the S-phase samples by the 2C reference data–which also normalized for sequencing biases and collapsed repeats (S3 Fig). The data were then Haar wavelet smoothed [24] to produce the final profiles for early, mid and late S-phase replication signals in the mitotic cycle and endocycle. Processed data files, formatted for the Integrative Genomics Viewer (IGV) [131], are available for download from CyVerse (formerly the iPlant Collaborative; [132]) via the information in S1 Spreadsheet.

Identifying regions of altered replication timing

The difference between normalized replication signal profiles of mitotic and endocycle Repli-seq data for early, mid, and late S was calculated in 3-kb windows, and the maximum negative and positive differences were then calculated for each chromosome and averaged. Regions showing a DRT of ≥ 25% (absolute difference in replication signal ≥ 1.0) or ≥ 10% (absolute difference in replication signal ≥ 0.4) of the total range of differences in each profile were identified (S1 Table; S5 Fig) using the data filter tool in SAS JMP Pro v14 (SAS Institute Inc.). Windows were kept in the analysis only if their DRT were “compensated” by opposite DRT of ≥ 25% or ≥ 10%, respectively, in one or both of the other two S-phase fractions. For example, a decrease in early replication signal in endocycling cells must be compensated by an increase in mid and/or late S-phase signal in the same cell population. Adjacent 3-kb windows with DRT that met either the ≥ 10% or ≥ 25% threshold were merged, keeping the two files separate, using mergeBED in the BEDTools suite, and allowing a 6 kb gap distance (parameter -d 6000) [130]. This initial step resulted in many very small regions being identified (S2 Table). As a second step, if ≥ 10% regions were immediately adjacent to ≥ 25% regions, they were merged together using mergeBED to highlight larger regions of contiguous change (S2 Table). Only regions that contained at least one ≥ 25% region were kept for further analysis, and termed regions of alternate timing (RATs). By requiring a ≥ 25% DRT core region to be included, all of the stand-alone, extremely small regions (< 24 kb) were effectively filtered out, without the requirement of an arbitrary size filter. RATs were categorized into three groups: 1) later in mitotic to earlier in endocycle (Later-to-Earlier), 2) earlier in mitotic to later in endocycle (Earlier-to-Later) and 3) a subset of the Earlier-to-Later RATS that were located in the previously identified functional centromeres (Earlier-to-Later-CEN) (coordinates from [46]). There were no Later-to-Earlier-CEN RATs. For a list of RAT regions, including genomic coordinates and genes within them, see S2 and S3 Spreadsheets.

ChIP-seq data analysis

ChIP-seq reads for H3K27me3, H3K4me3, H3K56ac (100-bp paired-end reads), H3K9me2 and CENH3 (150-bp paired-end reads) were trimmed, mapped to maize B73 RefGen_v4.33, and filtered to retain only properly-paired, uniquely-mapped reads (MAPQ score > 10) as described above for Repli-seq reads. The 2C ChIP and input data for H3K27me3, H3K4me3, H3K56ac is from [19], while the 4C and 8C ChIP data was generated for this study, see S1 Spreadsheet. For details on peak calling and analysis for H3K27me3, H3K4me3, H3K56ac, see S1 Text.

For visualization of CENH3 localization in 2C, 4C and 8C nuclei, read counts for individual biological replicates of CENH3 or input samples were scaled to 1× genome coverage using the reads per genomic content (RPGC) method. Biological replicate data had good agreement (Pearson’s correlation coefficient values between biological replicates of 0.97–0.99; S1 Spreadsheet), and were merged and scaled again to 1× coverage so the samples would be comparable. CENH3 scaled read counts in each 3-kb window were divided by the scaled read counts from the input sample for the corresponding ploidy level, resulting in CENH3 fold enrichment values relative to input.

To compare CENH3 enrichment relative to DNA content in 2C, 4C and 8C cells over entire centromeres, we calculated the percent of total CENH3 reads found in a given centromere (using coordinates from [46]), and divided by the percent of total input reads corresponding to that centromere. This was done separately for individual biological replicates; we then calculated the mean fold enrichment estimates. H3K9me2 fold enrichment over entire centromeres and pericentromeres was calculated in the same way.

Genomic features

The maize filtered gene set Zm00001d.2 annotation from B73 RefGen_v4 [46] was downloaded from Ensembl Plants (ftp://ftp.ensemblgenomes.org/pub/plants/release-33/gff3/zea_mays/). The updated B73 Refgen_v4 TEv2 disjoined annotation [47] was downloaded from http://mcstitzer.github.io/maize_TEs. Coordinates for mapped CentC satellite repeat regions are described in Gent et al. [57]. The percent AT content was calculated in 3-kb static windows across the genome.

Analysis of features in RATs and random permutation analysis

We tested the association of various genomic features with the non-CEN RAT categories by determining the overlap of a particular feature with each RAT type. The coordinates for genomic features (genes, expressed genes, TE superfamilies) were intersected with RAT coordinate intervals using intersectBED (parameters -wa -wb) in the BEDtools suite [130]. The percent of RATs containing a feature or the percent coverage of TE superfamilies were then computed. The number of genes per RAT was also determined using intersectBED (parameter -u).

For comparison, the coordinates for the non-CEN Earlier-to-Later and Later-to-Earlier RAT sets were randomly shuffled around the genome, excluding functional centromeres, using BEDTools shuffle [130]. These random sets preserved the number of regions and region size of the original RAT sets, and are labeled “EtoL shuffle1” and “LtoE shuffle1” for the Earlier-to-Later and Later-to-Earlier RATs, respectively (e.g. S7, S8 and S10 Figs). When there appeared to be differences in the observed overlap values with genomic features between non-CEN RATs and their corresponding random shuffle sets, a permutation or feature randomization test, as described in [19] was used to assess the statistical significance of the observed value. To do so, the coordinates for the non-CEN RAT sets were randomly shuffled around the genome 1000 times, as described above, to generate an “expected” distribution. Permutation P values, calculated as described in [19], that were ≤ 0.001 were considered evidence that the observed percent coverage value is significantly different than random expectation.

Analysis of features in centromeres and pericentromeres

For comparison to CEN regions (coordinates from [46]), pericentromeres were arbitrarily defined as the ± 1 Mb flanking each CEN. In the case of chromosome 9, the pericentromere included the ± 1 Mb flanking both CEN 9a and 9b. Replication signal values in CENs and pericentromeres were intersected with genes, CRM1 and CRM2 families and mapped CentC regions using intersectBED (parameters -wa -wb) in the BEDtools suite [130]. Only elements that covered at least half of a 3-kb window of Repli-seq data were included in Fig 5, while elements with any amount of overlap were included in S14 Fig. Additionally, if a single gene or CRM element spanned more than one of the 3-kb windows, the replication signals were averaged using mergeBED (parameter -o mean) to compute a single value for the entire gene or element.

Supporting information

S1 Text [docx]
Supplemental Methods.

S1 Fig [e]
(related to ) Assessment of purity of flow sorted endocycling nuclei.

S2 Fig [a]
(related to ) Genomic copy number analysis.

S3 Fig [a]
(related to Figs and ) Example of Repli-seq data processing with .

S4 Fig [br]
(related to ) Pearson’s correlation coefficient values between individual biological replicates of mitotic and endocycle Repli-seq data.

S5 Fig [drt]
(related to ) Boxplots of differences in early, mid and late replication signal profiles for each chromosome.

S6 Fig [a]
(related to ) Additional examples of non-CEN RATs.

S7 Fig [a]
(related to ) Activating and repressive histone marks in non-CEN RATs.

S8 Fig [p]
(related to ) Gene ontology analysis of genes in non-CEN RATs.

S9 Fig [pdf]
(related to ) Permutation analysis of the percent coverage of TE superfamilies in non-CEN RATs.

S10 Fig [a]
(related to ) AT content composition in non-CEN RATs.

S11 Fig [pdf]
(related to ) Uniquely mapping Repli-seq reads in centromeres.

S12 Fig [pdf]
(related to ) Replication signal profiles and RATs in complex and simple centromeres.

S13 Fig [pdf]
(related to ) RT differences in centromeres and pericentromeres.

S14 Fig [e]
(related to ) Replication times for all genomic features in complex and simple centromeres and corresponding pericentromeres.

S15 Fig [a]
(related to Figs and ) Activating and repressive histone mark peaks of enrichment in centromeres.

S16 Fig [pdf]
(related to Figs –) H3K9me2 fold enrichment relative to DNA content in complex and simple centromeres.

S17 Fig [a]
(related to ) CENH3 localization and enrichment in mitotic and endocycling centromeres.

S18 Fig [a]
(related to ) ChIP-qPCR antibody validations for anti-CENH3 and anti-H3K9me2 antibodies.

S1 Table [docx]
(related to ) Replication timing signal differences and thresholds.

S2 Table [s]
(related to ) Summary statistics of preliminary RAT calling steps.

S3 Table [docx]
(related to Figs and ) Gene summary in non-CEN RATs.

S4 Table [docx]
(related to Figs and ) Compensated differences in RT in complex centromeres and corresponding pericentromeres.

S5 Table [docx]
(related to ) CENH3 average fold enrichment relative to DNA content in centromeres.

S1 Spreadsheet [xlsx]
(related to Figs –) Mapping statistics and data availability for all included datasets.

S2 Spreadsheet [xlsx]
(related to Figs and ) RAT regions list.

S3 Spreadsheet [xlsx]
(related to Figs –) Genes found in RATs.