Heterochromatin

For the estimation of the level of LD, a total of 31,954 loci with minor allele frequency >0.10 and the number of missing data points less than 25% was used. Heterozygous alleles were treated as missing data. Only physically linked markers located within 20 Mbp distances were used for LD estimation. Haploview 4.2 [52 (link)] was used to make all pair-wise comparisons of the alleles to calculate r² (the squared allele frequency correlation between two loci), and to compute D’ (standardized disequilibrium coefficient) [53 (link)]. In addition, LD (D’ and r²) was estimated separately for euchromatic and heterochromatic regions. The euchromatin and the heterochromatin regions of each of the 20 chromosomes were determined as follows: The physical positions of 3,321 SNP and 862 simple sequence repeat (SSR) markers mapped in the soybean genome [21 (link),51 (link),54 (link)-56 (link)] were determined by BLAST analysis of the SNP and SSR-containing source sequences to the soybean whole genome sequence using the standalone Megablast software as previously described [57 (link)]. The cumulative genetic distances (cM) [21 (link)] were plotted against their cumulative physical distance (Mbp) to determine the base pair/centiMorgan relationship via the common SSR and SNP loci positions on the genetic linkage map and their physical position in the genome sequence along each chromosome. The region between the two inflection points of the cumulative genetic distance against cumulative physical distance on the plot was defined as the heterochromatin and the regions on each chromosome flanking the inflection points were defined as the euchromatin [45 (link)]. In order to provide an assessment of the difference in the extent of LD between euchromatic and heterochromatic regions, LD was calculated by a pairwise comparison of physically linked heterochromatic SNPs, and then two separate LD calculations were made for all pairs of markers from the two flanking euchromatic regions on each chromosome. The mean value of LD was estimated by calculating the mean LD of SNP pairs at distances of 0-200 Kbp, 200-400 Kbp, etc. to 19,800-20,000 Kbp in euchromatic and heterochromatic regions.

Sequencing was performed using standard protocols for the Illumina Genome Analyzer IIx. Initial data processing and quality analysis was performed using the standard Illumina pipeline. Sequence reads were deposited in the NIH Short Read Archive as project SRP005599. Alignments to the D. melanogaster reference genome (BDGP release 5) using BWA version 0.59 [31] (link) with default settings and the “-I” flag. Program defaults included a 32 bp seed length; reads could therefore map to the reference only if two or fewer reference differences were present within a seed. Although read lengths varied from 76 bp to 146 bp within this data set, only the first 76 bp of longer reads was used for the assemblies reported here. In order to exclude ambiguously mapping reads, those with a BWA mapping quality score less than 20 were eliminated from the assemblies.
Consensus sequences for each assembly were obtained using the SAMtools (version 0.1.16) pileup module [32] (link). These diploid consensus sequences generally included a few thousand heterozygous calls, scattered across the genome. Such sites are not expected to represent genuine heterozyosity in these haploid/homozygous samples (with the exception of ZK, in which large-scale heterozygosity was observed, presumbaly due to incomplete inbreeding). All putatively heterozygous sites were masked to ‘N’. Sites within 5 bp of a consensus indel were also masked to ‘N’ – this criterion was found to reduce errors associated with indel alignment; no appreciable benefit was observed if 10 bp was masked instead (data not shown).
Data were only considered for “target” chromosome arms, as defined in Table S1. These are chromosome arms expected to derive from the population sample of interest (as opposed to originating from laboratory balancer stocks), and observed to be free of heterozygous intervals. Chromosome arms were further defined as “focal” (the genomic regions analyzed here, namely the euchromatic portions of X, 2L, 2R, 3L, and 3R) or “non-focal” (the mitochondria and heterochromatin, including chromosomes 4 and Y). The assemblies analyzed here were defined as “release 2” data and are available for download at http://www.dpgp.org/dpgp2/DPGP2.html. Assemblies of mitochondrial and bacterial symbiont genomes are reported and analyzed separately [47] .

Stack images of the roots (1 µm section) were acquired using a ×40 oil objective with a NIKON A1R + confocal microscope. Tile-scanning (4 × 1 tiles) was used to ensure the imaging of the whole meristem and the elongation zone of the roots. The Z Project tool (Projection type: Sum slices) was used to sum the fluorescence intensity of the pixels corresponding to each nucleus present in the epidermal layer. The background fluorescence was subtracted for each color channel and the fluorescent intensity of chromocenters and nuclei was measured as the integrity density of a determined ROI. In all cases, independent measurements were taken for each color channel. Data acquisition and processing were done using FIJI. The Relative Heterochromatin Fraction (RHF) was calculated using the formula RHF = [Ac × (Ʃ Ic − Ib)] / [An × (In − Ib)]. Ac: the total chromocenters area; Ic: fluorescence intensity of chromocenter; Ib: fluorescence intensity of the background; An: area of the nucleus; In: fluorescence intensity of the nucleus.

