Supporting data

LINE retrotransposons characterize mammalian tissue-specific and evolutionarily dynamic regulatory regions

Posted on by

Description

To investigate the mechanisms driving regulatory evolution across tissues, we experimentally mapped promoters, enhancers, and gene expression in liver, brain, muscle, and testis from ten diverse mammals. The regulatory landscape around genes included both tissue-shared and tissue-specific regulatory regions, where tissue-specific promoters and enhancers evolved most rapidly. Genomic regions switching between promoters and enhancers were more common across species, and less common across tissues within a single species. Long Interspersed Nuclear Elements (LINEs) played recurrent evolutionary roles: LINE L1s were associated with tissue-specific regulatory regions, whereas more ancient LINE L2s were associated with tissue-shared regulatory regions and with those switching between promoter and enhancer signatures across species. Our analyses of the tissue-specificity and evolutionary stability among promoters and enhancers reveal how specific LINE families have helped shape the dynamic mammalian regulome.This work was published in  Genome Biology .


First posted as a preprint to  BioRxiv .

Data access

The raw and processed high-throughput sequencing data are available in ArrayExpress. The ChIP-seq datasets have accession number  E-MTAB-7127 , and matched RNA-seq experiments  E-MTAB-8122. For the reannotation of genomes, additional RNA-seq dataset were produced and are also available under the accession number  E-MTAB-8118.

Scripts

The scripts used to run all the analyses are available as  Additional file 9: Data S1  associated with the  Genome Biology  manuscript. The script for projecting .bed file coordinates across species using Ensembl EPO whole genome alignments is included in the manuscript supplementary, but also available on  GitHub  as an example Ensembl Compara API script.

Final regulatory region calls per species

The results of postprocessing histone enrichment peaks to final regulatory region definitions across tissues is available  here . These files combine all per-tissue histone peaks calls within each species to define cross-tissue activity of active promoters, active enhancers and primed enhancers. For more details, see the Materials and methods section of the manuscript published in  Genome Biology .

File naming convention: 

{Species}_regRegions_allTissue_parsed.txt


The columns in the tab-delimited file correspond to:

Column number Description
1 Chromosome (Ensembl convention)
2 Start coordinate (Ensembl convention)
3 End coordinate (Ensembl convention)
4 Comma separated list of all regulatory region signatures at this location; multiple regulatory regions if different between tissues.
5 Comma separated list of all tissues with regulatory signatures
6 Unique IDs (see more details below)


Unique ID naming convention: 

{Species}_{Tissue}_{HistoneCombination}_{uniqueNumber}


For example: 

Mouse_Testis_H3K4me1-H3K27ac-H3K4me3_1


The unique IDs include the combination of histone peaks used to call regulatory regions within each species.
H3K27ac+H3K4me3 – active promoters, i.e. regions with both H3K27ac and H3K4me3
H3K27ac-H3K4me3+H3K4me1 – active enhancers, i.e. regions with H3K27ac and H3K4me1, but not H3K4me3
H3K4me1-H3K27ac-H3K4me3 – primed enhancers, i.e. region with only H3K4me1

All evolutionarily maintained regulatory regions

The results of finding all maintained regulatory regions – i.e. those regulatory regions in a species that align to a regulatory region in another species – are available  here . Briefly, these files were generated by taking all the final regulatory region calls per species (see above) and using Ensembl EPO alignments to check whether a regulatory region call in another species aligns with at least one base overlap. Every possible reciprocal pairwise comparison is represented by a seperate file. For the definition of evolutionarily maintained regulatory regions and more method details, see the Materials and methods section of the manuscript published in  Genome Biology .

File naming convention: 

{Species1}_regRegions_allTissue_mainRegs_to_{Species2}_active.txt

This file reports all regulatory regions from Species 1 that align to a regulatory region in Species 2, as well as the regulatory signature, tissue of activity and position, in both species.

The columns in the tab-delimited files correspond to:

Column number Description
1 Alignment of query region from Species 1 (see Column 4) to Species 2. Species 2 chromosome (Ensembl convention)
2 Alignment of query region from Species 1 (see Column 4) to Species 2. Species 2 start coordinate (Ensembl convention)
3 Alignment of query region from Species 1 (see Column 4) to Species 2. Species 2 end coordinate (Ensembl convention)
4 Species 1 query region used to align to Species 2. Syntax Chromsome:Start-End (Ensembl convention)
5 Comma separated list of regulatory signature(s) in Species 1
6 Comma separated list of all tissues with regulatory signatures in Species 1
7 Unique IDs in Species 1 (see more details above under Final regulatory region calls)
8 Species 2 regulatory region chromosome (Ensembl convention)
9 Species 2 regulatory region start coordinate (Ensembl convention)
10 Species 2 regulatory region end coordinate (Ensembl convention)
11 Comma separated list of regulatory signature(s) in Species 2
12 Comma separated list of all tissues with regulatory signatures in Species 2
13 Unique IDs in Species 2 (see more details above under Final regulatory region calls)

Final normalised RNA-seq data

The results of RNA-seq normalisation across tissues and filtering is available  here . Briefly, the cufflinks suite was used to normalise RNA-seq across Ensembl annotated genes and transcripts and lowly covered genes/transcripts were removed using an FPKM cutoff. For more details, see the Materials and methods section of the manuscript published in  Genome Biology .  File naming convention:  {Species}_{genes/isoforms}.fpkm_table_FPKMfilter
Genes files correspond to all Ensembl annotated genes.
Isoform files correspond to all Ensembl annotated transcripts.

Repeat masking of genomes

Repeats were masked with RepeatMasker and the RepBase database, for more information please see the methods section published in  Genome Biology . The gtf formatted results of repeat masking are available  here .

Identification of male heterogametic sex-determining regions on the Atlantic herring Clupea harengus genome

Posted on by

bibtex key=32293027]

Description

The sex determination system of Atlantic herring  Clupea harengus  L., a commercially important fish, was investigated. Low coverage whole-genome sequencing of 48 females and 55 males and a genome-wide association study revealed two regions on chromosomes 8 and 21 associated with sex. The genotyping data of the single nucleotide polymorphisms associated with sex showed that 99.4% of the available female genotypes were homozygous, whereas 68.6% of the available male genotypes were heterozygous. This is close to the theoretical expectation of homo/heterozygous distribution at low sequencing coverage when the males are factually heterozygous. This suggested a male heterogametic sex determination system in  C. harengus , consistent with other species within the Clupeiformes group. There were 76 protein coding genes on the sex regions but none of these genes were previously reported master sex regulation genes, or obviously related to sex determination. However, many of these genes are expressed in testis or ovary in other species, but the exact genes controlling sex determination in  C. harengus  could not be identified.

Full details are provided in our  open access publication in the Journal of Fish Biology

Data access

The sequencing reads were submitted to the European Nucleotide Archive repository and are available in 103 samples with accession numbers from  ERS4329014  to  ERS4329116 .

SNPs were called using FreeBayes v1.1.0 and are available  here .

The results from the GWAS analysis can be found  here .

Clustered CTCF binding is an evolutionary mechanism to maintain topologically associating domains

