The Dynamic Regulatory Genome of Capsaspora and the Origin of Animal Multicellularity

Abstract

The unicellular ancestor of animals had a complex repertoire of genes linked to multicellular processes. This suggests that changes in the regulatory genome, rather than in gene innovation, were key to the origin of animals. Here, we carry out multiple functional genomic assays in Capsaspora owczarzaki, the unicellular relative of animals with the largest known gene repertoire for transcriptional regulation. We show that changing chromatin states, differential lincRNA expression, and dynamic cis-regulatory sites are associated with life cycle transitions in Capsaspora. Moreover, we demonstrate conservation of animal developmental transcription-factor networks and extensive network interconnection in this premetazoan organism. In contrast, however, Capsaspora lacks animal promoter types, and its regulatory sites are small, proximal, and lack signatures of animal enhancers. Overall, our results indicate that the emergence of animal multicellularity was linked to a major shift in genome cis-regulatory complexity, most notably the appearance of distal enhancer regulation.

Graphical abstract

Figure 1

Histone Modifications in Capsaspora Histone N-terminal tail sequences of Capsaspora with the identified posttranslational modifications are shown. Below: filled or empty circles indicate whether the particular histone mark is present or absent, respectively, in the different eukaryotic species represented in the phylogenetic tree (left). Above: the presence (green) or absence (red) of specific histone modifiers in the Capsaspora genome is shown; both enzymes that add the mark (writers) and enzymes that remove it (erasers) are indicated. Capsaspora-specific marks are highlighted in yellow. The repressive marks H3K9me and H3K27me3 are absent in Capsaspora and indicated separately in a box below the corresponding position. See also Figure S1.

Figure 2

Genome-Wide Chromatin Annotation in Capsaspora (A) Top: TSS-centered average normalized read coverage plots of hPTMs in the filopodial stage for genes with high (green), intermediate (yellow), and low (red) expression levels. The x axis spans −5 to +5 kb around the TSS. The shaded gray area represents the average size of Capsaspora genes. Bottom: scatterplots of hPTMs coverage (log2 normalized reads) compared to mRNA expression levels (log2 fragments per kilobase of transcript per million mapped reads [FPKMs]). (B) Heatmaps of ATAC mononucleosome-associated (left) and nuclesosome-free (right) reads centered around the TSS of genes sorted by level of expression in the filopodial stage. Right: histogram showing an example of the distribution of ATAC-seq fragment sizes obtained. (C) Boxplot representing the mean fuzziness score of the first four post-TSS nucleosomes of genes grouped by the level of expression in the filopodial stage. The p value is indicated for the Wilcoxon rank-sum test. (D) Heatmaps representing the emission (left) and transition (right) parameters of a seven-state hidden Markov model. In the left heatmap, the white-blue (0–1) scale represents the frequency with which a given mark is found at genomic positions corresponding to the chromatin state. In the right heatmap, the white-blue (0–1) scale represents the frequency with which a given state changes into another state at the neighboring location. (E) Chromatin signatures in active genes (>2 FPKMs) in the filopodial stage. The plot (left) represents the average normalized read coverage of histone modifications around the TSS of these active genes, and the heatmap (right) indicates the relative percentage of the genome represented by each chromatin state (first column) and relative fold enrichment for different genome features (other columns). (F) Chromatin signatures in silent genes in the filopodial stage (heatmap and plot as in C). (G) Boxplot representing the expression levels in the filopodial stage of genes (left) selected for having a significant peak of H3K27ac in the gene body (more than 800 bp from the TSS) and a significant peak of H3K4me1 after the TSS (within 800 bp), and vice versa (right). The p value is indicated for Wilcoxon the rank-sum test. See also Figures S2 and S4 and Data S1.

Figure 3

Dynamic Chromatin Modifications (A) Boxplots showing hPTMs coverage levels in differentially expressed genes between stages, as indicated above each boxplot. The p value is indicated for the Wilcoxon signed-rank test. (B) Illustrative examples of dynamic chromatin modifications in Capsaspora. Different genomic windows show normalized coverage for different chromatin features and their dynamic association with gene expression. For each feature, the top track corresponds to the filopodial stage, the middle track to the aggregative stage, and the bottom track to the cystic stage. (C) Histone deacetylase inhibition experiments. Pictures of Capsaspora cells at different time points of incubation with DMSO (negative control) and TSA 3 μM. Transition from cystic to filopodial stage is blocked in the TSA-treated cells. Scale bar, 10 μm. (D) Western blot against total H3 and H3K27ac on histone extracts from control cells (DMSO) and cells treated with 0.5 and 3 μM TSA. White line indicates a lane was removed. (E) Gene expression distributions from biological replicates of control (DMSO, gray colors) and TSA-treated (red colors) cells. Notice the decrease in the fraction of non-expressed genes and the general shift in the distribution of TSA-treated cells. See also Figures S2 and S3.

