On the estimation of genome-average recombination rates

The rate at which recombination events occur in a population is an indicator of its effective population size and the organism's reproduction mode. It determines the extent of linkage disequilibrium along the genome and, thereby, the efficacy of both purifying and positive selection. The population recombination rate can be inferred using models of genome evolution in populations. Classic methods based on the patterns of linkage-disequilibrium provide the most accurate estimates, providing large sample sizes are used and the demography of the population is properly accounted for. Here, the capacity of approaches based on the sequentially Markov coalescent (SMC) to infer the genome-average recombination rate from as little as a single diploid genome is examined. SMC approaches provide highly accurate estimates even in the presence of changing population sizes, providing that (1) within genome heterogeneity is accounted for and (2) classic maximum-likelihood optimization algorithms are employed to fit the model. SMC-based estimates proved sensitive to gene conversion, leading to an overestimation of the recombination rate if conversion events are frequent. Conversely, methods based on the correlation of heterozygosity succeed in disentangling the rate of crossing over from that of gene conversion events, but only when the population size is constant and the recombination landscape homogeneous. These results call for a convergence of these two methods to obtain accurate and comparable estimates of recombination rates between populations.

Adaptive evolution and loss of a putative magnetoreceptor in passerines

Migratory birds possess remarkable accuracy in orientation and navigation, which involves various compass systems including the magnetic compass. Identifying the primary magnetosensor remains a fundamental open question. Cryptochromes (Cry) have been shown to be magnetically sensitive, and Cry4a from a migratory songbird seems to show enhanced magnetic sensitivity in vitro compared to Cry4a from resident species. We investigate Cry and their potential involvement in magnetoreception in a phylogenetic framework, integrating molecular evolutionary analyses with protein dynamics modelling. Our analysis is based on 363 bird genomes and identifies different selection regimes in passerines. We show that Cry4a is characterized by strong positive selection and high variability, typical characteristics of sensor proteins. We identify key sites that are likely to have facilitated the evolution of an optimized sensory protein for night-time orientation in songbirds. Additionally, we show that Cry4 was lost in hummingbirds, parrots and Tyranni (Suboscines), and thus identified a gene deletion, which might facilitate testing the function of Cry4a in birds. In contrast, the other avian Cry (Cry1 and Cry2) were highly conserved across all species, indicating basal, non-sensory functions. Our results support a specialization or functional differentiation of Cry4 in songbirds which could be magnetosensation.

The Site/Group Extended Data Format and Tools

Comparative sequence analysis permits unraveling the molecular processes underlying gene evolution. Many statistical methods generate candidate positions within genes, such as fast or slowly evolving sites, coevolving groups of residues, sites undergoing positive selection, or changes in evolutionary rates. Understanding the functional causes of these evolutionary patterns requires combining the results of these analyses and mapping them onto molecular structures, a complex task involving distinct coordinate referential systems. To ease this task, we introduce the site/group extended data format, a simple text format to store (groups of) site annotations. We developed a toolset, the SgedTools, which permits site/group extended data file manipulation, creating them from various software outputs and translating coordinates between individual sequences, alignments, and three-dimensional structures. The package also includes a Monte-Carlo procedure to generate random site samples, possibly conditioning on site-specific features. This eases the statistical testing of evolutionary hypotheses, accounting for the structural properties of the encoded molecules.

Fine-Scale Map Reveals Highly Variable Recombination Rates Associated with Genomic Features in the Eurasian Blackcap

Recombination is responsible for breaking up haplotypes, influencing genetic variability, and the efficacy of selection. Bird genomes lack the protein PR domain-containing protein 9, a key determinant of recombination dynamics in most metazoans. Historical recombination maps in birds show an apparent stasis in positioning recombination events. This highly conserved recombination pattern over long timescales may constrain the evolution of recombination in birds. At the same time, extensive variation in recombination rate is observed across the genome and between different species of birds. Here, we characterize the fine-scale historical recombination map of an iconic migratory songbird, the Eurasian blackcap (Sylvia atricapilla), using a linkage disequilibrium-based approach that accounts for population demography. Our results reveal variable recombination rates among and within chromosomes, which associate positively with nucleotide diversity and GC content and negatively with chromosome size. Recombination rates increased significantly at regulatory regions but not necessarily at gene bodies. CpG islands are associated strongly with recombination rates, though their specific position and local DNA methylation patterns likely influence this relationship. The association with retrotransposons varied according to specific family and location. Our results also provide evidence of heterogeneous intrachromosomal conservation of recombination maps between the blackcap and its closest sister taxon, the garden warbler. These findings highlight the considerable variability of recombination rates at different scales and the role of specific genomic features in shaping this variation. This study opens the possibility of further investigating the impact of recombination on specific population-genomic features.

Pervasive incomplete lineage sorting illuminates speciation and selection in primates

Incomplete lineage sorting (ILS) causes the phylogeny of some parts of the genome to differ from the species tree. In this work, we investigate the frequencies and determinants of ILS in 29 major ancestral nodes across the entire primate phylogeny. We find up to 64% of the genome affected by ILS at individual nodes. We exploit ILS to reconstruct speciation times and ancestral population sizes. Estimated speciation times are much more recent than genomic divergence times and are in good agreement with the fossil record. We show extensive variation of ILS along the genome, mainly driven by recombination but also by the distance to genes, highlighting a major impact of selection on variation along the genome. In many nodes, ILS is reduced more on the X chromosome compared with autosomes than expected under neutrality, which suggests higher impacts of natural selection on the X chromosome. Finally, we show an excess of ILS in genes with immune functions and a deficit of ILS in housekeeping genes. The extensive ILS in primates discovered in this study provides insights into the speciation times, ancestral population sizes, and patterns of natural selection that shape primate evolution.

The rate of adaptive molecular evolution in wild and domesticated Saccharomyces cerevisiae populations