Posted on by
  • E Kentepozidou, SJ Aitken, C Feig, K Stefflova, X Ibarra-Soria, DT Odom, M Roller, P Flicek. Clustered CTCF binding is an evolutionary mechanism to maintain topologically associating domains. Genome Biol 2020;21(1):5. doi:10.1186/s13059-019-1894-x
    [BibTeX] [Abstract]

    BACKGROUND: CTCF binding contributes to the establishment of a higher-order genome structure by demarcating the boundaries of large-scale topologically associating domains (TADs). However, despite the importance and conservation of TADs, the role of CTCF binding in their evolution and stability remains elusive. RESULTS: We carry out an experimental and computational study that exploits the natural genetic variation across five closely related species to assess how CTCF binding patterns stably fixed by evolution in each species contribute to the establishment and evolutionary dynamics of TAD boundaries. We perform CTCF ChIP-seq in multiple mouse species to create genome-wide binding profiles and associate them with TAD boundaries. Our analyses reveal that CTCF binding is maintained at TAD boundaries by a balance of selective constraints and dynamic evolutionary processes. Regardless of their conservation across species, CTCF binding sites at TAD boundaries are subject to stronger sequence and functional constraints compared to other CTCF sites. TAD boundaries frequently harbor dynamically evolving clusters containing both evolutionarily old and young CTCF sites as a result of the repeated acquisition of new species-specific sites close to conserved ones. The overwhelming majority of clustered CTCF sites colocalize with cohesin and are significantly closer to gene transcription start sites than nonclustered CTCF sites, suggesting that CTCF clusters particularly contribute to cohesin stabilization and transcriptional regulation. CONCLUSIONS: Dynamic conservation of CTCF site clusters is an apparently important feature of CTCF binding evolution that is critical to the functional stability of a higher-order chromatin structure.

    @Article{31910870,
author = {Kentepozidou E and Aitken SJ and Feig C and Stefflova K and Ibarra-Soria X and Odom DT and Roller M and Flicek P},
title = {Clustered CTCF binding is an evolutionary mechanism to maintain topologically associating domains},
journal = {Genome Biol},
volume = {21},
number = {1},
pages = {5},
year = {2020},
doi = {10.1186/s13059-019-1894-x},
note = {First posted as a preprint: 12 June 2019},
abstract = {BACKGROUND: CTCF binding contributes to the establishment of a higher-order genome structure by demarcating the boundaries of large-scale topologically associating domains (TADs). However, despite the importance and conservation of TADs, the role of CTCF binding in their evolution and stability remains elusive. RESULTS: We carry out an experimental and computational study that exploits the natural genetic variation across five closely related species to assess how CTCF binding patterns stably fixed by evolution in each species contribute to the establishment and evolutionary dynamics of TAD boundaries. We perform CTCF ChIP-seq in multiple mouse species to create genome-wide binding profiles and associate them with TAD boundaries. Our analyses reveal that CTCF binding is maintained at TAD boundaries by a balance of selective constraints and dynamic evolutionary processes. Regardless of their conservation across species, CTCF binding sites at TAD boundaries are subject to stronger sequence and functional constraints compared to other CTCF sites. TAD boundaries frequently harbor dynamically evolving clusters containing both evolutionarily old and young CTCF sites as a result of the repeated acquisition of new species-specific sites close to conserved ones. The overwhelming majority of clustered CTCF sites colocalize with cohesin and are significantly closer to gene transcription start sites than nonclustered CTCF sites, suggesting that CTCF clusters particularly contribute to cohesin stabilization and transcriptional regulation. CONCLUSIONS: Dynamic conservation of CTCF site clusters is an apparently important feature of CTCF binding evolution that is critical to the functional stability of a higher-order chromatin structure.},}

Description

CTCF binding contributes to the establishment of a higher-order genome structure by demarcating the boundaries of large-scale topologically associating domains (TADs). However, despite the importance and conservation of TADs, the role of CTCF binding in their evolution and stability remains elusive.

We carry out an experimental and computational study that exploits the natural genetic variation across five closely related species to assess how CTCF binding patterns stably fixed by evolution in each species contribute to the establishment and evolutionary dynamics of TAD boundaries. We perform CTCF ChIP-seq in multiple mouse species to create genome-wide binding profiles and associate them with TAD boundaries. Our analyses reveal that CTCF binding is maintained at TAD boundaries by a balance of selective constraints and dynamic evolutionary processes. Regardless of their conservation across species, CTCF binding sites at TAD boundaries are subject to stronger sequence and functional constraints compared to other CTCF sites. TAD boundaries frequently harbor dynamically evolving clusters containing both evolutionarily old and young CTCF sites as a result of the repeated acquisition of new species-specific sites close to conserved ones. The overwhelming majority of clustered CTCF sites colocalize with cohesin and are significantly closer to gene transcription start sites than nonclustered CTCF sites, suggesting that CTCF clusters particularly contribute to cohesin stabilization and transcriptional regulation.

Dynamic conservation of CTCF site clusters is an apparently important feature of CTCF binding evolution that is critical to the functional stability of a higher-order chromatin structure.

Full details are available in the manuscript published in  Genome Biology  and the preprint submitted to BioRxiv .

Raw Data

All ChIP-seq and RNA-seq data generated in this study are available in the Array Express repository under the accession numbers E-MTAB-8014 , E-MTAB-8471 , and E-MTAB-8016 . Additional ChIP-seq and RNA-seq data that were reused in the study are available under the accession numbers E-MTAB-5769 , E-MTAB-1091 and E-MTAB-2483 .

Hi-C-derived TADs were retrieved from Table S1. Genomic Regions of Interest in: Vietri Rudan M, Barrington C, Henderson S, Ernst C, Odom DT, Tanay A, et al. Comparative Hi-C reveals that CTCF underlies evolution of chromosomal domain architecture . Cell Rep. 2015;10:1297–309 .

Transposable elements were retrieved from the supplementary page accompanying Thybert D, Roller M, Navarro FCP, Fiddes I, Streeter I, Feig C, et al. Repeat associated mechanisms of genome evolution and function revealed by the Mus caroli and Mus pahari genomes. Genome Res. 2018;28:448–59.

Conserved CTCF binding sites 

The results of the cross-species comparisons of CTCF binding conservation, centred on the mm10 genome, can be found here . Column descriptions:

CBS_ID Unique ID
coord_on_mmus Genome coordinates on the mm10 genome
m_musculus, m_castaneus, m_spretus, m_caroli, m_pahari These columns take values [0|1|2].
0 –> The site does not have any orthologous alignment on the genome of the corresponding species.
1–> the site has an orthologous alignment on the genome of the corresponding species, but is not bound by CTCF.
2 –> The site has an orthologous alignment on the genome of the corresponding species and is also bound by CTCF in the corresponding species (validated with ChIP-seq data from the corresponding species).
Cons_Score Number of species the site is conserved in (has orthologous alignment AND is bound by CTCF)
Align_Score Number of species with orthologous alignment but is NOT bound by CTCF

To find the genome coordinates of CTCF sites projected to all of the mouse species, a similar file exists here . It has all the columns of the above table, and additionally also coordinate columns  corresponding to the orthologous alignment of every site on the genome of each of the other species ( mmus, mcas, mspr, mcar  and  mpah ).

Using long and linked reads to improve an Atlantic herring (Clupea harengus) genome assembly