Figure 4

Comparative Proximal Distribution of Chromatin Marks across Opisthokonta Species For each species, a plot shows the average normalized read coverage of four different histone modifications around the TSS (±5 kb), and heatmaps represent the same coverage for all genes sorted by level of expression. ChIP-seq data were obtained from publicly available datasets: Homo sapiens, Drosophila melanogaster, Caenorhabditis elegans, Nematostella vectensis, and Saccharomyces cerevisiae.

Figure 5

The Genomic Landscape of cis-Regulatory Elements in Capsaspora (A) Distribution of the number of regulatory sites per gene. (B) Number of cis-regulatory elements associated with different gene categories. Highlighted in gray are those with a significant enrichment (Wilcoxon rank-sum test p value < 0.01) compared with all genes. (C) Capsaspora transcription factor families sorted by the number of cis-regulatory elements associated per gene. (D) Preferential distribution of cis-regulatory sites across genomic features. (E) Pie charts showing the distribution of the number of stages in which each site is occupied (left) and the stage distribution of the stage-specific fraction of regulatory sites (right). (F) Heatmaps of clustered cis-regulatory elements (±2 kb) showing dynamic normalized ATAC nucleosome-free read coverage between stages. Plots show the associated average coverage profiles of each cluster. See also Figure S6.

Figure 6

Capsaspora Brachyury and Myc Regulation (A) Plot of ATAC-seq nucleosome-free reads average density around Bra motifs (top) and heatmap of the signal around the individual sites (bottom). (B) Differential distribution of regulatory sites containing Bra motif compared with all sites according to genomic feature (top) and stage/s in which the site is active (bottom). (C) Enrichment of different histone modifications (ChIP versus input) at Bra sites across stages. (D) Enrichment of different histone modifications (ChIP versus input) at Bra motifs in ATAC-defined sites compared with motifs occurring randomly in the genome. (E) Western blot of recombinant Capsaspora-Brachyury protein and Capsaspora nuclear protein extract, using Capsaspora-Brachyury affinity-purified antibody from guinea pig. (F) Capsaspora filopodial stage cell stained with phalloidin (red, actin cytoskeleton), DAPI (blue, nucleus), and Capsaspora-Brachyury antibody (green). Notice Bra localization in the nucleus. (G) Boxplot showing the Capsaspora-Brachyury ChIP-qPCR signal for predicted Bra regulatory sites versus random Bra motifs in the genome. (H) Illustrative case example of a predicted Bra regulatory site (highlighted in blue). For each feature, the top track corresponds to the filopodial stage, the middle track to the aggregative stage, and the bottom track to the cystic stage. Notice the decreased ATAC signal in the putative Bra-regulatory site in the cystic stage. (I) Enriched gene ontology (GO) terms and KEGG pathways among genes associated with Bra regulatory sites. (J) Enriched GO terms and KEGG pathways among genes associated with Bra regulatory sites with shared orthologs regulated by Bra in mouse. (K–O) Same as (A–D) and (I) for Capsaspora Myc. See also Figure S7.

Figure 7

Comparative Analysis of Regulatory Sites between Capsaspora and Animals (A) Distribution of ATAC-defined regulatory site sizes (bp) in Capsaspora, Drosophila and Homo. (B) Genomic feature distribution of regulatory sites in Capsaspora, Drosophila and Homo. (C) Enrichment of different histone modifications (ChIP versus input) at regulatory sites in distal (top) and proximal (bottom) intergenic regions in Capsaspora, Drosophila and Homo. In each boxplot, p values are indicated for Wilcoxon signed rank tests between H3K4me3 and H3K4me1 and between H3K4me3 and K27ac (only for distal intergenic). See also Figure S7.

Figure S1

Identification of Histone Modifications in Capsaspora, Related to Figure 1 (A) Histone N-terminal tail sequences of Capsaspora with all identified post-translational modifications and their location. A quotation mark indicates the impossibility of reliably assigning a modification to one or another of a pair of neighboring residues. (B) Representative MSMS analysis of modified peptides from H3 (from top to bottom): K4me3 (TK(me3)QTAR); K4me1 (TK(me)QTAR); K27ac (K(ac)TAVTSGGVKKPHR); K36me3 (KTAVTSGGVK(me3)KPHR). The b- and y-ion series are represented in red and blue, respectively. Non-fragmented precursor peptides are shown in green.

Figure S2

