BACKGROUND
the primary processes that contribute to nucleotide substitution rates are mutation, selection, and population history, but their relative contributions can vary substantially among genes and genomic regions.
for example, selection varies across genes as a consequence of protein function and gene expression patterns  <cit> .
similarly, mutation rates vary across genomic regions as a function of base composition and recombination rate  <cit> .
population history may be major determinant of substitution rate in the presence of selection, but should not be a factor in the absence of selection  <cit> .
to date, our understanding about the evolutionary forces that contribute to nucleotide substitution rates has been based primarily on the study of coding regions.
inferences based on coding data reflect, in large part, the action of selection.
an obvious example is variation in nonsynonymous substitution rates among genes, which is determined primarily by differential selective constraint.
a subtler example is substitution rates at third codon positions, which are a function of both mutation and selection on codon usage.
highly biased genes evolve more slowly at synonymous sites  <cit> , but synonymous substitution rates are also correlated with gc content  <cit> .
the important point is that it can be difficult to disentangle the contribution of selection and mutation to rate variation among coding regions.
in contrast, non-coding regions should be under little selective constraint, and thus nucleotide substitution in these regions should be governed primarily by neutral processes like mutation.
studying non-coding regions can be difficult in practice, however, because orthology is not always clear.
one way to identify orthologous non-coding regions is to compare, for example, non-coding regions that are 5' and 3' to orthologous genes  <cit> .
the problem is that these regions are also expected to be enriched for functional elements, like promoter and enhancer sequences, and thus potentially under strong selection.
another way to identify orthologous regions is to compare non-coding regions among species via similarity  searches, but this approach is also likely to enrich for slow-evolving regions under selective constraint.
indeed, selective constraint on non-coding regions may be more pervasive in eukaryote genomes than previously assumed  <cit> .
for example, there is compelling evidence from drosophila and several mammalian species that some non-coding regions evolve more slowly than synonymous sites  <cit> , presumably due to selective constraint on non-coding nucleotide substitutions  <cit> .
thus, if non-coding data are not chosen carefully, it can be as difficult to disentangle the relative contribution of mutation and selection for non-coding data as it is for coding data.
to date, there have been few studies comparing evolutionary rates among plant nuclear genes  <cit> .
the dearth of studies reflects, in part, difficulties substantiating orthology relationships, which are complicated because plant genomes contain more duplicated genomic regions than animal genomes  <cit> .
this orthology problem is magnified for non-coding regions, which may evolve and rearrange more rapidly than coding regions.
as a result, rates and patterns of sequence evolution among plant non-coding regions have not been characterized in any detail.
here we contrast nucleotide substitution rates between orthologous coding and intergenic regions of arabidopsis thaliana and arabidopsis lyrata, two plant species that diverged ~ <dig> million years ago  <cit> .
their genomes are largely collinear, but they differ in chromosome number , in dna content , and by several translocations  <cit> .
the two species also differ in population history; a. thaliana is predominantly selfing, while a. lyrata is an obligate outcrosser.
these differences should affect differences in patterns and rates of molecular evolution between species, but the expected differences are not readily apparent  <cit> , perhaps because a. thaliana only recently shifted to a selfing mating system  <cit> .
to study rate variation in intergenic regions, we have generated sequence data in a. lyrata using pcr primers that flank remnants of transposable elements  in a. thaliana.
we reason that these regions are non-functional by virtue of te insertion, and thus comprise a data set that should be relatively free of selective constraint.
we verify their orthology both with comparative maps of the two species and by their molecular features.
the intergenic data are contrasted to a second data set consisting of large  exonic sequences from a. lyrata and a. thaliana.
with these two data sets, we address several questions about arabidopsis nucleotide substitution rates, such as: i) do intergenic sequences evolve at rates similar to synonymous sites in coding data?
ii) do any genomic features, like gc content or recombination, correlate with nucleotide substitution rate variation among loci?
iii) what can be inferred about the relative contribution of mutation and selection to nucleotide substitution?
and iv) do intergenic regions provide any hints to the mechanisms that contribute to genome size differences between a. lyrata and a. thaliana?
RESULTS
isolation and location of orthologous intergenic regions
we identified te remnants in the a. thaliana genome and designed pcr primers to flank  <dig> of these remnants.
three primers were designed: two flanked the te remnant, and a third  primer was specific to the te.
we attempted amplification in a. lyrata with two separate pcr reactions for each of the  <dig> intergenic regions.
the first used the two flanking primers, and the second used one flanking primer with the internal te-specific primer.
these two primer sets were also applied to a. thaliana ecotype colombia as a positive control .
the intergenic regions were distributed across the a. thaliana genome and ranged in size from 200– <dig> bases.
of  <dig> attempted a. lyrata pcr amplifications, 34%  were successful.
amplifications were considered a success when either of two cases occurred.
in case one, the flanking primers amplified a band near the expected size and the internal primer also amplified a band of the expected size.
the second successful case was when the flanking primers amplified a band consistent in size with the absence of the te, and the internal primer failed to amplify any band.
pcr failure was usually a failure of both primer combinations, but ~10% of pcr "failure" was unexpected pcr results.
examples include patterns in which flanking primers amplified a band greatly different than the expected size or in which flanking primers amplified a small band consistent with the loss of the te but amplification with the internal "te" primer was successful.
in order to understand sampling biases associated with data generation, we mapped the location of pcr successes and failures along the genome to determine whether successful pcr amplifications were clustered with respect to their relative positions on the a. thaliana genome.
we coded successful and unsuccessful pcrs as 1's and 0's, respectively, and created a binary string to represent the linear order of pcr results.
we then calculated the variance in the number of zeros  bounded by ones  or the end of chromosomes.
the observed variance in our data was compared to a distribution of variances based on  <dig>  random shufflings of the binary string.
the observed variance was much greater than  <dig> % of simulated variances, indicating that successful pcr amplifications were significantly clustered .
this clustering was related, in part, to gene density.
we found that attempted amplifications were more likely to be successful in low gene density regions of a. thaliana , when we calculated gene density in a window of  <dig>  mb around the a. thaliana te.
consistent with this observation, all attempted amplifications were successful in regions with gene density lower than  <dig> annotated genes per  <dig>  mb, although there were only  <dig> amplifications in regions with density this low .
this effect seemed to be relatively local, however, because larger  <dig>  mb windows retained the basic trend but the trend was no longer significant .
many of the amplicons were small, presumably due to the absence of the te in a. lyrata.
we cloned and sequenced a subset of  <dig> amplicons, yielding  <dig> alignable sequences.
in order to help establish orthology, we examined the position of these  <dig> amplicons with respect to comparative maps of a. thaliana and a. lyrata  <cit> .
we mapped each amplicon to the a. thaliana genome and found its place between markers that had been mapped on a. lyrata.
we considered an amplicon as belonging to a collinear region of the genome if the two markers flanking the marker were also neighboring markers in a. lyrata.
by this criterion, we could assign  <dig> of  <dig>  to unambiguously collinear regions.
another  <dig> were ambiguous either because they had only one flanking marker  or because a third marker had been rearranged between the two markers on the a. lyrata map.
only one amplicon was clearly in a non-collinear region.
although the resolution of comparative a. thaliana-a.
lyrata maps is limited and the possibility of paralogy due to segmental duplication cannot be dismissed entirely, most  of our intergenic regions were orthologous by the criterion of collinearity.
the dna sequence size of our  <dig> intergenic regions ranged from 218– <dig> bp, with an average length of  <dig> bp .
none of the sequences had extensive open reading frames; the longest orf was  <dig> bp.
based on the arabidopsis small rna project  <cit>  targetfinder, there were also no obvious conserved small rnas present in the sequences.
of the  <dig> regions,  <dig> a. lyrata sequences contained a homolog of the te remnant found in a. thaliana.
these  <dig> sequences were alignable both in the te remnant and in regions flanking the te remnant, providing further molecular evidence for orthology.
the remaining  <dig> a. lyrata sequences did not contain the te remnant, but the sequences were homologous to the regions flanking the a. thaliana te.
altogether, 91%  of our amplicons were either in a collinear region or contained the te remnant with associated flanking regions, providing strong evidence that our amplicons are orthologs.
substitution rates in intergenic vs. synonymous sites
given orthologous sequences, we estimated evolutionary distances for intergenic loci.
we report distances based on the k2p model  , but all results and analyses are qualitatively identical using the general time-reversible model .
the average dk2p  estimate over all  <dig> regions was  <dig>  substitutions per site .
the dk2p estimates ranged ~5-fold from  <dig>  to  <dig>  substitutions per site.
these numbers of course report the distances as if they were sampled from a single distribution, and a priori it might be thought  that sequences with tes evolve at different rates than those without.
there is, however, no evidence to this effect .
similarly, we examined whether evolutionary rates differed by collinearity, reasoning that paralogous sequences could be represented in the group for which collinearity was ambiguous.
there was no difference between groups based on collinearity .
furthermore, the six amplicons for which there was no evidence for orthology from either te presence or collinearity fell well within the extremes of the dk2p range .
because there is no obvious evidence for rate classes based on te presence or collinearity, we treated the  <dig> intergenic sequences as a single group.
in addition to our non-coding data set of  <dig> sequences, we generated an additional set of  <dig> large  exons.
these exons were distributed throughout the genome .
forty of these  <dig> genes are located in non-ambiguously collinear regions between a. thaliana and a. lyrata.
of the remaining  <dig>   <dig> had ambiguous collinearity and one was in a non-collinear region.
we calculated synonymous distances  between sequences for all  <dig> loci.
the  <dig> genes with uncertain collinearity did not differ from the remaining  <dig> loci either as group , but two of these exons had the highest ds values in our sample, with one value >  <dig> -fold higher than that of any other exon .
we chose to remove these two loci from further analyses because they may represent paralogs, although none of our overall results vary with their inclusion.
based on the remaining  <dig> exons, mean estimates of ds were  <dig>  substitutions per synonymous site, with a range from  <dig>  to  <dig>  .
our principle goal is to compare substitution patterns and rates between coding and intergenic regions.
loci within the two data types varied ~5-fold in genetic distance.
the coefficients of variation, which was  <dig>  for intergenic data and  <dig>  for coding data, indicate that distances vary to similar degrees for the two data types.
however, the mean genetic distance differed between the two classes , with synonymous sites evolving more rapidly than non-coding sites .
mean distances can be used to estimate nucleotide substitution rates.
assuming that the two lineages diverged roughly  <dig> million years ago  <cit> , the mean substitution rate for synonymous sites was  <dig>  × 10- <dig> substitutions per site per year.
in contrast, the mean non-coding rate was only two-thirds of that value, at  <dig>  × 10- <dig> substitutions per site per year .
genomic correlates with substitution rates
our results indicate that mean genetic distances vary between data types and also that genetic distances vary among loci.
we also sought to determine whether genetic distances are related to genomic features such as gc content and recombination rate.
we expect that such relationships will be clearer in intergenic data, because the patterns are less likely to be complicated by selective forces.
recombination may contribute to rate variation among loci by introducing mutations  <cit> .
recombination rates have been estimated for a. thaliana, based on comparisons of physical  and genetic  distances  <cit> , but there are no direct estimates of cm/mb recombination rates for a. lyrata.
we therefore used a. thaliana recombination rates to investigate relationships between recombination rate and genetic distances.
there was no strong correlation between recombination rate and genetic distance whether the data were combined  or examined separately as intergenic and genic.
the lack of correlation between recombination rates and genetic distances parallels a previous study that detected no correlation between recombination rates and a. thaliana sequence polymorphism  <cit> .
however, nordborg et al.
did detect a negative correlation between polymorphism and gene density.
given this observation, we too examined the relationship of genetic distance with gene density.
gene density is negatively correlated with substitution rate in both datasets, although this relationship is only significant in the intergenic data  and not in coding data  .
these results were generated with gene density measured as the number of genes in a  <dig>  mb window centered around each intergenic locus; the result is similar, but only borderline significant, with  <dig>  mb windows .
gc content is another factor known to correlate with substitution rate.
the gc content of our intergenic sequences did not differ markedly between a. lyrata and a. thaliana .
genic regions also had similar gc content between the two species, but average gc values for coding dna were  <dig> % higher than intergenic regions .
this contrast also held when only synonymous sites in coding regions were considered .
the observed transition:transversion ratio was also higher for coding data  than for intergenic regions , although statistical support for this trend was borderline .
cpg deamination is a common explanation for elevated transition:transversion ratios in high gc regions and could cause differences in rates between genic and intergenic regions.
we thus recalculated divergence, treating all cpg dinucleotides as non-variable in both coding and intergenic data sets.
although there were more cpg sites in genic regions, the cpg disparity alone did not explain the difference in rate between genic and intergenic regions.
genic synonymous divergence rates without cpg sites was  <dig>  × 10- <dig> substitutions per site per year, while intergenic divergence rates were  <dig>  × 10- <dig> substitutions per site per year.
these rates are still significantly different from each other .
with combined intergenic and genic data, there was a strong positive correlation between genetic distance and gc content .
this correlation also held for synonymous sites alone .
to better investigate the relationship between gc and genetic distances, we performed an analysis of covariance.
the ancova examined the effect of gc on genetic distance, with sequence type as a factor to determine if gc content alone accounted for the differences in genetic distance between our coding and intergenic data .
the ancova model indicated that gc content is a significant predictor for genetic distance , but also that there is an additional effect of sequence type .
patterns of indel variation
intergenic regions do not have constraints on coding frame, and thus accumulate indels.
these indels, in turn, may provide insights into processes that contribute to the  <dig> × size difference between a. lyrata and a. thaliana genomes.
we analyzed indel size distributions in intergenic data.
for these analyses, we ignored sequences in which the remnant te was present in a. thaliana but absent in a. lyrata.
we also ignored gaps at sequence ends.
for the purpose of clarity, we denoted the species with the non-gapped sequence in an alignment as containing an 'insertion'.
our intergenic data contained  <dig> distinct a. lyrata insertions totaling  <dig> bases.
a. thaliana, had more distinct insertions  but fewer inserted bases .
the mean size of insertions was  <dig>  bp in a. lyrata and  <dig>  bp in a. thaliana, and did not differ significantly.
in general, a. lyrata insertions tended to be longer and fewer .
although the distributions of insertion sizes appear somewhat different, with the a. lyrata distribution having a longer tail of large insertions and a. thaliana having a higher proportion of small insertions, statistical tests comparing the distributions of the samples were not significant .
furthermore, a sign test comparing the relative length of a. thaliana and a. lyrata sequences across loci was not significant , and thus there is no evidence that intergenic regions are systematically longer in a. lyrata.
discussion
we have generated a sample of orthologous intergenic sequences from a. lyrata to compare to a. thaliana.
these sequences were designed to include remnants of transposable elements , with the a priori thought that these regions are not under strong selection and therefore that the distribution of evolutionary rates provides insight into neutral genomic processes.
nonetheless, our analyses of the genomic location of pcr successes reveal subtle biases in our non-coding data that are important to keep in mind when discussing nucleotide substitution rates.
first, the data are not random with respect to genomic location and thus may be biased toward more slowly evolving regions of the genome.
second, they tend to originate from regions of low gene density.
this does superficially suggest, however, that they are not enriched for functional elements like enhancers or promoters.
third, most of the sequences originate from regions of genome collinearity.
coupled with molecular features, there is strong evidence that our final data set of  <dig> non-coding regions represents orthologous dna features.
this set of non-coding sequences represents one of the few – and perhaps only – multi-locus, intergenic divergence data in plants.
we compared our intergenic data to  <dig> long exons that were sequenced in a. lyrata.
the primary result is that the mean rate of intergenic nucleotide substitution is two-thirds that of the synonymous coding data, with an absolute rate estimated to be  <dig>  × 10- <dig> substitutions per site per year.
this result holds with alternative nucleotide models , and thus does not appear to be solely an issue of estimation procedures.
slower rates in non-coding regions relative to synonymous sites are becoming a surprisingly frequent observation.
for example, a recent study of drosophila demonstrated that non-coding dna evolves considerably slower than synonymous sites in terms of both divergence between species and polymorphism within species  <cit> .
by comparing studies, one can also make the case that pseudogenes  <cit>  and introns  <cit>  evolve more slowly than synonymous sites in apes and other mammals  <cit> .
studies of mammalian intergenic regions have also found slower rates than synonymous sites  <cit> .
although most of these studies encompass only a handful of genes, an overall picture of relatively slow non-coding rates is emerging.
why are non-coding rates slow relative to synonymous rates?
one potential reason discussed above  is methodological biases.
our observations were consistent across three different alignment methods and were qualitatively similar when we performed the analyses on synonymous and four-fold degenerate sites using the same k2p model as we applied to intergenic data.
thus, neither model choice nor alignment appears to contribute substantially to our inferences.
nonetheless, one cannot fully discount sampling biases in the noncoding data, particularly because only 34% of our intergenic pcr amplifications were successful.
some of this failure is attributable to the fact that some tes, particularly basho elements, do not appear to be shared between the two species and thus could not be amplified in a. lyrata.
we attempted to evaluate the magnitude and effect of pcr bias by performing a parallel experiment on coding data.
we designed three primers to amplify each of the  <dig> exons, with the outside primers anchored in flanking intronic or utr regions and an internal primer specific to the exon .
to mimic our amplifications of noncoding dna, we applied the same pcr procedures and same criteria as we applied to noncoding data .
although  <dig> of  <dig> amplifications were successful in a. thaliana, only  <dig> of  <dig>  were successful in a. lyrata, roughly mimicking the success rate  of noncoding amplifications.
the group of exons that were successfully amplified by this method did not differ from "unsuccessful" loci in mean genetic distance  has concluded that non-coding regions of drosophila are under selective constraint and also subject to bouts of adaptive selection.
it seems unlikely that this is a general phenomenon, however, particularly for plants.
plant genomes change rapidly in size, in large part due to the accumulation of repetitive dna  <cit> .
it is thus difficult to imagine that this repetitive dna is under strong selective constraint.
arabidopsis species may be an exception because they have relatively little repetitive dna for plant genomes.
nonetheless, there is little evidence for selection on non-coding regions that flank coding dna  <cit> .
further, the selective interpretation is particularly difficult to argue in this case because the intergenic regions in our study were a priori chosen for their apparent lack of function and tend to originate from genomic regions of low gene density.
consistent with the possibility of low functionality, the intergenic data were replete with indel variation.
a previous study found that intron lengths from  <dig> genes differed significantly between these two species  <cit> , suggesting that the  <dig> × difference in genome size between a. lyrata and a. thaliana is due, at least in part, to the accumulation of small sequence changes.
surprisingly, intergenic indel patterns and lengths do not differ substantially between a. lyrata and a. thaliana, and we find no evidence that differential indel events contribute to differences in genome size.
the third explanation for relatively slow rates of evolution in non-coding dna is differential mutation rates.
recombination could contribute to this mutational effect, because recombination can cause mutations during the resolution of double-strand breaks  <cit> .
if recombination and mutation are associated processes and most mutations are neutral, then both polymorphism and divergence should be positively correlated with recombination rate.
however, we find no correlation between genetic distances and a. thaliana recombination rates .
similarly, recombination does not correlate with polymorphism in ~ <dig> a. thaliana gene fragments  <cit>  or  <dig> a. lyrata genes  <cit> .
it is possible, of course, that the lack of correlation between evolutionary rate and recombination is a false negative.
a true positive correlation might go undetected if recombination rate estimates are imprecise  or if there have been shifts in genomic patterns of recombination between a. thaliana and a. lyrata.
to date, there is little evidence for the latter  <cit> .
overall, we have no evidence to suggest that recombination rates contribute to the differences in rates between intergenic and genic dna or to the observed ~5-fold variation in evolutionary rates among loci.
in contrast to recombination rates, there is a strong relationship between genetic distance and gc content in our data .
substitution rates have long been known to correlate positively with gc content  <cit> , presumably due to higher mutation rates in high gc regions, in part due to cpg effects  <cit> .
consistent with these previous observations, the synonymous sites of our coding data have a high average gc content relative to intergenic data, faster evolutionary rates, and a higher proportion of transitional  changes.
interestingly, gc could contribute to some of the non-coding effects observed in drosophila, too.
there is a striking difference in gc content between the coding and non-coding regions analyzed by andolfatto ; non-coding sequences have a mean gc content of 43%, compared to 58% gc in coding regions.
this disparity alone, rather than selective constraint, may account for the slow rate of divergence in drosophila non-coding dna.
gc content alone cannot explain, however, the skewed ratios of polymorphism to divergence found by andolfatto , unless there has been a recent shift in gc mutation biases.
nonetheless, gc content does not provide a complete explanation for evolutionary rate variation in our data.
one of the most interesting aspects of genome-wide analyses in a. thaliana has been the negative correlation between gene density and polymorphism in a. thaliana  <cit> .
this correlation has been interpreted to be the consequence of either selective sweeps or, more likely, background selection  <cit> .
we have uncovered a similar negative correlation between gene density and evolutionary rate for intergenic data , providing a potential neutral explanation for the observation – i.e., that mutation rates are higher in regions of low gene density.
note also that, if anything, our intergenic data are biased to low gene density regions of the genome where hitchhiking and background selection should not be particularly strong.
the negative correlation between genetic distance and gene density remains borderline significant with partial correlations that also consider gc content .
thus, the gene density effect does not appear to be solely an issue of gc, but the causes of this effect are elusive at this point.
more importantly, there is a "sequence-type" difference between genic and intergenic substitution rates that is not accounted for by differences in gc content .
differential selection of genomic regions is unlikely to explain the observed differences between sequence types, as we would expect the substitution rate effect to be in the opposite direction .
it is possible, though, that mutation varies between intergenic and coding regions.
for example, transcription-related mutations could increase synonymous substitution rates in coding regions over and above gc effects.
this possibility seems plausible because a. thaliana base composition varies as a function of gene expression, suggesting that mutation rates among genes scale with transcription rates  <cit> .
if this hypothesis is correct, then introns should evolve at a rate that is more similar to exons than to intergenic regions.
to examine this hypothesis, we aligned the  <dig> available a. lyrata introns from genbank to their a. thaliana ortholog.
the mean distance  between intron sequences was  <dig>  substitutions per site.
this distance does not differ significantly from our synonymous divergence estimates , but it is significantly higher than the mean genetic distance at our  <dig> intergenic loci .
thus, although the data are limited, intronic sequences are consistent with the hypothesis that transcription-related mutation contributes to differential substitution rates between exons and intergenic regions.
CONCLUSIONS
it is clear that gc content is a major determinant of evolutionary rate variation, not only between sequence types  but also among loci.
gc content may contribute to a. thaliana polymorphism, too, because there is a positive, borderline significant correlation  between gc content and snp polymorphism  for  <dig> polymorphic intergenic loci in the nordborg et al.
data panel of  <dig> individuals .
on the other hand, gc content does not fully explain intergenic rate variation; variation also correlates with gene density after correcting for gc content.
more importantly, there is a detectable effect of sequence type.
the explanation that we deem most reasonable for contributing to the sequence-type effect is transcription-related mutation.
transcription-related mutation predicts a pattern of higher evolutionary rates in transcribed regions, and this prediction is upheld with both exonic and intronic data.

 the primary processes that contribute to nucleotide substitution rates are mutation, selection, and population history, but their relative contributions can vary substantially among genes and genomic regions. for example, selection varies across genes as a consequence of protein function and gene expression patterns. similarly, mutation rates vary across genomic regions as a function of base composition and recombination rate. population history may be major determinant of substitution rate in the presence of selection, but should not be a factor in the absence of selection. to date, our understanding about the evolutionary forces that contribute to nucleotide substitution rates has been based primarily on the study of coding regions. inferences based on coding data reflect, in large part, the action of selection. an obvious example is variation in nonsynonymous substitution rates among genes, which is determined primarily by differential selective constraint. a subtler example is substitution rates at third codon positions, which are a function of both mutation and selection on codon usage. highly biased genes evolve more slowly at synonymous sites, but synonymous substitution rates are also correlated with gc content. the important point is that it can be difficult to disentangle the contribution of selection and mutation to rate variation among coding regions. in contrast, non-coding regions should be under little selective constraint, and thus nucleotide substitution in these regions should be governed primarily by neutral processes like mutation. studying non-coding regions can be difficult in practice, however, because orthology is not always clear. one way to identify orthologous non-coding regions is to compare, for example, non-coding regions that are 5' and 3' to orthologous genes. the problem is that these regions are also expected to be enriched for functional elements, like promoter and enhancer sequences, and thus potentially under strong selection. another way to identify orthologous regions is to compare non-coding regions among species via similarity searches, but this approach is also likely to enrich for slow-evolving regions under selective constraint. indeed, selective constraint on non-coding regions may be more pervasive in eukaryote genomes than previously assumed. for example, there is compelling evidence from drosophila and several mammalian species that some non-coding regions evolve more slowly than synonymous sites, presumably due to selective constraint on non-coding nucleotide substitutions. thus, if non-coding data are not chosen carefully, it can be as difficult to disentangle the relative contribution of mutation and selection for non-coding data as it is for coding data. to date, there have been few studies comparing evolutionary rates among plant nuclear genes. the dearth of studies reflects, in part, difficulties substantiating orthology relationships, which are complicated because plant genomes contain more duplicated genomic regions than animal genomes. this orthology problem is magnified for non-coding regions, which may evolve and rearrange more rapidly than coding regions. as a result, rates and patterns of sequence evolution among plant non-coding regions have not been characterized in any detail. here we contrast nucleotide substitution rates between orthologous coding and intergenic regions of arabidopsis thaliana and arabidopsis lyrata, two plant species that diverged ~ million years ago. their genomes are largely collinear, but they differ in chromosome number, in dna content, and by several translocations. the two species also differ in population history; a. thaliana is predominantly selfing, while a. lyrata is an obligate outcrosser. these differences should affect differences in patterns and rates of molecular evolution between species, but the expected differences are not readily apparent, perhaps because a. thaliana only recently shifted to a selfing mating system. to study rate variation in intergenic regions, we have generated sequence data in a. lyrata using pcr primers that flank remnants of transposable elements in a. thaliana. we reason that these regions are non-functional by virtue of te insertion, and thus comprise a data set that should be relatively free of selective constraint. we verify their orthology both with comparative maps of the two species and by their molecular features. the intergenic data are contrasted to a second data set consisting of large exonic sequences from a. lyrata and a. thaliana. with these two data sets, we address several questions about arabidopsis nucleotide substitution rates, such as: i) do intergenic sequences evolve at rates similar to synonymous sites in coding data? ii) do any genomic features, like gc content or recombination, correlate with nucleotide substitution rate variation among loci? iii) what can be inferred about the relative contribution of mutation and selection to nucleotide substitution? and iv) do intergenic regions provide any hints to the mechanisms that contribute to genome size differences between a. lyrata and a. thaliana?  isolation and location of orthologous intergenic regions we identified te remnants in the a. thaliana genome and designed pcr primers to flank of these remnants. three primers were designed: two flanked the te remnant, and a third primer was specific to the te. we attempted amplification in a. lyrata with two separate pcr reactions for each of the intergenic regions. the first used the two flanking primers, and the second used one flanking primer with the internal te-specific primer. these two primer sets were also applied to a. thaliana ecotype colombia as a positive control. the intergenic regions were distributed across the a. thaliana genome and ranged in size from 200– bases. of attempted a. lyrata pcr amplifications, 34% were successful. amplifications were considered a success when either of two cases occurred. in case one, the flanking primers amplified a band near the expected size and the internal primer also amplified a band of the expected size. the second successful case was when the flanking primers amplified a band consistent in size with the absence of the te, and the internal primer failed to amplify any band. pcr failure was usually a failure of both primer combinations, but ~10% of pcr "failure" was unexpected pcr results. examples include patterns in which flanking primers amplified a band greatly different than the expected size or in which flanking primers amplified a small band consistent with the loss of the te but amplification with the internal "te" primer was successful. in order to understand sampling biases associated with data generation, we mapped the location of pcr successes and failures along the genome to determine whether successful pcr amplifications were clustered with respect to their relative positions on the a. thaliana genome. we coded successful and unsuccessful pcrs as 1's and 0's, respectively, and created a binary string to represent the linear order of pcr results. we then calculated the variance in the number of zeros bounded by ones or the end of chromosomes. the observed variance in our data was compared to a distribution of variances based on random shufflings of the binary string. the observed variance was much greater than % of simulated variances, indicating that successful pcr amplifications were significantly clustered. this clustering was related, in part, to gene density. we found that attempted amplifications were more likely to be successful in low gene density regions of a. thaliana, when we calculated gene density in a window of mb around the a. thaliana te. consistent with this observation, all attempted amplifications were successful in regions with gene density lower than annotated genes per mb, although there were only amplifications in regions with density this low. this effect seemed to be relatively local, however, because larger mb windows retained the basic trend but the trend was no longer significant. many of the amplicons were small, presumably due to the absence of the te in a. lyrata. we cloned and sequenced a subset of amplicons, yielding alignable sequences. in order to help establish orthology, we examined the position of these amplicons with respect to comparative maps of a. thaliana and a. lyrata. we mapped each amplicon to the a. thaliana genome and found its place between markers that had been mapped on a. lyrata. we considered an amplicon as belonging to a collinear region of the genome if the two markers flanking the marker were also neighboring markers in a. lyrata. by this criterion, we could assign of to unambiguously collinear regions. another were ambiguous either because they had only one flanking marker or because a third marker had been rearranged between the two markers on the a. lyrata map. only one amplicon was clearly in a non-collinear region. although the resolution of comparative a. thaliana-a. lyrata maps is limited and the possibility of paralogy due to segmental duplication cannot be dismissed entirely, most of our intergenic regions were orthologous by the criterion of collinearity. the dna sequence size of our intergenic regions ranged from 218– bp, with an average length of bp. none of the sequences had extensive open reading frames; the longest orf was bp. based on the arabidopsis small rna project targetfinder, there were also no obvious conserved small rnas present in the sequences. of the regions, a. lyrata sequences contained a homolog of the te remnant found in a. thaliana. these sequences were alignable both in the te remnant and in regions flanking the te remnant, providing further molecular evidence for orthology. the remaining a. lyrata sequences did not contain the te remnant, but the sequences were homologous to the regions flanking the a. thaliana te. altogether, 91% of our amplicons were either in a collinear region or contained the te remnant with associated flanking regions, providing strong evidence that our amplicons are orthologs. substitution rates in intergenic vs. synonymous sites given orthologous sequences, we estimated evolutionary distances for intergenic loci. we report distances based on the k2p model, but all results and analyses are qualitatively identical using the general time-reversible model. the average dk2p estimate over all regions was substitutions per site. the dk2p estimates ranged ~5-fold from to substitutions per site. these numbers of course report the distances as if they were sampled from a single distribution, and a priori it might be thought that sequences with tes evolve at different rates than those without. there is, however, no evidence to this effect. similarly, we examined whether evolutionary rates differed by collinearity, reasoning that paralogous sequences could be represented in the group for which collinearity was ambiguous. there was no difference between groups based on collinearity. furthermore, the six amplicons for which there was no evidence for orthology from either te presence or collinearity fell well within the extremes of the dk2p range. because there is no obvious evidence for rate classes based on te presence or collinearity, we treated the intergenic sequences as a single group. in addition to our non-coding data set of sequences, we generated an additional set of large exons. these exons were distributed throughout the genome. forty of these genes are located in non-ambiguously collinear regions between a. thaliana and a. lyrata. of the remaining had ambiguous collinearity and one was in a non-collinear region. we calculated synonymous distances between sequences for all loci. the genes with uncertain collinearity did not differ from the remaining loci either as group, but two of these exons had the highest ds values in our sample, with one value > -fold higher than that of any other exon. we chose to remove these two loci from further analyses because they may represent paralogs, although none of our overall results vary with their inclusion. based on the remaining exons, mean estimates of ds were substitutions per synonymous site, with a range from to. our principle goal is to compare substitution patterns and rates between coding and intergenic regions. loci within the two data types varied ~5-fold in genetic distance. the coefficients of variation, which was for intergenic data and for coding data, indicate that distances vary to similar degrees for the two data types. however, the mean genetic distance differed between the two classes, with synonymous sites evolving more rapidly than non-coding sites. mean distances can be used to estimate nucleotide substitution rates. assuming that the two lineages diverged roughly million years ago, the mean substitution rate for synonymous sites was × 10- substitutions per site per year. in contrast, the mean non-coding rate was only two-thirds of that value, at × 10- substitutions per site per year. genomic correlates with substitution rates our results indicate that mean genetic distances vary between data types and also that genetic distances vary among loci. we also sought to determine whether genetic distances are related to genomic features such as gc content and recombination rate. we expect that such relationships will be clearer in intergenic data, because the patterns are less likely to be complicated by selective forces. recombination may contribute to rate variation among loci by introducing mutations. recombination rates have been estimated for a. thaliana, based on comparisons of physical and genetic distances, but there are no direct estimates of cm/mb recombination rates for a. lyrata. we therefore used a. thaliana recombination rates to investigate relationships between recombination rate and genetic distances. there was no strong correlation between recombination rate and genetic distance whether the data were combined or examined separately as intergenic and genic. the lack of correlation between recombination rates and genetic distances parallels a previous study that detected no correlation between recombination rates and a. thaliana sequence polymorphism. however, nordborg et al. did detect a negative correlation between polymorphism and gene density. given this observation, we too examined the relationship of genetic distance with gene density. gene density is negatively correlated with substitution rate in both datasets, although this relationship is only significant in the intergenic data and not in coding data. these results were generated with gene density measured as the number of genes in a mb window centered around each intergenic locus; the result is similar, but only borderline significant, with mb windows. gc content is another factor known to correlate with substitution rate. the gc content of our intergenic sequences did not differ markedly between a. lyrata and a. thaliana. genic regions also had similar gc content between the two species, but average gc values for coding dna were % higher than intergenic regions. this contrast also held when only synonymous sites in coding regions were considered. the observed transition:transversion ratio was also higher for coding data than for intergenic regions, although statistical support for this trend was borderline. cpg deamination is a common explanation for elevated transition:transversion ratios in high gc regions and could cause differences in rates between genic and intergenic regions. we thus recalculated divergence, treating all cpg dinucleotides as non-variable in both coding and intergenic data sets. although there were more cpg sites in genic regions, the cpg disparity alone did not explain the difference in rate between genic and intergenic regions. genic synonymous divergence rates without cpg sites was × 10- substitutions per site per year, while intergenic divergence rates were × 10- substitutions per site per year. these rates are still significantly different from each other. with combined intergenic and genic data, there was a strong positive correlation between genetic distance and gc content. this correlation also held for synonymous sites alone. to better investigate the relationship between gc and genetic distances, we performed an analysis of covariance. the ancova examined the effect of gc on genetic distance, with sequence type as a factor to determine if gc content alone accounted for the differences in genetic distance between our coding and intergenic data. the ancova model indicated that gc content is a significant predictor for genetic distance, but also that there is an additional effect of sequence type. patterns of indel variation intergenic regions do not have constraints on coding frame, and thus accumulate indels. these indels, in turn, may provide insights into processes that contribute to the × size difference between a. lyrata and a. thaliana genomes. we analyzed indel size distributions in intergenic data. for these analyses, we ignored sequences in which the remnant te was present in a. thaliana but absent in a. lyrata. we also ignored gaps at sequence ends. for the purpose of clarity, we denoted the species with the non-gapped sequence in an alignment as containing an'insertion'. our intergenic data contained distinct a. lyrata insertions totaling bases. a. thaliana, had more distinct insertions but fewer inserted bases. the mean size of insertions was bp in a. lyrata and bp in a. thaliana, and did not differ significantly. in general, a. lyrata insertions tended to be longer and fewer. although the distributions of insertion sizes appear somewhat different, with the a. lyrata distribution having a longer tail of large insertions and a. thaliana having a higher proportion of small insertions, statistical tests comparing the distributions of the samples were not significant. furthermore, a sign test comparing the relative length of a. thaliana and a. lyrata sequences across loci was not significant, and thus there is no evidence that intergenic regions are systematically longer in a. lyrata. discussion we have generated a sample of orthologous intergenic sequences from a. lyrata to compare to a. thaliana. these sequences were designed to include remnants of transposable elements, with the a priori thought that these regions are not under strong selection and therefore that the distribution of evolutionary rates provides insight into neutral genomic processes. nonetheless, our analyses of the genomic location of pcr successes reveal subtle biases in our non-coding data that are important to keep in mind when discussing nucleotide substitution rates. first, the data are not random with respect to genomic location and thus may be biased toward more slowly evolving regions of the genome. second, they tend to originate from regions of low gene density. this does superficially suggest, however, that they are not enriched for functional elements like enhancers or promoters. third, most of the sequences originate from regions of genome collinearity. coupled with molecular features, there is strong evidence that our final data set of non-coding regions represents orthologous dna features. this set of non-coding sequences represents one of the few – and perhaps only – multi-locus, intergenic divergence data in plants. we compared our intergenic data to long exons that were sequenced in a. lyrata. the primary result is that the mean rate of intergenic nucleotide substitution is two-thirds that of the synonymous coding data, with an absolute rate estimated to be × 10- substitutions per site per year. this result holds with alternative nucleotide models, and thus does not appear to be solely an issue of estimation procedures. slower rates in non-coding regions relative to synonymous sites are becoming a surprisingly frequent observation. for example, a recent study of drosophila demonstrated that non-coding dna evolves considerably slower than synonymous sites in terms of both divergence between species and polymorphism within species. by comparing studies, one can also make the case that pseudogenes and introns evolve more slowly than synonymous sites in apes and other mammals. studies of mammalian intergenic regions have also found slower rates than synonymous sites. although most of these studies encompass only a handful of genes, an overall picture of relatively slow non-coding rates is emerging. why are non-coding rates slow relative to synonymous rates? one potential reason discussed above is methodological biases. our observations were consistent across three different alignment methods and were qualitatively similar when we performed the analyses on synonymous and four-fold degenerate sites using the same k2p model as we applied to intergenic data. thus, neither model choice nor alignment appears to contribute substantially to our inferences. nonetheless, one cannot fully discount sampling biases in the noncoding data, particularly because only 34% of our intergenic pcr amplifications were successful. some of this failure is attributable to the fact that some tes, particularly basho elements, do not appear to be shared between the two species and thus could not be amplified in a. lyrata. we attempted to evaluate the magnitude and effect of pcr bias by performing a parallel experiment on coding data. we designed three primers to amplify each of the exons, with the outside primers anchored in flanking intronic or utr regions and an internal primer specific to the exon. to mimic our amplifications of noncoding dna, we applied the same pcr procedures and same criteria as we applied to noncoding data. although of amplifications were successful in a. thaliana, only of were successful in a. lyrata, roughly mimicking the success rate of noncoding amplifications. the group of exons that were successfully amplified by this method did not differ from "unsuccessful" loci in mean genetic distance has concluded that non-coding regions of drosophila are under selective constraint and also subject to bouts of adaptive selection. it seems unlikely that this is a general phenomenon, however, particularly for plants. plant genomes change rapidly in size, in large part due to the accumulation of repetitive dna. it is thus difficult to imagine that this repetitive dna is under strong selective constraint. arabidopsis species may be an exception because they have relatively little repetitive dna for plant genomes. nonetheless, there is little evidence for selection on non-coding regions that flank coding dna. further, the selective interpretation is particularly difficult to argue in this case because the intergenic regions in our study were a priori chosen for their apparent lack of function and tend to originate from genomic regions of low gene density. consistent with the possibility of low functionality, the intergenic data were replete with indel variation. a previous study found that intron lengths from genes differed significantly between these two species, suggesting that the × difference in genome size between a. lyrata and a. thaliana is due, at least in part, to the accumulation of small sequence changes. surprisingly, intergenic indel patterns and lengths do not differ substantially between a. lyrata and a. thaliana, and we find no evidence that differential indel events contribute to differences in genome size. the third explanation for relatively slow rates of evolution in non-coding dna is differential mutation rates. recombination could contribute to this mutational effect, because recombination can cause mutations during the resolution of double-strand breaks. if recombination and mutation are associated processes and most mutations are neutral, then both polymorphism and divergence should be positively correlated with recombination rate. however, we find no correlation between genetic distances and a. thaliana recombination rates. similarly, recombination does not correlate with polymorphism in ~ a. thaliana gene fragments or a. lyrata genes. it is possible, of course, that the lack of correlation between evolutionary rate and recombination is a false negative. a true positive correlation might go undetected if recombination rate estimates are imprecise or if there have been shifts in genomic patterns of recombination between a. thaliana and a. lyrata. to date, there is little evidence for the latter. overall, we have no evidence to suggest that recombination rates contribute to the differences in rates between intergenic and genic dna or to the observed ~5-fold variation in evolutionary rates among loci. in contrast to recombination rates, there is a strong relationship between genetic distance and gc content in our data. substitution rates have long been known to correlate positively with gc content, presumably due to higher mutation rates in high gc regions, in part due to cpg effects. consistent with these previous observations, the synonymous sites of our coding data have a high average gc content relative to intergenic data, faster evolutionary rates, and a higher proportion of transitional changes. interestingly, gc could contribute to some of the non-coding effects observed in drosophila, too. there is a striking difference in gc content between the coding and non-coding regions analyzed by andolfatto; non-coding sequences have a mean gc content of 43%, compared to 58% gc in coding regions. this disparity alone, rather than selective constraint, may account for the slow rate of divergence in drosophila non-coding dna. gc content alone cannot explain, however, the skewed ratios of polymorphism to divergence found by andolfatto, unless there has been a recent shift in gc mutation biases. nonetheless, gc content does not provide a complete explanation for evolutionary rate variation in our data. one of the most interesting aspects of genome-wide analyses in a. thaliana has been the negative correlation between gene density and polymorphism in a. thaliana. this correlation has been interpreted to be the consequence of either selective sweeps or, more likely, background selection. we have uncovered a similar negative correlation between gene density and evolutionary rate for intergenic data, providing a potential neutral explanation for the observation – i.e., that mutation rates are higher in regions of low gene density. note also that, if anything, our intergenic data are biased to low gene density regions of the genome where hitchhiking and background selection should not be particularly strong. the negative correlation between genetic distance and gene density remains borderline significant with partial correlations that also consider gc content. thus, the gene density effect does not appear to be solely an issue of gc, but the causes of this effect are elusive at this point. more importantly, there is a "sequence-type" difference between genic and intergenic substitution rates that is not accounted for by differences in gc content. differential selection of genomic regions is unlikely to explain the observed differences between sequence types, as we would expect the substitution rate effect to be in the opposite direction. it is possible, though, that mutation varies between intergenic and coding regions. for example, transcription-related mutations could increase synonymous substitution rates in coding regions over and above gc effects. this possibility seems plausible because a. thaliana base composition varies as a function of gene expression, suggesting that mutation rates among genes scale with transcription rates. if this hypothesis is correct, then introns should evolve at a rate that is more similar to exons than to intergenic regions. to examine this hypothesis, we aligned the available a. lyrata introns from genbank to their a. thaliana ortholog. the mean distance between intron sequences was substitutions per site. this distance does not differ significantly from our synonymous divergence estimates, but it is significantly higher than the mean genetic distance at our intergenic loci. thus, although the data are limited, intronic sequences are consistent with the hypothesis that transcription-related mutation contributes to differential substitution rates between exons and intergenic regions.  it is clear that gc content is a major determinant of evolutionary rate variation, not only between sequence types but also among loci. gc content may contribute to a. thaliana polymorphism, too, because there is a positive, borderline significant correlation between gc content and snp polymorphism for polymorphic intergenic loci in the nordborg et al. data panel of individuals. on the other hand, gc content does not fully explain intergenic rate variation; variation also correlates with gene density after correcting for gc content. more importantly, there is a detectable effect of sequence type. the explanation that we deem most reasonable for contributing to the sequence-type effect is transcription-related mutation. transcription-related mutation predicts a pattern of higher evolutionary rates in transcribed regions, and this prediction is upheld with both exonic and intronic data. 