Posted on by
  • S Kongsstovu Í, SO Mikalsen, EÍ Homrum, JA Jacobsen, P Flicek, HA Dahl. Using long and linked reads to improve an Atlantic herring (Clupea harengus) genome assembly. Sci Rep 2019;9(1):17716. doi:10.1038/s41598-019-54151-9
    [BibTeX] [Abstract]

    Atlantic herring (Clupea harengus) is one of the most abundant fish species in the world. It is an important economical and nutritional resource, as well as a crucial part of the North Atlantic ecosystem. In 2016, a draft herring genome assembly was published. Being a species of such importance, we sought to independently verify and potentially improve the herring genome assembly. We sequenced the herring genome generating paired-end, mate-pair, linked and long reads. Three assembly versions of the herring genome were generated based on a de novo assembly (A1), which was scaffolded using linked and long reads (A2) and then merged with the previously published assembly (A3). The resulting assemblies were compared using parameters describing the size, fragmentation, correctness, and completeness of the assemblies. Results showed that the A2 assembly was less fragmented, more complete and more correct than A1. A3 showed improvement in fragmentation and correctness compared with A2 and the published assembly but was slightly less complete than the published assembly. Thus, we here confirmed the previously published herring assembly, and made improvements by further scaffolding the assembly and removing low-quality sequences using linked and long reads and merging of assemblies.

    @Article{31776409,
author = {Í Kongsstovu S and Mikalsen SO and Homrum EÍ and Jacobsen JA and Flicek P and Dahl HA},
title = {Using long and linked reads to improve an Atlantic herring (Clupea harengus) genome assembly},
journal = {Sci Rep},
volume = {9},
number = {1},
pages = {17716},
year = {2019},
doi = {10.1038/s41598-019-54151-9},
abstract = {Atlantic herring (Clupea harengus) is one of the most abundant fish species in the world. It is an important economical and nutritional resource, as well as a crucial part of the North Atlantic ecosystem. In 2016, a draft herring genome assembly was published. Being a species of such importance, we sought to independently verify and potentially improve the herring genome assembly. We sequenced the herring genome generating paired-end, mate-pair, linked and long reads. Three assembly versions of the herring genome were generated based on a de novo assembly (A1), which was scaffolded using linked and long reads (A2) and then merged with the previously published assembly (A3). The resulting assemblies were compared using parameters describing the size, fragmentation, correctness, and completeness of the assemblies. Results showed that the A2 assembly was less fragmented, more complete and more correct than A1. A3 showed improvement in fragmentation and correctness compared with A2 and the published assembly but was slightly less complete than the published assembly. Thus, we here confirmed the previously published herring assembly, and made improvements by further scaffolding the assembly and removing low-quality sequences using linked and long reads and merging of assemblies.},}

Description

Atlantic herring (Clupea harengus) is one of the most abundant fish species in the world. It is an important economical and nutritional resource, as well as a crucial part of the North Atlantic ecosystem. In 2016, a draft herring genome assembly was published. Being a species of such importance, we sought to independently verify and potentially improve the herring genome assembly. We sequenced the herring genome generating paired-end, mate-pair, linked and long reads. Three assembly versions of the herring genome were generated based on a de novo assembly (A1), which was scaffolded using linked and long reads (A2) and then merged with the previously published assembly (A3). The resulting assemblies were compared using parameters describing the size, fragmentation, correctness, and completeness of the assemblies. Results showed that the A2 assembly was less fragmented, more complete and more correct than A1. A3 showed improvement in fragmentation and correctness compared with A2 and the published assembly but was slightly less complete than the published assembly. Thus, we here confirmed the previously published herring assembly, and made improvements by further scaffolding the assembly and removing low-quality sequences using linked and long reads and merging of assemblies.

Full details are provided in our  open access publication in Scientific Reports

Data access

The genome assemblies and sequencing reads were submitted to the European Nucleotide Archive and are available with the following accession numbers:

Assembly Accession number
A1 GCA_902175115.1
A2 GCA_902175115.2
A3 GCA_902175115.3
Sequencing reads Accession number
Paired-end library ERR2853087
Mate-pair 4500 library ERR2853088 ,  ERR2853089 ,  ERR2853090
Mate-pair 7000 library ERR2853091 ,  ERR2853092 ,  ERR2853093
Linked (10x Genomics) reads ERR2853094
MinION reads ERR2809166 ,  ERR2809167 ,  ERR2809168 ,  ERR2809169

Tables 2 and 6 show the QUAST results that were of most interest or relevance. The full QUAST results are available for the  reference  and  no reference  cases.

For connexin analysis the connexin sequences were aligned to the assemblies. Alignments are available for the  A1 ,  A2 ,  A3  and  draft  assemblies.

The features identified using FRC bam  are available here in gff format for the  A1 ,  A2  and  A3  assemblies. All reference data features are available  here .

Whole genome alignments were generated using the web tool D-Genies. Only the alignment between A3 and the draft assembly (Figure 2) are represented in the article. The files are available for alignments between  A1 and A2 ,  A2 and A3 , and  A3 and the draft assemblies .

Repeat associated mechanisms of genome evolution and function revealed by the Mus caroli and Mus pahari genomes

Posted on by
  • D Thybert, M Roller, FCP Navarro, I Fiddes, I Streeter, C Feig, D Martin-Galvez, M Kolmogorov, V Janoušek, W Akanni, B Aken, S Aldridge, V Chakrapani, W Chow, L Clarke, C Cummins, A Doran, M Dunn, L Goodstadt, K Howe, M Howell, AA Josselin, RC Karn, CM Laukaitis, L Jingtao, F Martin, M Muffato, S Nachtweide, MA Quail, C Sisu, M Stanke, K Stefflova, C Van Oosterhout, F Veyrunes, B Ward, F Yang, G Yazdanifar, A Zadissa, DJ Adams, A Brazma, M Gerstein, B Paten, S Pham, TM Keane, DT Odom, P Flicek. Repeat associated mechanisms of genome evolution and function revealed by the Mus caroli and Mus pahari genomes. Genome Res 2018;28(4):448–459. doi:10.1101/gr.234096.117
    [BibTeX] [Abstract]

    Understanding the mechanisms driving lineage-specific evolution in both primates and rodents has been hindered by the lack of sister clades with a similar phylogenetic structure having high-quality genome assemblies. Here, we have created chromosome-level assemblies of the Mus caroli and Mus pahari genomes. Together with the Mus musculus and Rattus norvegicus genomes, this set of rodent genomes is similar in divergence times to the Hominidae (human-chimpanzee-gorilla-orangutan). By comparing the evolutionary dynamics between the Muridae and Hominidae, we identified punctate events of chromosome reshuffling that shaped the ancestral karyotype of Mus musculus and Mus caroli between 3 and 6 million yr ago, but that are absent in the Hominidae. Hominidae show between four- and sevenfold lower rates of nucleotide change and feature turnover in both neutral and functional sequences, suggesting an underlying coherence to the Muridae acceleration. Our system of matched, high-quality genome assemblies revealed how specific classes of repeats can play lineage-specific roles in related species. Recent LINE activity has remodeled protein-coding loci to a greater extent across the Muridae than the Hominidae, with functional consequences at the species level such as reproductive isolation. Furthermore, we charted a Muridae-specific retrotransposon expansion at unprecedented resolution, revealing how a single nucleotide mutation transformed a specific SINE element into an active CTCF binding site carrier specifically in Mus caroli, which resulted in thousands of novel, species-specific CTCF binding sites. Our results show that the comparison of matched phylogenetic sets of genomes will be an increasingly powerful strategy for understanding mammalian biology.

    @Article{29563166,
author = {Thybert D and Roller M and Navarro FCP and Fiddes I and Streeter I and Feig C and Martin-Galvez D and Kolmogorov M and Janoušek V and Akanni W and Aken B and Aldridge S and Chakrapani V and Chow W and Clarke L and Cummins C and Doran A and Dunn M and Goodstadt L and Howe K and Howell M and Josselin AA and Karn RC and Laukaitis CM and Jingtao L and Martin F and Muffato M and Nachtweide S and Quail MA and Sisu C and Stanke M and Stefflova K and Van Oosterhout C and Veyrunes F and Ward B and Yang F and Yazdanifar G and Zadissa A and Adams DJ and Brazma A and Gerstein M and Paten B and Pham S and Keane TM and Odom DT and Flicek P},
title = {Repeat associated mechanisms of genome evolution and function revealed by the Mus caroli and Mus pahari genomes},
journal = {Genome Res},
volume = {28},
number = {4},
pages = {448--459},
year = {2018},
doi = {10.1101/gr.234096.117},
howpublished = {Advanced online publication: 21 March 2018},
note = {First posted as a preprint: 2 July 2017},
abstract = {Understanding the mechanisms driving lineage-specific evolution in both primates and rodents has been hindered by the lack of sister clades with a similar phylogenetic structure having high-quality genome assemblies. Here, we have created chromosome-level assemblies of the Mus caroli and Mus pahari genomes. Together with the Mus musculus and Rattus norvegicus genomes, this set of rodent genomes is similar in divergence times to the Hominidae (human-chimpanzee-gorilla-orangutan). By comparing the evolutionary dynamics between the Muridae and Hominidae, we identified punctate events of chromosome reshuffling that shaped the ancestral karyotype of Mus musculus and Mus caroli between 3 and 6 million yr ago, but that are absent in the Hominidae. Hominidae show between four- and sevenfold lower rates of nucleotide change and feature turnover in both neutral and functional sequences, suggesting an underlying coherence to the Muridae acceleration. Our system of matched, high-quality genome assemblies revealed how specific classes of repeats can play lineage-specific roles in related species. Recent LINE activity has remodeled protein-coding loci to a greater extent across the Muridae than the Hominidae, with functional consequences at the species level such as reproductive isolation. Furthermore, we charted a Muridae-specific retrotransposon expansion at unprecedented resolution, revealing how a single nucleotide mutation transformed a specific SINE element into an active CTCF binding site carrier specifically in Mus caroli, which resulted in thousands of novel, species-specific CTCF binding sites. Our results show that the comparison of matched phylogenetic sets of genomes will be an increasingly powerful strategy for understanding mammalian biology.},}