Through its fermentative capacities, Saccharomyces cerevisiae was central in the development of civilisation during the Neolithic period, and the yeast remains of importance in industry and biotechnology, giving rise to bona fide domesticated populations. Here, we conduct a population genomic study of domesticated and wild populations of S. cerevisiae. Using coalescent analyses, we report that the effective population size of yeast populations decreased since the divergence with S. paradoxus. We fitted models of distributions of fitness effects to infer the rate of adaptive ( ω a $$ {\omega}_a $$ ) and non-adaptive ( ω na $$ {\omega}_{na} $$ ) non-synonymous substitutions in protein-coding genes. We report an overall limited contribution of positive selection to S. cerevisiae protein evolution, albeit with higher rates of adaptive evolution in wild compared to domesticated populations. Our analyses revealed the signature of background selection and possibly Hill-Robertson interference, as recombination was found to be negatively correlated with ω na $$ {\omega}_{na} $$ and positively correlated with ω a $$ {\omega}_a $$ . However, the effect of recombination on ω a $$ {\omega}_a $$ was found to be labile, as it is only apparent after removing the impact of codon usage bias on the synonymous site frequency spectrum and disappears if we control for the correlation with ω na $$ {\omega}_{na} $$ , suggesting that it could be an artefact of the decreasing population size. Furthermore, the rate of adaptive non-synonymous substitutions is significantly correlated with the residue solvent exposure, a relation that cannot be explained by the population's demography. Together, our results provide a detailed characterisation of adaptive mutations in protein-coding genes across S. cerevisiae populations.

Being noisy in a crowd: Differential selective pressure on gene expression noise in model gene regulatory networks

Expression noise, the variability of the amount of gene product among isogenic cells grown in identical conditions, originates from the inherent stochasticity of diffusion and binding of the molecular players involved in transcription and translation. It has been shown that expression noise is an evolvable trait and that central genes exhibit less noise than peripheral genes in gene networks. A possible explanation for this pattern is increased selective pressure on central genes since they propagate their noise to downstream targets, leading to noise amplification. To test this hypothesis, we developed a new gene regulatory network model with inheritable stochastic gene expression and simulated the evolution of gene-specific expression noise under constraint at the network level. Stabilizing selection was imposed on the expression level of all genes in the network and rounds of mutation, selection, replication and recombination were performed. We observed that local network features affect both the probability to respond to selection, and the strength of the selective pressure acting on individual genes. In particular, the reduction of gene-specific expression noise as a response to stabilizing selection on the gene expression level is higher in genes with higher centrality metrics. Furthermore, global topological structures such as network diameter, centralization and average degree affect the average expression variance and average selective pressure acting on constituent genes. Our results demonstrate that selection at the network level leads to differential selective pressure at the gene level, and local and global network characteristics are an essential component of gene-specific expression noise evolution.

Strong evidence for the adaptive walk model of gene evolution in Drosophila and Arabidopsis

Understanding the dynamics of species adaptation to their environments has long been a central focus of the study of evolution. Theories of adaptation propose that populations evolve by “walking” in a fitness landscape. This “adaptive walk” is characterised by a pattern of diminishing returns, where populations further away from their fitness optimum take larger steps than those closer to their optimal conditions. Hence, we expect young genes to evolve faster and experience mutations with stronger fitness effects than older genes because they are further away from their fitness optimum. Testing this hypothesis, however, constitutes an arduous task. Young genes are small, encode proteins with a higher degree of intrinsic disorder, are expressed at lower levels, and are involved in species-specific adaptations. Since all these factors lead to increased protein evolutionary rates, they could be masking the effect of gene age. While controlling for these factors, we used population genomic data sets of Arabidopsis and Drosophila and estimated the rate of adaptive substitutions across genes from different phylostrata. We found that a gene’s evolutionary age significantly impacts the molecular rate of adaptation. Moreover, we observed that substitutions in young genes tend to have larger physicochemical effects. Our study, therefore, provides strong evidence that molecular evolution follows an adaptive walk model across a large evolutionary timescale.

This study uses population genomic datasets from Arabidopsis and Drosophila to show that young genes adapt faster and are subject to mutations of larger fitness effects, providing strong evidence that molecular evolution follows an adaptive walk model across a large evolutionary timescale.

Plant pathogens provide clues to the potential origin of bat white-nose syndrome Pseudogymnoascus destructans

White-nose syndrome has killed millions of bats, yet both the origins and infection strategy of the causative fungus, Pseudogymnoascus destructans, remain elusive. We provide evidence for a novel hypothesis that P. destructans emerged from plant-associated fungi and retained invasion strategies affiliated with fungal pathogens of plants. We demonstrate that P. destructans invades bat skin in successive biotrophic and necrotrophic stages (hemibiotrophic infection), a mechanism previously only described in plant fungal pathogens. Further, the convergence of hyphae at hair follicles suggests nutrient tropism. Tropism, biotrophy, and necrotrophy are often associated with structures termed appressoria in plant fungal pathogens; the penetrating hyphae produced by P. destructans resemble appressoria. Finally, we conducted a phylogenomic analysis of a taxonomically diverse collection of fungi. Despite gaps in genetic sampling of prehistoric and contemporary fungal species, we estimate an 88% probability the ancestral state of the clade containing P. destructans was a plant-associated fungus.

The Structural Determinants of Intra-Protein Compensatory Substitutions

Compensatory substitutions happen when one mutation is advantageously selected because it restores the loss of fitness induced by a previous deleterious mutation. How frequent such mutations occur in evolution and what is the structural and functional context permitting their emergence remain open questions. We built an atlas of intra-protein compensatory substitutions using a phylogenetic approach and a dataset of 1,630 bacterial protein families for which high-quality sequence alignments and experimentally derived protein structures were available. We identified more than 51,000 positions coevolving by the mean of predicted compensatory mutations. Using the evolutionary and structural properties of the analyzed positions, we demonstrate that compensatory mutations are scarce (typically only a few in the protein history) but widespread (the majority of proteins experienced at least one). Typical coevolving residues are evolving slowly, are located in the protein core outside secondary structure motifs, and are more often in contact than expected by chance, even after accounting for their evolutionary rate and solvent exposure. An exception to this general scheme is residues coevolving for charge compensation, which are evolving faster than noncoevolving sites, in contradiction with predictions from simple coevolutionary models, but similar to stem pairs in RNA. While sites with a significant pattern of coevolution by compensatory mutations are rare, the comparative analysis of hundreds of structures ultimately permits a better understanding of the link between the three-dimensional structure of a protein and its fitness landscape.

Population Genomics of the Maize Pathogen Ustilago maydis: Demographic History and Role of Virulence Clusters in Adaptation

