Both selective and neutral processes drive GC content evolution in the human genome

Data retrieval

Gene and intergenic sequences as well as intron/exon boundaries were obtained from the UCSC genome annotation database [26], assembly hg17. Gene selection was performed as previously described [27]. Isochore boundary coordinates were derived from a previous work [20]. Fine-scale recombination rates and recombination hotspot locations were obtained from the UCSC database; they are based on HapMap Phase I data [28]. Pseudogene sequences and genomic locations derive from Pseudogene.org [29, 30]; only duplicated pseudogenes were selected and genes that generated more than one pseudogene were discarded (this procedure limits the number of observations but avoids multiple ties in statistical analysis). Also we retained only gene-pseudogenes pairs located in the same isochore type (for example, both gene and pseudogene located in isochores H1). The final data set consisted of 364 gene-pseudogene pairs. Duplicated pseudogenes often represent gene fragments; we therefore aligned gene-pseudogene couples using ClustalW [31] and corresponding intron-pseudointron pairs were retained only if they were both longer than 25 bp. Expression data were obtained as previously described [27] and derive from microarray data on 72 healthy human tissues. Mean expression level was calculated as the mean averaged over all tissues (counting as zero all tissues in which there is no detectable expression). Peak expression was calculated as the maximum expression level across all tissues and expression breadth was the number of distinct tissues expressing the gene.

Polymorphism data

Biallelic SNP locations and allele frequencies were downloaded from the HapMap web site [28] (non-redundant dataset, release 21a). Since previous authors [32] have indicated African populations as having genetic variation patterns most compatible with a constant population size, SNP allele frequencies were obtained for Yoruba (YRI), and derive from the genotyping of 60 individuals. The ancestral allele was inferred by alignment with the chimpanzee sequence (UCSC genome browser, assembly panTro1); SNPs were discarded when orthologous chimpanzee regions were unavailable or did not match either human allele. A total of about 2.2 million GC->AT and 1.7 million AT->GC SNPs were retained. We next purged SNPs at CpG sites, as well as those with no associated allele frequencies: the final dataset comprised more than 2 million SNPs.

For the analysis of substitution rates and stationary GC content (GC*), SNPs deriving from the Seattle SNP database [33], which derive from resequencing experiments, were used; for 206 human genes in the Seattle SNP dataset both chimpanzee and macaque orthologous loci could be retrieved.

Data on polymorphic repeat insertions were obtained through the UCSC genome database (RIPs track) and derive from the dbRIP database [34]; RIPs which have been associated with a human genetic disease were discarded. Also, polymorphic insertions were not included in the study if less than 10 instances were described for the same retrotransposon subfamily. Fixed transposon instances were identified and categorized using the UCSC annotation tables that rely on RepeatMasker. Since fixed and polymorphic repeat instances derive from different sources, we verified that no systematic bias occurs in the detection of either insertion events by calculating correlation between polymorphic and fixed chromosomal frequencies; significant correlations were retrieved for Alus, SVAs and L1s (Spearman rho = 0.854, 0.439 and 0.923, respectively; all p values < 0.05). Reference sequences for different retrotransposon subtypes were derived from Repbase Update [35, 36].

Analysis of allele frequency spectra

Introns/intergenic spacers were divided in 1 kb windows (1 kbseqs) starting from the most 5' nucleotide position (with respect to the chromosome orientation) and extending through the intron/intergenic region in 1000 bp non-overlapping steps (residual nucleotides in 3' were discarded). The following features were then calculated (or retrieved) for all 1 kbseqs: (1) fine scale recombination rate, (2) GC content, (3) allele frequencies of comprised SNPs, (4) expression parameters (peak, mean level and breadth) of the corresponding genes. In order to analyze allele frequency spectra after controlling for recombination rate, we applied the following procedure: starting from all 1 kbseqs, we identified couples of 1 kbseqs that differed less than 10% in recombination rate but displayed extremely different GC contents; in particular, we asked one partner of the recombination-coupled 1 kbseqs to be located below the 30^th percentile in the distribution of 1 kbseqs GC content and the other one above the 70^th percentile. This approach yielded two groups of sequences having extremely similar recombination rates (the equality of medians was checked using the Wilcoxon Rank Sum Test) but very different GC contents. A similar procedure was applied to analyze allele frequency spectra after controlling for recombination rate; in particular, 1 kbseq couples were created having similar GC content (a difference lower than 5% was required) but extremely different recombination rates. Again, two groups of sequences were obtained and used for comparisons.

The same approach described above can be extended to control for two variables: for example, in order to compare allele frequency spectra between highly and lowly expressed sequences, 1 kbseqs couples were identified that displayed both similar GC content and recombination rate (less than 5% and 10% difference, respectively) but extremely different expression levels. To allow comparisons between introns and intergenic spacers, percentile values were calculated over the complete set of 1 kbseqs, irrespective of their location.