Description

Understanding the mechanisms driving lineage-specific evolution in both primates and rodents has been hindered by the lack of sister clades with a similar phylogenetic structure having high-quality genome assemblies. Here, we have created chromosome-level assemblies of the  Mus caroli  and  Mus pahari  genomes. Together with the  Mus musculus  and  Rattus norvegicus  genomes, this set of rodent genomes is similar in divergence times to the Hominidae (human-chimpanzee-gorilla-orangutan). By comparing the evolutionary dynamics between the Muridae and Hominidae, we identified punctate events of chromosome reshuffling that shaped the ancestral karyotype of  Mus musculus  and  Mus caroli  between 3 to 6 MYA, but that are absent in the Hominidae. In fact, Hominidae show between four- and seven-fold lower rates of nucleotide change and feature turnover in both neutral and functional sequences suggesting an underlying coherence to the Muridae acceleration. Our system of matched, high-quality genome assemblies revealed how specific classes of repeats can play lineage-specific roles in related species. For example, recent LINE activity has remodeled protein-coding loci to a greater extent across the Muridae than the Hominidae, with functional consequences at the species level such as reproductive isolation. Furthermore, we charted a Muridae-specific retrotransposon expansion at unprecedented resolution, revealing how a single nucleotide mutation transformed a specific SINE element into an active CTCF binding site carrier specifically in  Mus caroli . This process resulted in thousands of novel, species-specific CTCF binding sites. Our results demonstrate that the comparison of matched phylogenetic sets of genomes will be an increasingly powerful strategy for understanding mammalian biology.

Full details are provided in our open access publication in  Genome Research . 

Data access

The genome assemblies of Mus caroli and Mus pahari were submitted to the European Nucleotide Archive and are available with accession numbers  GCA_900094665  for  Mus caroli  and  GCA_900095145  for  Mus pahari . All reads from the ChIP-seq and RNA-seq experiments in this study were submitted to ArrayExpress and are available with accession numbers  E-MTAB-5768  (RNA-seq) and  E-MTAB-5769  (ChIP-seq).

Transposible element annotation

The results of RepeatMasker identifications were postprocessed to merge fragmented hits and remove non-transposable repeats to create the final set used in this study. The results of postprocessing are available  here . The columns in the postprocessed files correspond to:

Column number Description
1 Toplevel genomic segment, i.e. chromosome, scaffold
2 Start position of match in genomic segment
3 End position of match in genomic segment
4 RepeatMasker result: % substitutions in matching region compared to the consensus
5 RepeatMasker result: % of bases opposite a gap in the query sequence (deleted bp)
6 RepeatMasker result: % of bases opposite a gap in the repeat consensus (inserted bp)
7 Transposable element classification: transposable element class
8 Transposable element classification: transposable element family
9 Transposable element classification: transposable element subfamily
10 Unique id created for this study

Whole genome alignments

The whole genome alignments used in this study are available  here  in Ensembl Multi Format (EMF) and multiple alignment format (MAF). Pairwise whole genome alignments were generated using LastZ and multiple whole genome alignments with the Enredo-Pecan-Ortheus (EPO) pipeline. For more information please refer to the Methods in our  paper .

CTCF occupancy sites

The peak sets per biological replicate are available in ArrayExpress under the accession number  E-MTAB-5769 . The peaks present in at least two biological replicates are available  here . The B2_Mm1 transposable elements used as input for building the neighbour joining tree are available  here . For more information please refer to the Methods in our  paper .

CTCF occupancy sites associated with repetitive elements are available  here . The columns in the postprocessed files describe:

Column name Description
PeakChr Toplevel genomic segment containing the CTCF peak, i.e. chromosome, scaffold
PeakStart Start position of CTCF peak in genomic segment
PeakEnd End position of CTCF peak in genomic segment
RepeatChr Toplevel genomic segment containing the repeat element, i.e. chromosome, scaffold consensus
RepeatStart Start position of repeat element in genomic segment
RepeatEnd End position of repeat element in genomic segment
RepeatClass Transposable element classification: transposable element class
RepeatFamily Transposable element classification: transposable element family
RepeatSubfamily Transposable element classification: transposable element subfamily

Results of multiple alignments of CTCF sites between the rodents are available  here . The name of the file corresponds to the species which was used as query to align to other species. The first column contains all CTCF sites from the query species, and each following column an alignment to another species. If a column contains genomic coordinates, a bound CTCF site is aligned. A lack of aligned CTCF sites is marked with an “X”. A row with “X” in all columns except the first represents a species-specific CTCF binding site. The position of peaks is in “Chr:Start-End” format.

Mus pahari breakpoints