The tight interaction between pathogens and their hosts results in reciprocal selective forces that impact the genetic diversity of the interacting species. The footprints of this selection differ between pathosystems because of distinct life-history traits, demographic histories, or genome architectures. Here, we studied the genome-wide patterns of genetic diversity of 22 isolates of the causative agent of the corn smut disease, Ustilago maydis, originating from five locations in Mexico, the presumed center of origin of this species. In this species, many genes encoding secreted effector proteins reside in so-called virulence clusters in the genome, an arrangement that is so far not found in other filamentous plant pathogens. Using a combination of population genomic statistical analyses, we assessed the geographical, historical, and genome-wide variation of genetic diversity in this fungal pathogen. We report evidence of two partially admixed subpopulations that are only loosely associated with geographic origin. Using the multiple sequentially Markov coalescent model, we inferred the demographic history of the two pathogen subpopulations over the last 0.5 Myr. We show that both populations experienced a recent strong bottleneck starting around 10,000 years ago, coinciding with the assumed time of maize domestication. Although the genome average genetic diversity is low compared with other fungal pathogens, we estimated that the rate of nonsynonymous adaptive substitutions is three times higher in genes located within virulence clusters compared with nonclustered genes, including nonclustered effector genes. These results highlight the role that these singular genomic regions play in the evolution of this pathogen.

The insertion of a mitochondrial selfish element into the nuclear genome and its consequences

Homing endonucleases (HE) are enzymes capable of cutting DNA at highly specific target sequences, the repair of the generated double-strand break resulting in the insertion of the HE-encoding gene ("homing" mechanism). HEs are present in all three domains of life and viruses; in eukaryotes, they are mostly found in the genomes of mitochondria and chloroplasts, as well as nuclear ribosomal RNAs. We here report the case of a HE that accidentally integrated into a telomeric region of the nuclear genome of the fungal maize pathogen Ustilago maydis. We show that the gene has a mitochondrial origin, but its original copy is absent from the U. maydis mitochondrial genome, suggesting a subsequent loss or a horizontal transfer from a different species. The telomeric HE underwent mutations in its active site and lost its original start codon. A potential other start codon was retained downstream, but we did not detect any significant transcription of the newly created open reading frame, suggesting that the inserted gene is not functional. Besides, the insertion site is located in a putative RecQ helicase gene, truncating the C-terminal domain of the protein. The truncated helicase is expressed during infection of the host, together with other homologous telomeric helicases. This unusual mutational event altered two genes: The integrated HE gene subsequently lost its homing activity, while its insertion created a truncated version of an existing gene, possibly altering its function. As the insertion is absent in other field isolates, suggesting that it is recent, the U. maydis 521 reference strain offers a snapshot of this singular mutational event.

On Variant Discovery in Genomes of Fungal Plant Pathogens

Comparative genome analyses of eukaryotic pathogens including fungi and oomycetes have revealed extensive variability in genome composition and structure. The genomes of individuals from the same population can exhibit different numbers of chromosomes and different organization of chromosomal segments, defining so-called accessory compartments that have been shown to be crucial to pathogenicity in plant-infecting fungi. This high level of structural variation confers a methodological challenge for population genomic analyses. Variant discovery from population sequencing data is typically achieved using established pipelines based on the mapping of short reads to a reference genome. These pipelines have been developed, and extensively used, for eukaryote genomes of both plants and animals, to retrieve single nucleotide polymorphisms and short insertions and deletions. However, they do not permit the inference of large-scale genomic structural variation, as this task typically requires the alignment of complete genome sequences. Here, we compare traditional variant discovery approaches to a pipeline based on de novo genome assembly of short read data followed by whole genome alignment, using simulated data sets with properties mimicking that of fungal pathogen genomes. We show that the latter approach exhibits levels of performance comparable to that of read-mapping based methodologies, when used on sequence data with sufficient coverage. We argue that this approach further allows additional types of genomic diversity to be explored, in particular as long-read third-generation sequencing technologies are becoming increasingly available to generate population genomic data.

A Population Genomics Lexicon

Population genomics is a growing field stemming from soon a 100 years of developments in population genetics. Here, we summarize the main concepts and terminology underlying both theoretical and empirical statistical population genomics studies. We provide the reader with pointers toward the original literature as well as methodological and historical reviews.

Correction to: Processing and Analyzing Multiple Genomes Alignments with MafFilter

Chapter 2, “Processing and Analyzing Multiple Genomes Alignments with MafFilter,” was previously published without including the Electronic Supplementary Material. This has now been included in the revised version of this book.

Variation of the adaptive substitution rate between species and within genomes

The importance of adaptive mutations in molecular evolution is extensively debated. Recent developments in population genomics allow inferring rates of adaptive mutations by fitting a distribution of fitness effects to the observed patterns of polymorphism and divergence at sites under selection and sites assumed to evolve neutrally. Here, we summarize the current state-of-the-art of these methods and review the factors that affect the molecular rate of adaptation. Several studies have reported extensive cross-species variation in the proportion of adaptive amino-acid substitutions (α) and predicted that species with larger effective population sizes undergo less genetic drift and higher rates of adaptation. Disentangling the rates of positive and negative selection, however, revealed that mutations with deleterious effects are the main driver of this population size effect and that adaptive substitution rates vary comparatively little across species. Conversely, rates of adaptive substitution have been documented to vary substantially within genomes. On a genome-wide scale, gene density, recombination and mutation rate were observed to play a role in shaping molecular rates of adaptation, as predicted under models of linked selection. At the gene level, it has been reported that the gene functional category and the macromolecular structure substantially impact the rate of adaptive mutations. Here, we deliver a comprehensive review of methods used to infer the molecular adaptive rate, the potential drivers of adaptive evolution and how positive selection shapes molecular evolution within genes, across genes within species and between species.

Inference of recombination maps from a single pair of genomes and its application to ancient samples

Understanding the causes and consequences of recombination landscape evolution is a fundamental goal in genetics that requires recombination maps from across the tree of life. Such maps can be obtained from population genomic datasets, but require large sample sizes. Alternative methods are therefore necessary to research organisms where such datasets cannot be generated easily, such as non-model or ancient species. Here we extend the sequentially Markovian coalescent model to jointly infer demography and the spatial variation in recombination rate. Using extensive simulations and sequence data from humans, fruit-flies and a fungal pathogen, we demonstrate that iSMC accurately infers recombination maps under a wide range of scenarios-remarkably, even from a single pair of unphased genomes. We exploit this possibility and reconstruct the recombination maps of ancient hominins. We report that the ancient and modern maps are correlated in a manner that reflects the established phylogeny of Neanderthals, Denisovans, and modern human populations.

The genomic determinants of adaptive evolution in a fungal pathogen