In order to quantify the displacement of GC vs AT derived allele frequency distributions observed in Quantile-Quantile plots, differences between corresponding percentiles in the two distributions were summed. These measures were used to compare different groups of sequences selected on the basis of relevant variables (for example high and low GC content or recombination rate). We used bootstrapping procedures to assess the statistical significance of differences in allele frequencies shifts. In particular, 1000 permutations were performed and p values were calculated after normality assessment through the Shapiro-Wilk Test.

Multispecies alignments, substitution rates and stationary GC content

Orthologous human-chimpanzee-macaque regions were retrieved using the liftOver utility from UCSC (assemblies: panTro1 and rheMac2) with a cutoff of at least 70% remapping bases. Three-way species alignments were performed using MAVID [37].

In order to calculate substitution rates and GC* after controlling for ancestral GC content or recombination rate, a procedure similar to the one described above for SNP allele frequencies was applied, with the only difference that the inferred ancestral GC content was used instead of human GC content. In particular, 1 kbseqs couples were created (on the basis of either recombination rate or ancestral GC content) and their position subsequently mapped onto the 3-way species (human/chimpanzee/macaque) alignments; at this stage windows containing less than 600 perfectly aligning bases (i.e. the same nucleotide in the 3 species) were discarded and, for the remaining ones, the ancestral sequence was reconstructed by parsimony (only positions where the macaque was identical to either human or chimpanzee were considered).

The number of 1 kbseqs couples and the corresponding number of sites (in MB) that were analyzed for each comparison are reported in table notes. The number of sites does not exactly correspond to the number of sequences multiplied by 1000 because the presence of gaps in the human sequence (as compared to the two primates) can result in alignments longer than 1000 bp.

Substitution rates and stationary GC content were calculated using a previously developed neighbor-dependent substitution model [38, 39]. For each comparison, the two 1 kbseqs groups were then divided in 20 paired sub-samples of equal size; GC* and substitution rates were calculated for each sub-sample; average values are reported in the tables, together with p values obtained from two tailed Wilcoxon Rank Sum Tests for paired samples.

For the analysis of intron GC content in relation to size, we discarded first introns (due to their increased sequence constraints) and introns shorter than 750 bp (in order to spare constrained sequences at splice sites).

For the analysis, of recombination rates in long and short introns, for each gene, two introns were selected so that one was longer than 80^th and the other shorter than 20^th size percentile of introns length distribution. If no introns satisfied the criteria, the gene was not analyzed. Recombination rates were calculated for 500 bp centered around the median position of each intron. Differences in recombination rates were evaluated using the Wilcoxon Rank Sum Test.

For the analysis of fixed variations in recombination hotspots, we selected 897 hotspot on the basis of their size (smaller than 5 kb) and recombination rate (above the 80^th percentile of the distribution of all hotspots); in 790 cases both chimpanzee and macaque orthologous regions could be retrieved.

Statistical analysis

All statistical analysis were performed using R [40]. For loess fittings [41] a smoothing span of 0.5 was used.

Analysis of SNP allele frequencies

The analysis of SNP allele frequencies is a convenient strategy to study GC content evolution for two main reasons. First, when SNP allele frequencies are analyzed, no requirement for base composition stationarity is needed; this is relevant since base composition has been shown not to be at equilibrium in mammals [42, 43]. Second, given the fast evolution of recombination rates and hotspots [44, 45], allele frequencies of SNPs, which represent relatively recent variations, should carry the most evident signature of recombination-associated fixation biases.

Starting from our gene set, we therefore used the chimpanzee sequence to infer the ancestral allele so that variations could be classified as either GC->AT or AT->GC (SNPs at CpG sites were excluded). As previously noted [43] treating SNPs as independent data, despite the extensive presence of linkage disequilibrium in the human genome [42], introduces no bias since linkage is expected to be independent from the GC/AT status of individual SNPs.

BGC and selection are both expected to result in AT->GC mutations segregating at higher frequencies compared to GC->AT. Yet, this effect is expected to be stronger in highly recombining regions if BGC is involved. Conversely, selection on GC content should be acting differentially depending on the background GC content of SNP flanking regions; in particular, AT->GC variations are expected to segregate at higher frequency in GC rich regions, irrespective of recombination rate. In order to disentangle a possible selective effect from BGC, we analyzed SNP allele frequencies in noncoding genomic sequences after correcting for either GC content or recombination rate. As further detailed in methods, genomic regions were divided in 1 kb sub-sequences. These latter were arranged in couples having very similar recombination rate and extremely different GC content for the comparison of allele frequencies between GC-rich and -poor sequences. Similarly, for the comparison between high- and low-recombining sequences, sequences were arranged in couples showing very similar GC content but extremely different recombination rates. The results of SNP frequency spectra analysis are reported in figure 1 as quantile-quantile plots; in agreement with previous findings [14] and consistent with the action of BGC, GC derived alleles display higher frequencies than AT alleles but the effect is significantly (p < 10^-17) stronger in highly recombining regions for both introns (Figure 1A) and intergenic spacers (see Additional file 1); yet, when allele frequencies were compared after fixing recombination rate, a residual effect of GC content was observed: derived GC alleles segregate at significantly higher frequencies in regions showing a high GC content (p = 2.07 × 10^-10) compared to AT-rich regions (Figure 1B and Additional file 1 for intergenic spacers).