chr1: mm_7:27,524,252-end + mm_19:start-end
                chr2: mm_5:start-30,341,378 + mm_6:strt-end
                chr3: mm_2:23,133,000-end
                chr4: mm_3:start-end
                chr5: mm_1:18,015,577-end
                chr6: mm_4:45,383,207-end
                chr7: mm_12:start-end
                chr8: mm_14:start-end
                chr9: mm_10:33,415,951-end
                chr10: mm_9:start-end
                chr11: mm_13:[66,989,676/67,139,392]-end + mm_15:start-32,682,454
                chr12: mm_16:start-end
                chr13: mm_11:start-31,004,372+mm_5:33,090,849-[109,638,275/110,356,165]
                chr14: mm_11:31,014,735-end
                chr15: mm_18:start-end
                chr16: mm_2:start-23,116,788 + mm_13:start-[66,989,676/67,139,392]
                chr17: mm_15:32,682,454-end
                chr18: mm_17:31,811,330-end
                chr19: mm_7:start-27,524,252+ mm_8:start-75,294,883
                chr20: mm_8:75,299,631-end
                chr21: mm_17:start-27,052,643 + mm_10:start-33,381,753
                chr22: mm_1:start-17,966,305 + mm_4:start-45,383,207-end
                chr23: mm_5:[109,638,275/110,356,165]-end

The above lines refer to the composition of the Music pahari chromosomes. For example, chr1 of Mus pahari is composed of chr7 of Mus musculus starting from the position 27,524,252 until the end and then chr19 of Mus musculus from start to end. This gives a break point from the Mus musculus perspective.

For data points such as [66,989,676/67,139,392], this results from the array CGH being done twice. We have the mean between these two values for the analysis of repeat enrichment in the paper.

Mouse-rat breakpoints

chr1 rn_5:start-13Mb + rn_9:24-104Mb + rn_13:start-end
                chr2:rn_17:70Mb-end + rn_3:start-end
                chr3:rn_2:87Mb-end
                chr4:rn_5:16Mb- end
                chr5:rn_4:start-28Mb + rn_14:start-83Mb + rn_chr12:start-end
                chr6: rn_chr4:28Mb-end
                chr7:rn_1:63Mb-218Mb
                chr8:rn_16:19Mb-end+ rn_19:start-end
                chr9:rn_8:start-end
                chr10:rn_1:start-44Mb+ rn_20:10Mb-end+rn_7:start-67Mb
                chr11:rn_14:83Mb-end + rn_10:15Mb-end
                chr12:rn_6:27Mb-end
                chr13:rn_17:start-70Mb + rn_2:start-51Mb
                chr14:rn_15:start-23Mb + rn_16:start-18Mb + rn_15:23Mb-end
                chr15:rn_2:52Mb-86Mb + rn_7:71Mb-end
                chr16:rn_10:start-12Mb + rn_11:start-end
                chr17:rn_1:44Mb-62Mb + rn_20:start-10Mb +  rn_9:start-24Mb +rn_9:104Mb-end + rn_6:start-25Mb
                chr18:rn_17:54Mb-62Mb + rn_18:start-end
                chr19:rn_1:218Mb-end

The mouse rat synteny breaks have been defined using the  Ensembl synteny tool . In order to match the chromosome painting resolution we only retained breaks involved in inter-chromosome rearrangement involved in rearrangement with a size of 5Mb or more.

Complexity and conservation of regulatory landscapes underlie evolutionary resilience of mammalian gene expression

Posted on by
  • C Berthelot, D Villar, JE Horvath, DT Odom, P Flicek. Complexity and conservation of regulatory landscapes underlie evolutionary resilience of mammalian gene expression. Nat Ecol Evol 2018;2(1):152–163. doi:10.1038/s41559-017-0377-2
    [BibTeX] [Abstract]

    To gain insight into how mammalian gene expression is controlled by rapidly evolving regulatory elements, we jointly analysed promoter and enhancer activity with downstream transcription levels in liver samples from 15 species. Genes associated with complex regulatory landscapes generally exhibit high expression levels that remain evolutionarily stable. While the number of regulatory elements is the key driver of transcriptional output and resilience, regulatory conservation matters: elements active across mammals most effectively stabilize gene expression. In contrast, recently evolved enhancers typically contribute weakly, consistent with their high evolutionary plasticity. These effects are observed across the entire mammalian clade and are robust to potential confounders, such as the gene expression level. Using liver as a representative somatic tissue, our results illuminate how the evolutionary stability of gene expression is profoundly entwined with both the number and conservation of surrounding promoters and enhancers.

    @Article{29180706,
author = {Berthelot C and Villar D and Horvath JE and Odom DT and Flicek P},
title = {Complexity and conservation of regulatory landscapes underlie evolutionary resilience of mammalian gene expression},
journal = {Nat Ecol Evol},
volume = {2},
number = {1},
pages = {152--163},
year = {2018},
doi = {10.1038/s41559-017-0377-2},
howpublished = {Advanced online publication: 27 November 2017},
note = {First posted as a preprint: 7 April 2017},
abstract = {To gain insight into how mammalian gene expression is controlled by rapidly evolving regulatory elements, we jointly analysed promoter and enhancer activity with downstream transcription levels in liver samples from 15 species. Genes associated with complex regulatory landscapes generally exhibit high expression levels that remain evolutionarily stable. While the number of regulatory elements is the key driver of transcriptional output and resilience, regulatory conservation matters: elements active across mammals most effectively stabilize gene expression. In contrast, recently evolved enhancers typically contribute weakly, consistent with their high evolutionary plasticity. These effects are observed across the entire mammalian clade and are robust to potential confounders, such as the gene expression level. Using liver as a representative somatic tissue, our results illuminate how the evolutionary stability of gene expression is profoundly entwined with both the number and conservation of surrounding promoters and enhancers.},}

Description

To gain insight into how mammalian gene expression is controlled by rapidly evolving regulatory elements, we jointly analysed promoter and enhancer activity with downstream transcription levels in liver samples from fifteen species. Genes associated with complex regulatory landscapes generally exhibit high expression levels that remain evolutionarily stable. While the number of regulatory elements is the key driver of transcriptional output and resilience, regulatory conservation matters: elements active across mammals most effectively stabilise gene expression. In contrast, recently-evolved enhancers typically contribute weakly, consistent with their high evolutionary plasticity. These effects are observed across the entire mammalian clade and robust to potential confounders, such as gene expression level. Using liver as a representative somatic tissue, our results illuminate how the evolutionary stability of gene expression is profoundly entwined with both the number and conservation of surrounding promoters and enhancers.

Full details are available in our paper published in  Nature Ecology and Evolution  with full text freely available at  EuropePMC .

Raw Data

The raw RNA-seq data from livers of 25 mammalian species can be found in ArrayExpress with the accession number  E-MTAB-4550 , with the exception of three human and four mouse datasets, previously reported in  E-MTAB-4052 .

The processed RNA-seq datasets are also available from ArrayExpress with   E-MTAB-4550  (after read alignment with TopHat2 and transcript quantifications with Cufflinks – please refer to Methods in the  bioRxiv preprint ).

Average gene expression levels per species

The gene expression summaries for all replicates in each species are accessible  here . Expression levels are provided both as FPKM (as output by Cufflinks) and after TPM transform for each replicate.

Promoter and enhancer datasets

Please refer to the companion webpage for  Villar et al. 2015  for the raw datasets, consensus ChIP-seq peaks and evolutionary conservation analyses.

Putative target genes of the active promoters and enhancers in each species are accessible  here  (please refer to the Methods in the  bioRxiv preprint  for details on the putative target assignation).

Gene expression levels across species

The set of orthologous genes used in this study is accessible  here .

The table comparing normalized average expression levels across species is accessible  here .
This table additionally includes the following information:

orthtype Whether the gene is a strict 1-to-1 ortholog across all study species (0=no, 1=yes)
meanexp Mean expression across species
stdexp Standard deviation of expression across species
cvexp Coefficient of variation across species
core Whether the gene is a core liver gene (0=no, 1=yes)
hk Whether the gene is a housekeeping gene (0=no, 1=yes)
cvstab Whether the gene was classified as stable, variable, not expressed or unmatched (see Methods).