Unravelling the strength, frequency, and distribution of selective variants along the genome as well as the underlying factors shaping this distribution are fundamental goals of evolutionary biology. Antagonistic host-pathogen coevolution is thought to be a major driver of genome evolution between interacting species. While rapid evolution of pathogens has been documented in several model organisms, the genetic mechanisms of their adaptation are still poorly understood and debated, particularly the role of sexual reproduction. Here, we apply a population genomic approach to infer genome-wide patterns of selection among 13 isolates of Zymoseptoria tritici, a fungal pathogen characterized by extremely high genetic diversity, gene density, and recombination rates. We report that the genome of Z. tritici undergoes a high rate of adaptive substitutions, with 44% of nonsynonymous substitutions being adaptive on average. This fraction reaches 68% in so-called effector genes encoding determinants of pathogenicity, and the distribution of fitness effects differs in this class of genes as they undergo adaptive mutations with stronger positive fitness effects, but also more slightly deleterious mutations. Besides the globally high rate of adaptive substitutions, we report a negative relationship between pN/pS and the fine-scale recombination rate and a strong positive correlation between the rate of adaptive nonsynonymous substitutions (ωa) and recombination rate. This result suggests a pervasive role of both background selection and Hill-Robertson interference even in a species with an exceptionally high recombination rate (60 cM/Mb on average). While transposable elements (TEs) have been suggested to contribute to adaptation by creating compartments of fast-evolving genomic regions, we do not find a significant effect of TEs on the rate of adaptive mutations. Overall our study suggests that sexual recombination is a significant driver of genome evolution, even in rapidly evolving organisms subject to recurrent mutations with large positive effects.

Ancestral Population Genomics

Borrowing both from population genetics and phylogenetics, the field of population genomics emerged as full genomes of several closely related species were available. Providing we can properly model sequence evolution within populations undergoing speciation events, this resource enables us to estimate key population genetics parameters such as ancestral population sizes and split times. Furthermore we can enhance our understanding of the recombination process and investigate various selective forces. With the advent of resequencing technologies, genome-wide patterns of diversity in extant populations have now come to complement this picture, offering an increasing power to study more recent genetic history.We discuss the basic models of genomes in populations, including speciation models for closely related species. A major point in our discussion is that only a few complete genomes contain much information about the whole population. The reason being that recombination unlinks genomic regions, and therefore a few genomes contain many segments with distinct histories. The challenge of population genomics is to decode this mosaic of histories in order to infer scenarios of demography and selection. We survey modeling strategies for understanding genetic variation in ancestral populations and species. The underlying models build on the coalescent with recombination process and introduce further assumptions to scale the analyses to genomic data sets.

Fine-Scale Recombination Maps of Fungal Plant Pathogens Reveal Dynamic Recombination Landscapes and Intragenic Hotspots

Meiotic recombination is an important driver of evolution. Variability in the intensity of recombination across chromosomes can affect sequence composition, nucleotide variation, and rates of adaptation. In many organisms, recombination events are concentrated within short segments termed recombination hotspots. The variation in recombination rate and positions of recombination hotspot can be studied using population genomics data and statistical methods. In this study, we conducted population genomics analyses to address the evolution of recombination in two closely related fungal plant pathogens: the prominent wheat pathogen Zymoseptoria tritici and a sister species infecting wild grasses Z. ardabiliae We specifically addressed whether recombination landscapes, including hotspot positions, are conserved in the two recently diverged species and if recombination contributes to rapid evolution of pathogenicity traits. We conducted a detailed simulation analysis to assess the performance of methods of recombination rate estimation based on patterns of linkage disequilibrium, in particular in the context of high nucleotide diversity. Our analyses reveal overall high recombination rates, a lack of suppressed recombination in centromeres, and significantly lower recombination rates on chromosomes that are known to be accessory. The comparison of the recombination landscapes of the two species reveals a strong correlation of recombination rate at the megabase scale, but little correlation at smaller scales. The recombination landscapes in both pathogen species are dominated by frequent recombination hotspots across the genome including coding regions, suggesting a strong impact of recombination on gene evolution. A significant but small fraction of these hotspots colocalize between the two species, suggesting that hotspot dynamics contribute to the overall pattern of fast evolving recombination in these species.

Positively Selected Effector Genes and Their Contribution to Virulence in the Smut Fungus Sporisorium reilianum

Plants and fungi display a broad range of interactions in natural and agricultural ecosystems ranging from symbiosis to parasitism. These ecological interactions result in coevolution between genes belonging to different partners. A well-understood example is secreted fungal effector proteins and their host targets, which play an important role in pathogenic interactions. Biotrophic smut fungi (Basidiomycota) are well-suited to investigate the evolution of plant pathogens, because several reference genomes and genetic tools are available for these species. Here, we used the genomes of Sporisorium reilianum f. sp. zeae and S. reilianum f. sp. reilianum, two closely related formae speciales infecting maize and sorghum, respectively, together with the genomes of Ustilago hordei, Ustilago maydis, and Sporisorium scitamineum to identify and characterize genes displaying signatures of positive selection. We identified 154 gene families having undergone positive selection during species divergence in at least one lineage, among which 77% were identified in the two investigated formae speciales of S. reilianum. Remarkably, only 29% of positively selected genes encode predicted secreted proteins. We assessed the contribution to virulence of nine of these candidate effector genes in S. reilianum f. sp. zeae by deleting individual genes, including a homologue of the effector gene pit2 previously characterized in U. maydis. Only the pit2 deletion mutant was found to be strongly reduced in virulence. Additional experiments are required to understand the molecular mechanisms underlying the selection forces acting on the other candidate effector genes, as well as the large fraction of positively selected genes encoding predicted cytoplasmic proteins.

The Evolution of Gene-Specific Transcriptional Noise Is Driven by Selection at the Pathway Level

Biochemical reactions within individual cells result from the interactions of molecules, typically in small numbers. Consequently, the inherent stochasticity of binding and diffusion processes generates noise along the cascade that leads to the synthesis of a protein from its encoding gene. As a result, isogenic cell populations display phenotypic variability even in homogeneous environments. The extent and consequences of this stochastic gene expression have only recently been assessed on a genome-wide scale, owing, in particular, to the advent of single-cell transcriptomics. However, the evolutionary forces shaping this stochasticity have yet to be unraveled. Here, we take advantage of two recently published data sets for the single-cell transcriptome of the domestic mouse Mus musculus to characterize the effect of natural selection on gene-specific transcriptional stochasticity. We show that noise levels in the mRNA distributions (also known as transcriptional noise) significantly correlate with three-dimensional nuclear domain organization, evolutionary constraints on the encoded protein, and gene age. However, the position of the encoded protein in a biological pathway is the main factor that explains observed levels of transcriptional noise, in agreement with models of noise propagation within gene networks. Because transcriptional noise is under widespread selection, we argue that it constitutes an important component of the phenotype and that variance of expression is a potential target of adaptation. Stochastic gene expression should therefore be considered together with the mean expression level in functional and evolutionary studies of gene expression.