These data suggest that GC content or other related features affect SNP allele segregation independently of recombination rates, although we cannot formally rule out the possibility that extinct recombination hotspots have played a role in the allele frequency spectra we observe. Indeed, as reported above, recombination hotspots are fast evolving [44, 45] and, therefore, the observed increased segregation of GC alleles in GC-rich regions might have been caused by the presence of an hotspot which is now inactive. Yet, if this latter were the case, given the relatively small effect that recombination has played in recent primate history on GC variation (see below), and given that most SNPs are specific to humans (and, therefore, relatively young), a direct role for GC content in promoting recombination events must be postulated to explain our results.

Since GC content has been shown to increase transcriptional activity [46] and some authors detected a positive correlation between gene expression parameters and GC content [16, 18, 19, 46], we wished to determine whether expression level, rather than GC content per se, was responsible for increased segregation of GC alleles. Yet, after controlling for both GC content and recombination rates (as described in methods, we used a similar approach to the one described above) we detected no significant difference in SNP allele frequencies between high- and low-level expressed genes (Figure 1c and Additional file 1). These data are not consistent with selection acting on highly or broadly expressed human genes to increase (or maintain) their GC content, although we cannot exclude that such a selection has acted during vertebrate evolution and subsequently relaxed in humans (further data on gene expression level and GC content evolution are reported below).

Analysis of retrotransposon insertion polymorphisms

It has been suggested [3] that, if selection is acting on base composition, it should also affect the fixation probabilities of transposable elements; indeed, fixed Alu and LINE-1 (L1) elements (average GC content of reference sequences = 0.53 and 0.41, respectively) are differentially represented in the human genome depending on GC content [2], despite both having a preference for AT-rich integration sites [47, 48]. Yet, a better estimation of fixation versus integration probabilities might be obtained by the comparison of polymorphic and fixed transposable elements. We retrieved all available instances of retrotransopson insertion polymorphisms. For both Alu and L1 repeats, we restricted the analysis of fixed repeats to the same subfamilies showing at least 10 instances of polymorphic insertions. Such subfamilies represent relatively young insertion events, yet, given the previously reported preference of older Alus for GC rich regions [2, 49] we further purged all fixed Alu elements showing a divergence higher than 5%. As shown in figure 2A, both SVA (GC content of reference sequence = 0.63) and Alu fixed elements are located in regions with significantly higher GC content compared to their polymorphic counterparts (Wilcoxon Rank Sum Test, two-tailed, p = 0.028 and p < 10^-21, respectively); conversely, fixed and polymorphic L1 flanking regions do not show different average GC contents, being both relatively GC poor (as L1 sequences are).

We next wished to verify whether polymorphic and fixed repetitive elements were differently distributed depending on isochore type. Isochores were classified according to a recent [20] description and are referred to as L1, L2, H1, H2, and H3, in order of increasing GC levels. The results of transposable element distribution are reported in figure 2B and indicate that fixed Alus are significantly enriched (Chi Square Test, p < 10^-5) within heavy isochores compared to polymorphic instances, while no different isochore distribution of fixed vs polymorphic repeats was evident for SVAs (possibly because of the small number of polymorphic insertions, n = 60) or L1s. For further confirmation we performed this same analysis using IsoFinder isochores [50] and the same results were obtained (see Additional file 2). These data confirm the preferential integration of Alus and L1s in AT-rich regions (polymorphic L1 and Alu distributions are relatively similar, Figure 2B), but indicate that additional forces, which relate to GC content, drive their fixation.

We believe that our results differ from previous reports showing no different GC content surrounding polymorphic and fixed Alus [51, 52] because of the larger sample of polymorphic elements we analyzed.

Analysis of substitution rates and stationary GC content

Still, the data we report here are consistent with selection acting to maintain GC content but also with the presence of mutation biases operating in different GC content regions. In order to evaluate this latter possibility we calculated substitution rates and GC* using either fixed variations or SNPs; while SNPs can reasonably be thought to reflect mutation rates, fixed variations depend on both mutation rates and fixation probabilities. In this case, in order to avoid biases towards high frequency variants, the analysis was restricted to intronic regions deriving from 206 fully resequenced genes (see methods). Also, given the influence, documented above, of recombination on mutation rates, we used only gene regions (1 kb windows) showing low crossover rates.