Meta-genes and their regulatory landscapes

The integrated summary of the regulatory landscape over 20 species (meta-promoters and meta-enhancers) is accessible  here .
Important:  This file is an overview of cross-alignable regulatory regions based on whole-genome alignments, and their putative target genes. It was designed to investigate global trends (rather than individual loci) and may be locally affected by alignment anomalies. Please pay careful attention to the ‘flag’ column (see below) if interested in specific loci.

This table includes the following information:

id The species identifiers of the active regulatory elements that form the meta-element
activeAll Number of species where the meta-element is active (out of all twenty)
activeRefs Number of reference species where the meta-element is active (out of ten reference species)
sequenceFound Number of reference species with an identifiable orthologous sequence (out of ten reference species)
averageLength Average length of the meta-element across the species where it is active
gene Putative target gene (Ensembl human Gene ID)
flag “Warning” indicates instances of regulatory elements that are included in the meta-element by sequence alignment, but are not in the vicinity of the putative target gene. Most of the time this is due to assembly fragmentation – but in some instances this may result from erroneous alignment of paralogous sequences. For these cases, the species identifiers of affected regulatory elements are given (e.g. Warning: sarHarH3K4me33870).
type Whether the meta-element is annotated as a promoter or an enhancer, based on its histone marking across species

Interplay of cis and trans mechanisms driving transcription factor binding and gene expression evolution

Posted on by
  • ES Wong, BM Schmitt, A Kazachenka, D Thybert, A Redmond, F Connor, TF Rayner, C Feig, AC Ferguson-Smith, JC Marioni, DT Odom, P Flicek. Interplay of cis and trans mechanisms driving transcription factor binding and gene expression evolution. Nat Commun 2017;8(1):1092. doi:10.1038/s41467-017-01037-x
    [BibTeX] [Abstract]

    Noncoding regulatory variants play a central role in the genetics of human diseases and in evolution. Here we measure allele-specific transcription factor binding occupancy of three liver-specific transcription factors between crosses of two inbred mouse strains to elucidate the regulatory mechanisms underlying transcription factor binding variations in mammals. Our results highlight the pre-eminence of cis-acting variants on transcription factor occupancy divergence. Transcription factor binding differences linked to cis-acting variants generally exhibit additive inheritance, while those linked to trans-acting variants are most often dominantly inherited. Cis-acting variants lead to local coordination of transcription factor occupancies that decay with distance; distal coordination is also observed and may be modulated by long-range chromatin contacts. Our results reveal the regulatory mechanisms that interplay to drive transcription factor occupancy, chromatin state, and gene expression in complex mammalian cell states.

    @Article{29061983,
author = {Wong ES and Schmitt BM and Kazachenka A and Thybert D and Redmond A and Connor F and Rayner TF and Feig C and Ferguson-Smith AC and Marioni JC and Odom DT and Flicek P},
title = {Interplay of cis and trans mechanisms driving transcription factor binding and gene expression evolution},
journal = {Nat Commun},
volume = {8},
number = {1},
pages = {1092},
year = {2017},
doi = {10.1038/s41467-017-01037-x},
note = {First posted as a preprint: 19 June 2016},
abstract = {Noncoding regulatory variants play a central role in the genetics of human diseases and in evolution. Here we measure allele-specific transcription factor binding occupancy of three liver-specific transcription factors between crosses of two inbred mouse strains to elucidate the regulatory mechanisms underlying transcription factor binding variations in mammals. Our results highlight the pre-eminence of cis-acting variants on transcription factor occupancy divergence. Transcription factor binding differences linked to cis-acting variants generally exhibit additive inheritance, while those linked to trans-acting variants are most often dominantly inherited. Cis-acting variants lead to local coordination of transcription factor occupancies that decay with distance; distal coordination is also observed and may be modulated by long-range chromatin contacts. Our results reveal the regulatory mechanisms that interplay to drive transcription factor occupancy, chromatin state, and gene expression in complex mammalian cell states.},}

Description 

Noncoding regulatory variants play a central role in the genetics of human diseases and in evolution. Here we measure allele-specific transcription factor binding occupancy of three liver-specific transcription factors  between crosses of two inbred mouse strains to elucidate the regulatory mechanisms underlying transcription factor  binding variations in mammals. Our results highlight the pre-eminence of cis-acting variants on transcription factor occupancy divergence. Transcription factor binding differences linked to cis-acting variants generally exhibit additive inheritance, while those linked to trans-acting variants are most often dominantly inherited. Cis-acting variants lead to local coordination of transcription factor occupancies that decay with distance; distal coordination is also observed and may be modulated by long-range chromatin contacts. Our results reveal the regulatory mechanisms that interplay to drive transcription factor occupancy, chromatin state, and gene expression in complex mammalian cell states.

Full details have been published in  Nature Communications  .

Raw Data

The raw ChIP-seq data for can be found in ArrayExpress with the accession number  E-MTAB-4089.

Processed Data

Processed data is available for FOXA1, CEBFA, HNF4A and H3K4me3 as a binary format Excel workbook (xlsb) linked  here .  Column headings for the data spreadsheets are:

chr_position Location of SNV underlying TFBS
suffix of .b Normalized counts for BL6 F0 individuals
suffix of .c Normalized counts for CAST F0 individuals
r Dispersion estimate
suffix of .bi Counts for BL6 allele in BL6xCAST individuals – where total counts (sum of allelic counts) have been adjusted for sequencing depth disparities across F1 libraries 
suffix of .br Counts for BL6 allele in CASTxBL6 individuals –  where total counts (sum of allelic counts) have been adjusted for sequencing depth disparities across F1 libraries 
suffix of .ci Counts for CAST allele in BL6xCAST individuals –  where total counts (sum of allelic counts) have been adjusted for sequencing depth disparities across F1 libraries 
suffix of .cr Counts for CAST allele in CASTxBL6 individuals –  where total counts (sum of allelic counts) have been adjusted for sequencing depth disparities across F1 libraries 
cons likelihood Llikelihood of data fitting the cons model
cis likelihood Likelihood of data fitting the cis model
trans likelihood Likelihood of data fitting the trans model
cistrans likelihood Llikelihood of data fitting the cistrans model
cons_bic BIC for cons  
cis_bic BIC for cis  
trans_bic BIC for trans
cistrans_bic BIC for cistrans
cat Category with lowest BIC

Mitochondrial heteroplasmy in vertebrates using ChIP-sequencing data