hotspot: software to support sperm-typing for investigating recombination hotspots

MOTIVATION

In many organisms, including humans, recombination clusters within recombination hotspots. The standard method for de novo detection of recombinants at hotspots is sperm typing. This relies on allele-specific PCR at single nucleotide polymorphisms. Designing allele-specific primers by hand is time-consuming. We have therefore written a package to support hotspot detection and analysis.

RESULTS

hotspot consists of four programs: asp looks up SNPs and designs allele-specific primers; aso constructs allele-specific oligos for mapping recombinants; xov implements a maximum-likelihood method for estimating the crossover rate; six, finally, simulates typing data.

AVAILABILITY AND IMPLEMENTATION

hotspot is written in C. Sources are freely available under the GNU General Public License from http://github.com/evolbioinf/hotspot/

CONTACT

haubold@evolbio.mpg.de

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

The WOPR Protein Ros1 Is a Master Regulator of Sporogenesis and Late Effector Gene Expression in the Maize Pathogen Ustilago maydis

The biotrophic basidiomycete fungus Ustilago maydis causes smut disease in maize. Hallmarks of the disease are large tumors that develop on all aerial parts of the host in which dark pigmented teliospores are formed. We have identified a member of the WOPR family of transcription factors, Ros1, as major regulator of spore formation in U. maydis. ros1 expression is induced only late during infection and hence Ros1 is neither involved in plant colonization of dikaryotic fungal hyphae nor in plant tumor formation. However, during late stages of infection Ros1 is essential for fungal karyogamy, massive proliferation of diploid fungal cells and spore formation. Premature expression of ros1 revealed that Ros1 counteracts the b-dependent filamentation program and induces morphological alterations resembling the early steps of sporogenesis. Transcriptional profiling and ChIP-seq analyses uncovered that Ros1 remodels expression of about 30% of all U. maydis genes with 40% of these being direct targets. In total the expression of 80 transcription factor genes is controlled by Ros1. Four of the upregulated transcription factor genes were deleted and two of the mutants were affected in spore development. A large number of b-dependent genes were differentially regulated by Ros1, suggesting substantial changes in this regulatory cascade that controls filamentation and pathogenic development. Interestingly, 128 genes encoding secreted effectors involved in the establishment of biotrophic development were downregulated by Ros1 while a set of 70 “late effectors” was upregulated. These results indicate that Ros1 is a master regulator of late development in U. maydis and show that the biotrophic interaction during sporogenesis involves a drastic shift in expression of the fungal effectome including the downregulation of effectors that are essential during early stages of infection.

The fungus Ustilago maydis is a pathogen of maize which induces tumor formation in the infected tissue. In these tumors huge amounts of fungal spores develop. As a biotrophic pathogen, U. maydis establishes itself in the plant with the help of a large number of secreted effector proteins. Many effector proteins are important for virulence because they counteract plant defense reactions. In this manuscript we have identified and characterized Ros1, a master regulator for the late stages of U. maydis development. This transcription factor is expressed late during infection and controls nuclear fusion, hyphal aggregation and late proliferation. ros1 mutants are still able to induce tumor formation but these are a dead end because they do not contain any spores. We show that Ros1 interferes with the early regulatory cascade controlled by a complex of two homeodomain proteins. In addition, Ros1 triggers a major switch in the effector repertoire, suggesting that different sets of effectors are needed for different stages of fungal development inside the plant.

A Tale of Genome Compartmentalization: The Evolution of Virulence Clusters in Smut Fungi

Smut fungi are plant pathogens mostly parasitizing wild species of grasses as well as domesticated cereal crops. Genome analysis of several smut fungi including Ustilago maydis revealed a singular clustered organization of genes encoding secreted effectors. In U. maydis, many of these clusters have a role in virulence. Reconstructing the evolutionary history of clusters of effector genes is difficult because of their intrinsically fast evolution, which erodes the phylogenetic signal and homology relationships. Here, we describe the use of comparative evolutionary analyses of quality draft assemblies of genomes to study the mechanisms of this evolution. We report the genome sequence of a South African isolate of Sporisorium scitamineum, a smut fungus parasitizing sugar cane with a phylogenetic position intermediate to the two previously sequenced species U. maydis and Sporisorium reilianum. We show that the genome of S. scitamineum contains more and larger gene clusters encoding secreted effectors than any previously described species in this group. We trace back the origin of the clusters and find that their evolution is mainly driven by tandem gene duplication. In addition, transposable elements play a major role in the evolution of the clustered genes. Transposable elements are significantly associated with clusters of genes encoding fast evolving secreted effectors. This suggests that such clusters represent a case of genome compartmentalization that restrains the activity of transposable elements on genes under diversifying selection for which this activity is potentially beneficial, while protecting the rest of the genome from its deleterious effect.

Strong Selective Sweeps on the X Chromosome in the Human-Chimpanzee Ancestor Explain Its Low Divergence

The human and chimpanzee X chromosomes are less divergent than expected based on autosomal divergence. We study incomplete lineage sorting patterns between humans, chimpanzees and gorillas to show that this low divergence can be entirely explained by megabase-sized regions comprising one-third of the X chromosome, where polymorphism in the human-chimpanzee ancestral species was severely reduced. We show that background selection can explain at most 10% of this reduction of diversity in the ancestor. Instead, we show that several strong selective sweeps in the ancestral species can explain it. We also report evidence of population specific sweeps in extant humans that overlap the regions of low diversity in the ancestral species. These regions further correspond to chromosomal sections shown to be devoid of Neanderthal introgression into modern humans. This suggests that the same X-linked regions that undergo selective sweeps are among the first to form reproductive barriers between diverging species. We hypothesize that meiotic drive is the underlying mechanism causing these two observations.

Optimization of sequence alignments according to the number of sequences vs. number of sites trade-off

Background

Comparative analysis of homologous sequences enables the understanding of evolutionary patterns at the molecular level, unraveling the functional constraints that shaped the underlying genes. Bioinformatic pipelines for comparative sequence analysis typically include procedures for (i) alignment quality assessment and (ii) control of sequence redundancy. An additional, underassessed step is the control of the amount and distribution of missing data in sequence alignments. While the number of sequences available for a given gene typically increases with time, the site-specific coverage of each alignment position remains highly variable because of differences in sequencing and annotation quality, or simply because of biological variation. For any given alignment-based analysis, the selection of sequences thus defines a trade-off between the species representation and the quantity of sites with sufficient coverage to be included in the subsequent analyses.