Finally, we wished to verify whether analysis of substitution rates and GC* confirmed our above indication that gene expression levels have not been influencing base composition evolution in recent human history. In addition to serving as a useful confirmation, this approach allows analysis of fixed variations at CpGs, which is not feasible using SNP allele frequency spectra (due to recurrent mutations at these dinucleotides); this is relevant to the topic we are addressing since previous authors have indicated that both gene GC content and CpG level correlate with gene expression parameters [19]. Again, we analyzed substitution rates and GC* in genes displaying narrow and wide expression breadth, after controlling for both GC content and recombination rates: we found no significant differences in either substitution rates (including CpG->TpG) or GC* between the two groups of sequences (not shown).

Local excess of AT->GC fixed variations at recombination hotspots

Intron GC distribution deviates from neutral expectations

Finally, we wished to determine whether GC content in human introns conforms to neutral expectations. As shown in figure 3A, human introns located in light and heavy isochores yield two relatively distinct distributions when their GC content is plotted against size; the effect is not due to the presence of transposable elements, since similar trends are observed when GC content is calculated after masking repetitive sequences (see Additional file 4). Yet, a similar relationship is somehow expected: shorter introns are likely to display more extreme GC values due to sampling biases. In order to verify that this is not the sole explanation for our findings, for each intron and after masking for repetitive sequences, we calculated the GC content in a 200 bp window (GC₂₀₀) centered around its median position. In particular, only introns longer than 500 bp were analyzed (in order to avoid splice site constraints) and GC content was calculated only if the 200 central nucleotides were covered by repeats for less than 20% of their sequence. This procedure assures that the same number of intronic nucleotides is used for GC content calculations so that sampling biases (due to extreme variations in intron size) are avoided. Introns from either light or heavy isochores were then grouped in 6 percentile size classes and their GC₂₀₀ analyzed: a significant decrease of GC₂₀₀ with increasing residual intron size (intron size calculated after repeat removal) is observed for introns located in heavy isochores while an increasing trend is evident for those located in GC-poor isochores (Figure 3B).

We speculated that these results might originate from the preferential location of genes with short introns in regions displaying an extreme GC content. Yet, we verified that this is not the case, since introns belonging to the same gene tend to recapitulate the distributions observed above; in particular, introns belonging to genes located in light isochores tend to display an increase in GC content with size; those located in heavy isochores behave in the opposite manner. This is shown in figure 3: we selected genes having more than 15 introns and calculated, for each one, the correlation coefficient between the masked GC content of its intervening regions and their residual size; the distributions of correlation coefficients are shifted to positive and negative values for genes located in light and heavy isochores (Figure 3C), respectively. The significance of this finding was assessed by re-sampling (GC content and intron size were randomly assorted 1000 times for each gene).

All these analyzes have been performed after removal of transposable elements from both GC and size calculations; still, it might be argued that old, unrecognizable transposable elements have contributed to both intronic GC content and size, therefore explaining the observed distributions. In order to verify that this is not the sole explanation for our findings, we analyzed nonrepetitive GC content and residual intron length in intron-pseudointron pairs: old transposable elements gave the same contribution to both intron and pseudointrons (as their insertion predated pseudogene duplication) and therefore, once recognizable transposable elements have been masked, any difference in GC distribution is expected to be accounted for by repeat-independent events. Data are reported in figure 3D and show the homogenization of GC content in short pseudointrons (compared to real ones) located in light or heavy isochores.

It should be noted that many different isochore-identification methods have been described. We therefore verified that the results above were also obtained using IsoFinder [50] isochore definition (see additional file 4 for figures and details); also, the same results are obtained when the gene GC content (rather than isochore attribution) is used to define "light" (average GC content < 0.41) and "heavy" genes (see additional file 4 for figures and details).

In summary, these data indicate that intron GC content and size do not evolve independently; even when possible confounding effects such as size variation, presence of transposable elements and skewed genomic location are taken into account, isochore-specific correlations exist between intron size and GC content. Although there is no theoretic basis to expect it, we verified that no significant difference exists between recombination rates of long and short introns in both heavy and light isochores (not shown, see methods for details). Therefore, the data we report here can hardly be reconciled with a vision whereby BGC alone drives GC content evolution; rather, these finding might be consistent with a role of both base composition and intron size in gene regulation mediated by nucleosome positioning or chromatin conformation, as previously proposed [18, 23]. In agreement with this view, it has recently been shown [55] that a considerable amount of human intronic sequence is weakly selected, possibly due to its functioning in chromatin structure and transcription regulation.