Posted on by
  • T Rensch, D Villar, J Horvath, DT Odom, P Flicek. Mitochondrial heteroplasmy in vertebrates using ChIP-sequencing data. Genome Biol 2016;17(1):139. doi:10.1186/s13059-016-0996-y
    [BibTeX] [Abstract]

    \textbf{BACKGROUND:} Mitochondrial heteroplasmy, the presence of more than one mitochondrial DNA (mtDNA) variant in a cell or individual, is not as uncommon as previously thought. It is mostly due to the high mutation rate of the mtDNA and limited repair mechanisms present in the mitochondrion. Motivated by mitochondrial diseases, much focus has been placed into studying this phenomenon in human samples and in medical contexts. To place these results in an evolutionary context and to explore general principles of heteroplasmy, we describe an integrated cross-species evaluation of heteroplasmy in mammals that exploits previously reported NGS data. Focusing on ChIP-seq experiments, we developed a novel approach to detect heteroplasmy from the concomitant mitochondrial DNA fraction sequenced in these experiments.
    \textbf{RESULTS:} We first demonstrate that the sequencing coverage of mtDNA in ChIP-seq experiments is sufficient for heteroplasmy detection. We then describe a novel detection method for accurate detection of heteroplasmies, which also accounts for the error rate of NGS technology. Applying this method to 79 individuals from 16 species resulted in 107 heteroplasmic positions present in a total of 45 individuals. Further analysis revealed that the majority of detected heteroplasmies occur in intergenic regions.
    \textbf{CONCLUSION:} In addition to documenting the prevalence of mtDNA in ChIP-seq data, the results of our mitochondrial heteroplasmy detection method suggest that mitochondrial heteroplasmies identified across vertebrates share similar characteristics as found for human heteroplasmies. Although largely consistent with previous studies in individual vertebrates, our integrated cross-species analysis provides valuable insights into the evolutionary dynamics of mitochondrial heteroplasmy

    @Article{27349964,
author = {Rensch T and Villar D and Horvath J and Odom DT and Flicek P},
title = {Mitochondrial heteroplasmy in vertebrates using ChIP-sequencing data},
journal = {Genome Biol},
volume = {17},
number = {1},
pages = {139},
year = {2016},
doi = {10.1186/s13059-016-0996-y},
abstract = {\textbf{BACKGROUND:} Mitochondrial heteroplasmy, the presence of more than one mitochondrial DNA (mtDNA) variant in a cell or individual, is not as uncommon as previously thought. It is mostly due to the high mutation rate of the mtDNA and limited repair mechanisms present in the mitochondrion. Motivated by mitochondrial diseases, much focus has been placed into studying this phenomenon in human samples and in medical contexts. To place these results in an evolutionary context and to explore general principles of heteroplasmy, we describe an integrated cross-species evaluation of heteroplasmy in mammals that exploits previously reported NGS data. Focusing on ChIP-seq experiments, we developed a novel approach to detect heteroplasmy from the concomitant mitochondrial DNA fraction sequenced in these experiments.
    \textbf{RESULTS:} We first demonstrate that the sequencing coverage of mtDNA in ChIP-seq experiments is sufficient for heteroplasmy detection. We then describe a novel detection method for accurate detection of heteroplasmies, which also accounts for the error rate of NGS technology. Applying this method to 79 individuals from 16 species resulted in 107 heteroplasmic positions present in a total of 45 individuals. Further analysis revealed that the majority of detected heteroplasmies occur in intergenic regions.
    \textbf{CONCLUSION:} In addition to documenting the prevalence of mtDNA in ChIP-seq data, the results of our mitochondrial heteroplasmy detection method suggest that mitochondrial heteroplasmies identified across vertebrates share similar characteristics as found for human heteroplasmies. Although largely consistent with previous studies in individual vertebrates, our integrated cross-species analysis provides valuable insights into the evolutionary dynamics of mitochondrial heteroplasmy},}

Description

Mitochondrial heteroplasmy, the presence of more than one mtDNA variant in a cell or individual is not as uncommon as previously thought. It is mostly due to the high mutation rate of the mtDNA and limited repair mechanisms present in the mitochondrion. Motivated by mitochondrial diseases, much focus has been placed into studying this phenomenon in human samples and in medical contexts. To place these results in an evolutionary context and to explore general principles of heteroplasmy, we describe an integrated cross-species evaluation of heteroplasmy in mammals that exploits previously reported NGS data. Focusing on ChIP-seq experiments, we developed a novel approach to detect mitochondrial heteroplasmy from the concomitant mitochondrial DNA fraction sequenced in these experiments.
We first demonstrate that the sequencing coverage of mtDNA in ChIP-sequencing experiments is sufficient for heteroplasmy detection. We then describe a novel detection method for accurate detection of heteroplasmies, which also accounts for the error rate of NGS technology. Applying this method to 79 individuals from 16 species resulted in 107 heteroplasmic positions present in a total of 45 individuals. Further analysis revealed that the majority of detected heteroplasmies occur in intergenic regions.
In addition to documenting the prevalence of mtDNA in ChIP-sequencing data, the results of our mitochondrial heteroplasmy detection method suggest that mitochondrial heteroplasmies identified across vertebrates share similar characteristics as found for human heteroplasmies. Although largely consistent with previous studies in individual vertebrates, our integrated cross-species analysis provides valuable insights into the evolutionary dynamics of mitochondrial heteroplasmy.

Raw Data

Newly created raw ChIP-seq data for CEBPA, H3K4me1, H3K27ac, and Histone3 can be found in ArrayExpress with the accession number  E-MTAB-3933.
We also used the following previously published ChIP-seq datasets for this study:
HNF4A and CEBPA data from  Schmidt, Wilson, Ballester, et al.  with accession number  E-TABM-722 .
CTCF, SA1, NRSF/REST and H2AK5ac data from  Schmidt, Schwalie, et al.  with accession number  E-MTAB-437 .
CTCF and YY1 ChIP-seq data  Schwalie, Ward, et al.  with accession number  E-MTAB-1511 .
CEBPA, FOXA1, ONECUT1, and HNF4A ChIP-seq data from Ballester, et al. with accession number  E-MTAB-1509 .
H3K4me3 and H3K27ac ChIP-seq from  Villar, Berthelot, et al.  with accession number  E-MTAB-2633 .

Enhancer evolution across twenty mammalian species

Posted on by
  • D Villar, C Berthelot, S Aldridge, TF Rayner, M Lukk, M Pignatelli, TJ Park, R Deaville, JT Erichsen, AJ Jasinska, JMA Turner, MF Bertelsen, EP Murchison, P Flicek, DT Odom. Enhancer Evolution across 20 Mammalian Species. Cell 2015;160(3):554–566. doi:10.1016/j.cell.2015.01.006
    [BibTeX] [Abstract]

    The mammalian radiation has corresponded with rapid changes in noncoding regions of the genome, but we lack a comprehensive understanding of regulatory evolution in mammals. Here, we track the evolution of promoters and enhancers active in liver across 20 mammalian species from six diverse orders by profiling genomic enrichment of H3K27 acetylation and H3K4 trimethylation. We report that rapid evolution of enhancers is a universal feature of mammalian genomes. Most of the recently evolved enhancers arise from ancestral DNA exaptation, rather than lineage-specific expansions of repeat elements. In contrast, almost all liver promoters are partially or fully conserved across these species. Our data further reveal that recently evolved enhancers can be associated with genes under positive selection, demonstrating the power of this approach for annotating regulatory adaptations in genomic sequences. These results provide important insight into the functional genetics underpinning mammalian regulatory evolution

    @Article{25635462,
author = {Villar D and Berthelot C and Aldridge S and Rayner TF and Lukk M and Pignatelli M and Park TJ and Deaville R and Erichsen JT and Jasinska AJ and Turner JMA and Bertelsen MF and Murchison EP and Flicek P and Odom DT},
title = {Enhancer Evolution across 20 Mammalian Species},
journal = {Cell},
volume = {160},
number = {3},
pages = {554--566},
year = {2015},
doi = {10.1016/j.cell.2015.01.006},
abstract = {The mammalian radiation has corresponded with rapid changes in noncoding regions of the genome, but we lack a comprehensive understanding of regulatory evolution in mammals. Here, we track the evolution of promoters and enhancers active in liver across 20 mammalian species from six diverse orders by profiling genomic enrichment of H3K27 acetylation and H3K4 trimethylation. We report that rapid evolution of enhancers is a universal feature of mammalian genomes. Most of the recently evolved enhancers arise from ancestral DNA exaptation, rather than lineage-specific expansions of repeat elements. In contrast, almost all liver promoters are partially or fully conserved across these species. Our data further reveal that recently evolved enhancers can be associated with genes under positive selection, demonstrating the power of this approach for annotating regulatory adaptations in genomic sequences. These results provide important insight into the functional genetics underpinning mammalian regulatory evolution},}