The genomic regions were defined as chromosome arms and centers based on previous works^{12 (link),83 (link),84 (link)}. Briefly, recombination rates and H3K9 methylation enrichment were used to estimate the physical boundaries between the arms and central regions of each chromosome. Chromosome arms were characterized by high levels of recombination rate and H3K9 methylation. Coordinates of arms and center regions of each chromosome were shown in Supplementary Table 4.
Heterochromatin was defined by occupancy of H3K9me2, H3K9me3, and H3K27me3, and euchromatin was determined by covering of H3K4me3, H3K36me3, and H3K79me3^{12 (link),85 (link)–88 (link)}. Association between chromodomain proteins and histone modifications was detected by correlation, heatmap, and peak overlapping analysis.
For correlation analysis, combined peaks were determined by peaks called from any of the chromodomain protein or histone modification ChIP-seq datasets. The overlapping peaks were merged (n = 42,496). Then the merged peaks were divided into chromosome arm (n = 24,284) and center peaks (n = 18,212). The 95th percentile values were extracted from bigWig files using the Python package pyBigWig over each chromosome arm and center peak. The values were normalized to input signals, logarithm-transformed, and then standardized to Z scores. Correlation coefficients were calculated by the cor [Pearson] function in R using the treated values. A two-sided t test was performed.
For heatmap analysis in Fig. 2b–e, Supplementary Figs. 8d–g, 9a, the called peaks of each chromodomain proteins were used. Deeptools subcommand computeMatrix (version 3.4.3) was used to calculate score matrix for heatmap with defined parameters (reference-point --referencePoint center -b 3000 -a 3000 --skipZeros). The heatmap was plotted with Deeptools subcommand plotHeatmap (version 3.5.0). To identify heterogeneity of chromatin states of chromodomain protein targets, heatmap of each chromodomain protein was clustered. Parameter --kmeans was used to define the number of clusters. For each chromodomain protein, --kmeans 1, 2, and 3 were all applied. Plots displaying no heterogeneity of chromatin states within each cluster were selected. A plot with minimum kmeans number for each chromodomain protein from the selected plots was shown.
Significant overlapped peaks of chromodomain proteins and histone modifications were defined by using IntervalStats software package^{89 (link)}. The method compared each single peak region from a ‘query’ experiment to the set of peak regions in a ‘reference’ experiment. P-values represents the significance of query peak proximity to the reference peak. Pairs of peaks were considered significantly overlapped with P < 0.05. Both chromodomain proteins and histone modifications were used as queries and references (Fig. 2a, Supplementary Fig. 8a–c). The percentage of overlap was listed in Supplementary Table 5. We used the percentage of PTM-covered peaks of each chromodomain protein for “graded histone occupied” evaluation. For many chromodomain proteins, PTM covered more than 30% of their peaks. However, for PTMs, chromodomain proteins rarely covered more than 30% (Supplementary Fig. 8b, c, Supplementary Table 5). We speculated that for a given PTM, there are a number of additional readers besides chromodomain proteins, such as Tudor and MBT, etc., may function redundantly to recognize the PTM.
“Graded histone modification occupied” for a pair of chromodomain protein and histone modification fell into 3 categories: grade 3 = “prominent”, grade 2 = “detectable”, and grade 1 = “weak”. Grade 3 was defined as pairs meeting all of the 3 criteria: 1. “Pearson correlation coefficient (r) ≥ 0.300”. 2. “ ≥30% of significant peaks occupied by a certain histone modification (using IntervalStats software package with P < 0.05)”. 3. “Overlap detected by heatmap analysis”. Grade 2 meets 2 out of the 3 criteria and grade 1 meets any of the 3. For example, on chromosome arms, Pearson correlation coefficient of CEC-5 GFP(A) and H3K9me2, r = 0.333 (Supplementary Fig. 10); 44.3% peaks of CEC-5 were significantly occupied by H3K9me2 (Fig. 2a and Supplementary Table 5); Heatmap analysis found 1585 out of 1854 targets of CEC-5 were covered by H3K9me2 (Fig. 2b). Therefore, the association of CEC-5 and H3K9me2 on chromosome arms was classified as “grade 3 = prominent” (Fig. 2f and Supplementary Table 6).

To identify actively transcribed genes and silent genes, we analyzed 2 repeats of mRNA-seq datasets of adult N2 worms. The average of Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) of each gene was calculated. Actively transcribed genes were defined by “FPKM $\ge$ 1” and being occupied by any of the euchromatin marks (H3K4me3, H3K36me3, or H3K79me3). Silent genes were defined by “FPKM < 1” and being occupied by any of the heterochromatin marks (H3K9me2, H3K9me3, or H3K27me3).