Runaway GC Evolution in Gerbil Genomes

Abstract

Recombination increases the local GC-content in genomic regions through GC-biased gene conversion (gBGC). The recent discovery of a large genomic region with extreme GC-content in the fat sand rat Psammomys obesus provides a model to study the effects of gBGC on chromosome evolution. Here, we compare the GC-content and GC-to-AT substitution patterns across protein-coding genes of four gerbil species and two murine rodents (mouse and rat). We find that the known high-GC region is present in all the gerbils, and is characterized by high substitution rates for all mutational categories (AT-to-GC, GC-to-AT, and GC-conservative) both at synonymous and nonsynonymous sites. A higher AT-to-GC than GC-to-AT rate is consistent with the high GC-content. Additionally, we find more than 300 genes outside the known region with outlying values of AT-to-GC synonymous substitution rates in gerbils. Of these, over 30% are organized into at least 17 large clusters observable at the megabase-scale. The unusual GC-skewed substitution pattern suggests the evolution of genomic regions with very high recombination rates in the gerbil lineage, which can lead to a runaway increase in GC-content. Our results imply that rapid evolution of GC-content is possible in mammals, with gerbil species providing a powerful model to study the mechanisms of gBGC.

Keywords: GC-biased gene conversion; GC-content; biased substitution, fixation bias; gBGC; genome evolution; recombination.

Fig. 1.

GC-skewed synonymous evolution in protein-coding genes located in the known high-GC region. (A) Topology of the phylogenetic relationships between the species analyzed in this study; (B) GC-content in the third-codon position (GC3) per species per gene; (C) the rate of synonymous substitution (dS) per mutational category: weak-to-strong (WS), strong-to-weak (SW), weak-to-weak (WW), and strong-to-strong (SS); (D) comparison between weak-to-strong dS (dS_WS) and strong-to-weak dS (dS_SW).

Fig. 2.

GC-skewed nonsynonymous evolution in protein-coding genes located in the known high-GC region. (A) GC-content in the first and second codon positions (GC12) per species per gene, (B) and nonsynonymous substitution rate (dN) per mutational category: weak-to-strong (WS), strong-to-weak (SW), weak-to-weak (WW), and strong-to-strong (SS), (C) comparison between weak-to-strong dN (dN_WS) and strong-to-weak dN (dN_SW).

Fig. 3.

Outliers in the rate of synonymous substitution (dS) measured from the gerbil-murine split for 8,809 groups of orthologous genes outside the known high-GC region. (A) Topography of the tree of the species included in the alignments of each group of orthologous genes; our focal species include two gerbil species (in pink) and two murine species (in green); dS was measured in the colored branches, that is, from the point of gerbil-murine divergence to the tip of each of the four focal species. (B) Comparison between the rate of weak-to-strong synonymous substitution (dS_WS) and the strong-to-weak synonymous substitution (dS_SW) in the two gerbil and the two murine species, normalized by dividing each value by the average rate of the respective species and mutation category; colored points represent genes above the dS threshold chosen to define outliers (dS >2.5 times the average for the respective species and mutational category); only values under 6 are shown. (C) Euler diagram of the dS_WS outliers, where the area of each segment is approximately proportional to the number of overlapping outlier genes per species (only gene numbers >1 included), shown also as a percentage of the 524 genes that were outliers in at least one species. (D) Pairwise comparisons of the GC-content at the third-codon position (GC3) between the gerbil Psammomys obesus and the mouse Mus musculus for 8,797 groups of orthologous genes (out of 8,809) for which both species have a GC3 measurement.

Fig. 4.

Clustering of outliers in the weak-to-strong dS. (A) Normalized rate of synonymous substitution (dS) for two mutational categories, strong-to-weak (SW) and weak-to-strong (WS), for 8,809 genes; (B) average normalized weak-to-strong dS (dS_WS) in sliding windows of 1 Mb with a step of 0.25 Mb, showing only windows with >3 genes, and with outlying regions marked with a vertical gray bar. The genes are mapped by row to the chromosomes of the Mus musculus reference genome assembly. The horizontal gray line represents the dS threshold chosen to define outlier genes (dS >2.5 times the average for the respective species and mutational category). The dark blue horizontal bars represent the location of the previously known high-GC region.