Description

The mammalian radiation has corresponded with rapid changes in the noncoding genome, but we lack a comprehensive understanding of regulatory evolution in mammals. Here, we track the evolution of promoters and enhancers active in liver across twenty mammals from six diverse orders by profiling genomic enrichment of H3K27 acetylation and H3K4 trimethylation. We report rapid evolution of enhancers as a universal feature of mammalian genomes: across all study species half of all liver enhancers are unique to a single species. Most of these recently-evolved enhancers arise from exaptation of ancestral DNA, and not from lineage-specific expansions of repeat elements. In contrast, almost all liver promoters are partially or fully conserved across our study species. Recently-evolved enhancers can be significantly associated with genes under positive selection, demonstrating a powerful approach to annotating regulatory adaptations in newly sequenced genomes. These results provide unprecedented insight into the functional genetics underpinning mammalian regulatory evolution.

Raw Data

The raw chip-seq data for H3K4me3 and H3K27ac from the 20 mammalian species species can be found in ArrayExpress with the accession number  E-MTAB-2633.

ChIP-seq peaks

Consensus MACS peaks reproducible across two or more biological replicates for H3K4me3 and H3K27ac can be found in ArrayExpress with the accession number  E-MTAB-2633.
The format used for the consensus peak files is a custom bed file, with columns as detailed below (tab-separated):
[Chromosome] [Start] [End] [Unique Peak Identifier]

Combined peak calls for H3K4me3, H3K4me3&H3K27ac or H3K27ac in each species can be found at  Combined Peak Calls  with columns as follows (tab-separated):
[Chromosome] [Start] [End] [Identifiers of combined peaks, separated by slashes]

Highly conserved promoters and enhancers

The folder below contains BED files of highly conserved promoters and enhancers identified in the study:

Highly Conserved Elements

Columns are as follows (tab-separated; see header for exact column-matching to species):
[Human chromosome] [Human Start] [Human End] [Human Identifier] [Identifier in species 1] … [Identifier in species 20]

When the reference region in human could be  mapped  in another species but the activity of the region was  not conserved  in this species, the species’ column contains a dash (“-“).
When the reference region in human could  not be mapped  in another species (no alignment), this species’ column is marked as not available (“NA”).

Classification into promoters and enhancers was based on the histone marking observed across the majority of species, and may not always be in agreement with the marking found in the reference human region.

Lineage-specific promoters and enhancers

The folder below contains BED files of lineage-specific promoters and enhancers identified in primates, rodents, ungulates and carnivores, with column headers as above:

Lineage-specific Elements

Different lineages use different reference species, and the reference species is listed in the file name. Column counts vary depending on available alignments between the reference species and the other study species.

Recently evolved promoters and enhancers

The folder below contains BED files of recently-evolved promoters and enhancers identified in human, mouse, cow, dog, naked mole rat and dolphin, with column headers as above:

Recently Evolved Elements

Conservation of human, mouse, cow and dog promoters and enhancers

The folder below contains BED files of all promoters and enhancers identified in human, mouse, cow and dog, and their conservation in other species, with column headers as above:

Conservation

Decoupling of evolutionary changes in transcription factor binding and gene expression in mammals

Posted on by
  • ES Wong, D Thybert, BM Schmitt, K Stefflova, DT Odom, P Flicek. Decoupling of evolutionary changes in transcription factor binding and gene expression in mammals. Genome Res 2015;25(2):167–178. doi:10.1101/gr.177840.114
    [BibTeX] [Abstract]

    To understand the evolutionary dynamics between transcription factor (TF) binding and gene expression in mammals, we compared transcriptional output and the binding intensities for three tissue-specific TFs in livers from four closely related mouse species. For each transcription factor, TF dependent genes and the TF binding sites most likely to influence mRNA expression were identified by comparing mRNA expression levels between wildtype and TF knockout mice. Independent evolution was observed genome-wide between the rate of change in TF binding and the rate of change in mRNA expression across taxa, with the exception of a small number of TF dependent genes. We also found that binding intensities are preferentially conserved near genes whose expression is dependent on the TF, and the conservation is shared among binding peaks in close proximity to each other near the TSS. Expression of TF dependent genes typically showed an increased sensitivity to changes in binding levels, as measured by mRNA abundance. Taken together, these results highlight a significant tolerance to evolutionary changes in TF binding intensity in mammalian transcriptional networks, and suggest that some TF dependent genes may be largely regulated by a single TF across evolution

    @Article{25394363,
author = {Wong ES and Thybert D and Schmitt BM and Stefflova K and Odom DT and Flicek P},
title = {Decoupling of evolutionary changes in transcription factor binding and gene expression in mammals},
journal = {Genome Res},
volume = {25},
number = {2},
pages = {167--178},
year = {2015},
doi = {10.1101/gr.177840.114},
howpublished = {Advanced online publication: 13 November 2014},
abstract = {To understand the evolutionary dynamics between transcription factor (TF) binding and gene expression in mammals, we compared transcriptional output and the binding intensities for three tissue-specific TFs in livers from four closely related mouse species. For each transcription factor, TF dependent genes and the TF binding sites most likely to influence mRNA expression were identified by comparing mRNA expression levels between wildtype and TF knockout mice. Independent evolution was observed genome-wide between the rate of change in TF binding and the rate of change in mRNA expression across taxa, with the exception of a small number of TF dependent genes. We also found that binding intensities are preferentially conserved near genes whose expression is dependent on the TF, and the conservation is shared among binding peaks in close proximity to each other near the TSS. Expression of TF dependent genes typically showed an increased sensitivity to changes in binding levels, as measured by mRNA abundance. Taken together, these results highlight a significant tolerance to evolutionary changes in TF binding intensity in mammalian transcriptional networks, and suggest that some TF dependent genes may be largely regulated by a single TF across evolution},}

Raw Data

The raw gene expression data for mouse species can be found in ArrayExpress  with the accesion number  E-MTAB-2483.
The raw gene expression data for HNF4A knockout can be found in ArrayExpress with the accession number  E-MTAB-2484.

Processed data

Processed gene expression data and list of target genes can be found  here .

R code

Snippets of R code for analyses detailed in the manuscript can be found  here .

Complete peak calls

Datasets are from  “Cooperativity and rapid evolution of cobound transcription factors in closely related mammals” , Stefflova and Thybert et al. Cell 2013, and processed as described in their manuscript.
See the description below for the meaning of each data field:

chr:  chromosome
start_bl6 :  start of binding region wiht C57BL/6J coordinate system
end_bl6 :  start of binding region wiht C57BL/6J coordinate system
summit_bl6 :  summit of binding region with C57BL/6J coordinate system
intensity :  intensity in normalized read count
intensity_class :  intensity class
start_own :  start in species coordinate system
end_own :  end in species coordinate system
summit_own :  summit position in specie coordinate system

CEBPA peak calls for 5 mouse species
HNF4A peak calls for 5 mouse species
FOXA1 peak calls for 5 mouse species