Results

We introduce an algorithm for the optimization of sequence alignments according to the number of sequences vs. number of sites trade-off. The algorithm uses a guide tree to compute scores for each bipartition of the alignment, allowing the recursive selection of sequence subsets with optimal combinations of sequence and site numbers. By applying our methods to two large data sets of several thousands of gene families, we show that significant site-specific coverage increases can be achieved while controlling for the species representation.

Conclusions

The algorithm introduced in this work allows the control of the distribution of missing data in any sequence alignment by removing sequences to increase the number of sites with a defined minimum coverage. We advocate that our missing data optimization procedure in an important step which should be considered in comparative analysis pipelines, together with alignment quality assessment and control of sampled diversity. An open source C++ implementation is available at http://bioweb.me/physamp.

Electronic supplementary material

The online version of this article (doi:10.1186/s12859-015-0619-8) contains supplementary material, which is available to authorized users.

Mitochondrial DNA as a tool for reconstructing past life-history traits in mammals

Reconstructing the ancestral characteristics of species is a major goal in evolutionary and comparative biology. Unfortunately, fossils are not always available and sufficiently informative, and phylogenetic methods based on models of character evolution can be unsatisfactory. Genomic data offer a new opportunity to estimate ancestral character states, through (i) the correlation between DNA evolutionary processes and species life-history traits and (ii) available reliable methods for ancestral sequence inference. Here, we assess the relevance of mitochondrial DNA--the most popular molecular marker in animals--as a predictor of ancestral life-history traits in mammals, using the order of Cetartiodactyla as a benchmark. Using the complete set of 13 mitochondrial protein-coding genes, we show that the lineage-specific nonsynonymous over synonymous substitution rate ratio (dN/dS) is closely correlated with the species body mass, longevity and age of sexual maturity in Cetartiodactyla and can be used as a marker of ancestral traits provided that the noise introduced by short branches is appropriately dealt with. Based on ancestral dN/dS estimates, we predict that the first cetartiodactyls were relatively small animals (around 20 kg). This finding is in accordance with Cope's rule and the fossil record but could not be recovered via continuous character evolution methods.

MafFilter: a highly flexible and extensible multiple genome alignment files processor

BACKGROUND

Sequence alignments are the starting point for most evolutionary and comparative analyses. Full genome sequences can be compared to study patterns of within and between species variation. Genome sequence alignments are complex structures containing information such as coordinates, quality scores and synteny structure, which are stored in Multiple Alignment Format (MAF) files. Processing these alignments therefore involves parsing and manipulating typically large MAF files in an efficient way.

RESULTS

MafFilter is a command-line driven program written in C++ that enables the processing of genome alignments stored in the Multiple Alignment Format in an efficient and extensible manner. It provides an extensive set of tools which can be parametrized and combined by the user via option files. We demonstrate the software's functionality and performance on several biological examples covering Primate genomics and fungal population genomics. Example analyses involve window-based alignment filtering, feature extractions and various statistics, phylogenetics and population genomics calculations.

CONCLUSIONS

MafFilter is a highly efficient and flexible tool to analyse multiple genome alignments. By allowing the user to combine a large set of available methods, as well as designing his/her own, it enables the design of custom data filtering and analysis pipelines for genomic studies. MafFilter is an open source software available at http://bioweb.me/maffilter.

A fine-scale recombination map of the human-chimpanzee ancestor reveals faster change in humans than in chimpanzees and a strong impact of GC-biased gene conversion

Recombination is a major determinant of adaptive and nonadaptive evolution. Understanding how the recombination landscape has evolved in humans is thus key to the interpretation of human genomic evolution. Comparison of fine-scale recombination maps of human and chimpanzee has revealed large changes at fine genomic scales and conservation over large scales. Here we demonstrate how a fine-scale recombination map can be derived for the ancestor of human and chimpanzee, allowing us to study the changes that have occurred in human and chimpanzee since these species diverged. The map is produced from more than one million accurately determined recombination events. We find that this new recombination map is intermediate to the maps of human and chimpanzee but that the recombination landscape has evolved more rapidly in the human lineage than in the chimpanzee lineage. We use the map to show that recombination rate, through the effect of GC-biased gene conversion, is an even stronger determinant of base composition evolution than previously reported.

A new isolation with migration model along complete genomes infers very different divergence processes among closely related great ape species

We present a hidden Markov model (HMM) for inferring gradual isolation between two populations during speciation, modelled as a time interval with restricted gene flow. The HMM describes the history of adjacent nucleotides in two genomic sequences, such that the nucleotides can be separated by recombination, can migrate between populations, or can coalesce at variable time points, all dependent on the parameters of the model, which are the effective population sizes, splitting times, recombination rate, and migration rate. We show by extensive simulations that the HMM can accurately infer all parameters except the recombination rate, which is biased downwards. Inference is robust to variation in the mutation rate and the recombination rate over the sequence and also robust to unknown phase of genomes unless they are very closely related. We provide a test for whether divergence is gradual or instantaneous, and we apply the model to three key divergence processes in great apes: (a) the bonobo and common chimpanzee, (b) the eastern and western gorilla, and (c) the Sumatran and Bornean orang-utan. We find that the bonobo and chimpanzee appear to have undergone a clear split, whereas the divergence processes of the gorilla and orang-utan species occurred over several hundred thousands years with gene flow stopping quite recently. We also apply the model to the Homo/Pan speciation event and find that the most likely scenario involves an extended period of gene flow during speciation.

Fast and robust characterization of time-heterogeneous sequence evolutionary processes using substitution mapping

Genes and genomes do not evolve similarly in all branches of the tree of life. Detecting and characterizing the heterogeneity in time, and between lineages, of the nucleotide (or amino acid) substitution process is an important goal of current molecular evolutionary research. This task is typically achieved through the use of non-homogeneous models of sequence evolution, which being highly parametrized and computationally-demanding are not appropriate for large-scale analyses. Here we investigate an alternative methodological option based on probabilistic substitution mapping. The idea is to first reconstruct the substitutional history of each site of an alignment under a homogeneous model of sequence evolution, then to characterize variations in the substitution process across lineages based on substitution counts. Using simulated and published datasets, we demonstrate that probabilistic substitution mapping is robust in that it typically provides accurate reconstruction of sequence ancestry even when the true process is heterogeneous, but a homogeneous model is adopted. Consequently, we show that the new approach is essentially as efficient as and extremely faster than (up to 25 000 times) existing methods, thus paving the way for a systematic survey of substitution process heterogeneity across genes and lineages.