RNApolII ChIP-Seq Experiments, Related to Figure 2 (A) TSS-centered average normalized read coverage plots for RNApolII in the filopodial stage, using three different antibodies: 8WG16 (which preferentially recognizes unphosphorylated RNApolII), CTD4H8 (which recognizes both phospho- and unphosphorylated RNApolII) and S2P (which recognizes S2P-CTD phosphorylated RNApolII, the form associated to transcriptional elongation). The x axis spans −5Kb to +5Kb around the TSS. Shaded gray area represents the average size of Capsaspora genes. (B) Scatterplots of RNApolII coverage (log2 normalized reads) compared to mRNA expression levels (log2 FPKMs) in the filopodial stage. (C) Illustrative examples of RNApolII dynamic changes. Different genomic windows showing normalized coverage for different chromatin features. For each feature, the top track corresponds to filopodial stage, middle track to aggregative stage and bottom track to cystic stage.

Figure S3

Additional Illustrative Examples of Dynamic Chromatin Modifications in Capsaspora, Related to Figure 3 Different genomic windows showing normalized coverage for different chromatin features and their dynamic association with gene expression. For each feature, the top track corresponds to filopodial stage, middle track to aggregative stage and bottom track to cystic stage.

Figure S4

Capsaspora lincRNA Populations Defined by Chromatin Marks, Related to Figure 2 (A) Heatmap showing clustered lincRNA expression (RPKMs) across replicates of each stage. Only significantly differentially expressed lincRNAs (DESeq FDR < 0.05) are represented. (B) Characteristics of lincRNA loci compared with coding protein genes, including exon number distribution (top left), GC content (bottom left), length (top right), level of expression (middle right) and coefficient of variation in expression between stages and replicates (bottom right). (C) Heatmaps showing average read normalized coverage of four different histone modifications along lincRNA loci. (D) Characteristics of H3K4me1 versus H3K4me3 marked lincRNA loci. (E) Illustration of the genomic location of a lincRNA locus and normalized read coverage of histone modifications ChIP-seq. (F) RT-PCR validation of CUFF.777 lincRNA, revealing the existence of 3 isoforms. The minus sign indicates the negative control performed using RNA without reverse transcription to check for genomic DNA contamination.

Figure S5

RT-PCR Validation of lincRNAs, Related to Figure S6 Each panel shows the result of a PCR over poly-A-selected cDNA for each stage, as well as a control sample for genomic contamination (a cross-stages pool of the original non-retrotranscribed RNAs). The name of the lincRNA locus and the predicted size is indicated below each image. Cases of splicing-event validation are indicated schematically. An asterisk indicates that the observed sizes are smaller than predicted.

Figure S6

Examples of ATAC Profiling of Regulatory Sites, Related to Figure 5 (A and B) Different genomic windows (size indicated above) showing normalized coverage for different chromatin features and their dynamic association with gene expression. For each feature, the top track corresponds to filopodial stage, middle track to aggregative stage and bottom track to cystic stage. Significant peaks of ATAC nucleosome-free reads are highlighted.

Figure S7

Capsaspora Runx and NFAT/NFkappaB Regulatory Networks, Related to Figures 6 and 7 (A) Plot of ATAC-seq nucleosome-free average signal density around Runx motifs (top) and heatmap of signal around individual sites (bottom). (B) Differential distribution of regulatory sites containing Runx motif compared with all sites according to genomic feature (top) and stage/s where the site is active (bottom). (C) Enriched GO terms and KEGG pathways among genes associated with Runx regulatory sites. (D) Enrichment of different histone modifications (ChIP versus input) at Runx motifs in ATAC-defined sites compared with motifs occurring randomly in the genome (left) and at Runx sites across stages (right). (E) Plot of ATAC-seq nucleosome-free average signal density around NFAT/NFkappaB motifs (top) and heatmap of signal around individual sites (bottom). (F) Differential distribution of regulatory sites containing NFAT/NFkappaB motif compared with all sites according to genomic feature (top) and stage/s where the site is active (bottom). (G) Enriched GO terms and KEGG pathways among genes associated with NFAT/NFkappaB regulatory sites. (H) Enrichment of different histone modifications (ChIP versus input) at NFAT/NFkappaB motifs in ATAC-defined sites compare with motifs occurring randomly in the genome (left) and at NFAT/NFkappaB sites across stages (right). (I) Regulatory Site Distribution in Capsaspora and yeast. Genomic feature distribution of ATAC-defined regulatory sites in Capsaspora compared with the distribution of transcription factor binding sites (TFBS) extracted from UCSC (with evidence support > 2). In the left pie chart, TFBS at < 50bp of distance were merged.