as a result, evolutionary rates in intergenic regions may more accurately reflect genomic features, like recombination and gc content, that contribute to nucleotide substitution. our chief observations were that the average rate of nucleotide substitution is slower in intergenic regions than in synonymous sites, that rate variation in both intergenic and coding regions correlate with gc content, that gc content alone is not sufficient to explain differences in rates between intergenic and coding regions, and that rates of evolution in intergenic regions correlate negatively with gene density. the negative correlation between nucleotide substitution rate and gene density provides a potential neutral explanation for a previously documented correlation between gene density and polymorphism levels within a. thaliana. we generated a set of intergenic sequences in arabidopsis lyrata, a close relative of arabidopsis thaliana. our observations indicated that mutation rates vary among genomics regions as a function of base composition, suggesting that previous observations of "selective constraint" on non-coding regions could more accurately be attributed to a gc effect instead of selection. orthology is even more problematic for intergenic regions of plant nuclear genomes, because plant genomes generally harbor a wealth of repetitive dna. the intergenic regions included transposable element remnants and regions flanking the tes. finally, we discuss potential forces that could contribute to rapid synonymous rates, and provide evidence to suggest that transcription-related mutation contributes to rate differences between intergenic and synonymous sites. there has been remarkably little study of nucleotide substitution rate variation among plant nuclear genes, in part because orthology is difficult to establish. in theory orthologous intergenic data is valuable for studying rate variation because nucleotide substitutions in these regions should be under little selective constraint compared to coding regions. we compared substitution rates among the intergenic loci, which exhibit ~5-fold rate variation, and compared intergenic rates to a set of orthologous coding sequences. we verified orthology of these amplified regions both by comparison of existing a. lyrata – a. thaliana genetic maps and by using molecular features.
 there has been remarkably little study of nucleotide substitution rate variation among plant nuclear genes, in part because orthology is difficult to establish. orthology is even more problematic for intergenic regions of plant nuclear genomes, because plant genomes generally harbor a wealth of repetitive dna. in theory orthologous intergenic data is valuable for studying rate variation because nucleotide substitutions in these regions should be under little selective constraint compared to coding regions. as a result, evolutionary rates in intergenic regions may more accurately reflect genomic features, like recombination and gc content, that contribute to nucleotide substitution.  we generated a set of intergenic sequences in arabidopsis lyrata, a close relative of arabidopsis thaliana. the intergenic regions included transposable element remnants and regions flanking the tes. we verified orthology of these amplified regions both by comparison of existing a. lyrata – a. thaliana genetic maps and by using molecular features. we compared substitution rates among the intergenic loci, which exhibit ~5-fold rate variation, and compared intergenic rates to a set of orthologous coding sequences. our chief observations were that the average rate of nucleotide substitution is slower in intergenic regions than in synonymous sites, that rate variation in both intergenic and coding regions correlate with gc content, that gc content alone is not sufficient to explain differences in rates between intergenic and coding regions, and that rates of evolution in intergenic regions correlate negatively with gene density.  our observations indicated that mutation rates vary among genomics regions as a function of base composition, suggesting that previous observations of "selective constraint" on non-coding regions could more accurately be attributed to a gc effect instead of selection. the negative correlation between nucleotide substitution rate and gene density provides a potential neutral explanation for a previously documented correlation between gene density and polymorphism levels within a. thaliana. finally, we discuss potential forces that could contribute to rapid synonymous rates, and provide evidence to suggest that transcription-related mutation contributes to rate differences between intergenic and synonymous sites. 