Insights into hominid evolution from the gorilla genome sequence

Gorillas are humans’ closest living relatives after chimpanzees, and are of comparable importance for the study of human origins and evolution. Here we present the assembly and analysis of a genome sequence for the western lowland gorilla, and compare the whole genomes of all extant great ape genera. We propose a synthesis of genetic and fossil evidence consistent with placing the human-chimpanzee and human-chimpanzee-gorilla speciation events at approximately 6 and 10 million years ago (Mya). In 30% of the genome, gorilla is closer to human or chimpanzee than the latter are to each other; this is rarer around coding genes, indicating pervasive selection throughout great ape evolution, and has functional consequences in gene expression. A comparison of protein coding genes reveals approximately 500 genes showing accelerated evolution on each of the gorilla, human and chimpanzee lineages, and evidence for parallel acceleration, particularly of genes involved in hearing. We also compare the western and eastern gorilla species, estimating an average sequence divergence time 1.75 million years ago, but with evidence for more recent genetic exchange and a population bottleneck in the eastern species. The use of the genome sequence in these and future analyses will promote a deeper understanding of great ape biology and evolution.

Accelerated evolution and coevolution drove the evolutionary history of AGPase sub-units during angiosperm radiation

BACKGROUND AND AIMS

ADP-glucose pyrophosphorylase (AGPase) is a key enzyme of starch biosynthesis. In the green plant lineage, it is composed of two large (LSU) and two small (SSU) sub-units encoded by paralogous genes, as a consequence of several rounds of duplication. First, our aim was to detect specific patterns of molecular evolution following duplication events and the divergence between monocotyledons and dicotyledons. Secondly, we investigated coevolution between amino acids both within and between sub-units.

METHODS

A phylogeny of each AGPase sub-unit was built using all gymnosperm and angiosperm sequences available in databases. Accelerated evolution along specific branches was tested using the ratio of the non-synonymous to the synonymous substitution rate. Coevolution between amino acids was investigated taking into account compensatory changes between co-substitutions.

KEY RESULTS

We showed that SSU paralogues evolved under high functional constraints during angiosperm radiation, with a significant level of coevolution between amino acids that participate in SSU major functions. In contrast, in the LSU paralogues, we identified residues under positive selection (1) following the first LSU duplication that gave rise to two paralogues mainly expressed in angiosperm source and sink tissues, respectively; and (2) following the emergence of grass-specific paralogues expressed in the endosperm. Finally, we found coevolution between residues that belong to the interaction domains of both sub-units.

CONCLUSIONS

Our results support the view that coevolution among amino acid residues, especially those lying in the interaction domain of each sub-unit, played an important role in AGPase evolution. First, within SSU, coevolution allowed compensating mutations in a highly constrained context. Secondly, the LSU paralogues probably acquired tissue-specific expression and regulatory properties via the coevolution between sub-unit interacting domains. Finally, the pattern we observed during LSU evolution is consistent with repeated sub-functionalization under 'Escape from Adaptive Conflict', a model rarely illustrated in the literature.

Comparing fungal genomes: insight into functional and evolutionary processes

Large amount of genome data are being generated by second- and now also third-generation sequencing technologies. The challenge no longer lies in the generation of the data, but in the analyses of it. We present an overview of approaches and methods to compare complete sequences of related fungal genomes. We focus on evolutionary analyses of genome alignments to describe species divergence and to identify footprints of demography and natural selection within and between species.

Ancestral population genomics

The full genomes of several closely related species are now available, opening an emerging field of investigation borrowing both from population genetics and phylogenetics. Providing we can properly model sequence evolution within populations undergoing speciation events, this resource enables us to estimate key population genetics parameters, such as ancestral population sizes and split times. Furthermore, we can enhance our understanding of the recombination process and investigate various selective forces. We discuss the basic speciation models for closely related species, including the isolation and isolation-with-migration models. A major point in our discussion is that only a few complete genomes contain much information about the whole population. The reason being that recombination unlinks genomic regions, and therefore a few genomes contain many segments with distinct histories. The challenge of population genomics is to decode this mosaic of histories in order to infer scenarios of demography and selection. We survey different approaches for understanding ancestral species from analyses of genomic data from closely related species. In particular, we emphasize core assumptions and working hypothesis. Finally, we discuss computational and statistical challenges that arise in the analysis of population genomics data sets.

The making of a new pathogen: insights from comparative population genomics of the domesticated wheat pathogen Mycosphaerella graminicola and its wild sister species

The fungus Mycosphaerella graminicola emerged as a new pathogen of cultivated wheat during its domestication ~11,000 yr ago. We assembled 12 high-quality full genome sequences to investigate the genetic footprints of selection in this wheat pathogen and closely related sister species that infect wild grasses. We demonstrate a strong effect of natural selection in shaping the pathogen genomes with only ~3% of nonsynonymous mutations being effectively neutral. Forty percent of all fixed nonsynonymous substitutions, on the other hand, are driven by positive selection. Adaptive evolution has affected M. graminicola to the highest extent, consistent with recent host specialization. Positive selection has prominently altered genes encoding secreted proteins and putative pathogen effectors supporting the premise that molecular host-pathogen interaction is a strong driver of pathogen evolution. Recent divergence between pathogen sister species is attested by the high degree of incomplete lineage sorting (ILS) in their genomes. We exploit ILS to generate a genetic map of the species without any crossing data, document recent times of species divergence relative to genome divergence, and show that gene-rich regions or regions with low recombination experience stronger effects of natural selection on neutral diversity. Emergence of a new agricultural host selected a highly specialized and fast-evolving pathogen with unique evolutionary patterns compared with its wild relatives. The strong impact of natural selection, we document, is at odds with the small effective population sizes estimated and suggest that population sizes were historically large but likely unstable.

Detecting coevolving positions in a molecule: why and how to account for phylogeny

Positions in a molecule that share a common constraint do not evolve independently, and therefore leave a signature in the patterns of homologous sequences. Exhibiting such positions with a coevolution pattern from a sequence alignment has great potential for predicting functional and structural properties of molecules through comparative analysis. This task is complicated by the existence of additional correlation sources, leading to false predictions. The nature of the data is a major source of noise correlation: sequences are taken from individuals with different degrees of relatedness, and who therefore are intrinsically correlated. This has led to several method developments in different fields that are potentially confusing for non-expert users interested in these methodologies. It also explains why coevolution detection methods are largely unemployed despite the importance of the biological questions they address. In this article, I focus on the role of shared ancestry for understanding molecular coevolution patterns. I review and classify existing coevolution detection methods according to their ability to handle shared ancestry. Using a ribosomal RNA benchmark data set, for which detailed knowledge of the structure and coevolution patterns is available, I demonstrate and explain why taking the underlying evolutionary history of sequences into account is the only way to extract the full coevolution signal in the data. I also evaluate, using rigorous statistical procedures, the best approaches to do so, and discuss several important biological aspects to consider when performing coevolution analyses.

IMPACT: integrated multiprogram platform for analyses in ConTest

In this paper, we introduce a new Graphical User Interface that estimates evolutionary rates on protein sequences by assessing changes in biochemical constraints. We describe IMPACT, a platform-independent (tested in Linux, Windows, and MacOS), easy to install software written in Java. IMPACT integrates the use of a built-in multiple sequence alignment editor, with programs that perform phylogenetic and protein structure analyses (ConTest, PhyML, ATV, and Jmol) allowing the user to quickly and efficiently perform evolutionary analyses on protein sequences, including the detection of selection (negative and positive) signatures at the amino acid scale, which can provide fundamental insight about species evolution and ecological fitness. IMPACT provides the user with a working platform that combines a number of bioinformatics tools and utilities in one place, transferring information directly among the various programs and therefore increasing the overall performance of evolutionary analyses on proteins.

Estimating divergence time and ancestral effective population size of Bornean and Sumatran orangutan subspecies using a coalescent hidden Markov model

Due to genetic variation in the ancestor of two populations or two species, the divergence time for DNA sequences from two populations is variable along the genome. Within genomic segments all bases will share the same divergence—because they share a most recent common ancestor—when no recombination event has occurred to split them apart. The size of these segments of constant divergence depends on the recombination rate, but also on the speciation time, the effective population size of the ancestral population, as well as demographic effects and selection. Thus, inference of these parameters may be possible if we can decode the divergence times along a genomic alignment. Here, we present a new hidden Markov model that infers the changing divergence (coalescence) times along the genome alignment using a coalescent framework, in order to estimate the speciation time, the recombination rate, and the ancestral effective population size. The model is efficient enough to allow inference on whole-genome data sets. We first investigate the power and consistency of the model with coalescent simulations and then apply it to the whole-genome sequences of the two orangutan sub-species, Bornean (P. p. pygmaeus) and Sumatran (P. p. abelii) orangutans from the Orangutan Genome Project. We estimate the speciation time between the two sub-species to be thousand years ago and the effective population size of the ancestral orangutan species to be , consistent with recent results based on smaller data sets. We also report a negative correlation between chromosome size and ancestral effective population size, which we interpret as a signature of recombination increasing the efficacy of selection.

We present a hidden Markov model that uses variation in coalescence times between two distantly related populations, or closely related species, to infer population genetics parameters in ancestral population or species. The model infers the divergence times in segments along the alignment. Using coalescent simulations, we show that the model accurately estimates the divergence time between the two populations and the effective population size of the ancestral population. We apply the model to the recently sequenced orangutan sub-species and estimate their divergence time and the effective population size of their ancestor population.

Incomplete lineage sorting patterns among human, chimpanzee, and orangutan suggest recent orangutan speciation and widespread selection

We search the complete orangutan genome for regions where humans are more closely related to orangutans than to chimpanzees due to incomplete lineage sorting (ILS) in the ancestor of human and chimpanzees. The search uses our recently developed coalescent hidden Markov model (HMM) framework. We find ILS present in ∼1% of the genome, and that the ancestral species of human and chimpanzees never experienced a severe population bottleneck. The existence of ILS is validated with simulations, site pattern analysis, and analysis of rare genomic events. The existence of ILS allows us to disentangle the time of isolation of humans and orangutans (the speciation time) from the genetic divergence time, and we find speciation to be as recent as 9-13 million years ago (Mya; contingent on the calibration point). The analyses provide further support for a recent speciation of human and chimpanzee at ∼4 Mya and a diverse ancestor of human and chimpanzee with an effective population size of about 50,000 individuals. Posterior decoding infers ILS for each nucleotide in the genome, and we use this to deduce patterns of selection in the ancestral species. We demonstrate the effect of background selection in the common ancestor of humans and chimpanzees. In agreement with predictions from population genetics, ILS was found to be reduced in exons and gene-dense regions when we control for confounding factors such as GC content and recombination rate. Finally, we find the broad-scale recombination rate to be conserved through the complete ape phylogeny.

Base pairing constraints drive structural epistasis in ribosomal RNA sequences

It has long been accepted that the structural constraints stemming from the 3D structure of ribosomal RNA (rRNA) lead to coevolution through compensating mutations between interacting sites. State-of-the-art methods for detecting coevolving sites, however, while reaching high levels of specificity and sensitivity for Watson-Crick (WC) pairs of the helices defining the secondary structure, only scarcely reveal tertiary interactions occurring at the level of the 3D structure. In order to understand the relative failure of coevolutionary methods to detect such interactions, we analyze 2,682 interacting sites derived from high-resolution structures, which include a comprehensive data set of rRNA sequences from Archaea and Bacteria. We report a striking difference in the amount of coevolution between WC and non-WC pairs. In order to understand this pattern, we derive fitness landscapes from the geometry of base pairing interactions and construct neutral networks of substitutions for each type of interaction. These networks show that coevolution is a property of WC pairs because, unlike non-WC pairs, their landscapes exhibit fitness valleys, a single mutation in a WC pair resulting in a fitness drop. Second, we used the inferred neutral networks to estimate the level of constraint acting on each type of base pair and show that it correlates negatively with the observed rate of substitutions for all non-WC pairs. WC pairs appear as outliers, fixing more substitutions than expected according to their level of constraint. We here propose that the rate of substitution in WC pairs is due to coevolution resulting from constraints acting at intermediate levels of organization, namely the one of the helical stem with its forming WC pairs. In agreement with this hypothesis, we report a significant excess of intrahelical, inter-WC pairs coevolution compared with interhelices pairs. Altogether, these results show that detailed biochemical knowledge is required and has to be incorporated into evolutionary reasoning in order to understand the fine patterns of variation at the molecular level. They also demonstrate that coevolutionary analysis provides almost exclusively 2D information and only little 3D signal.