BACKGROUND
most traits vary quantitatively between individuals and are significantly influenced by genetic variation and its interaction with the environment.
in general, transcript abundance of a gene can be considered as a quantitative trait since it differs between individuals with respect to genetic factors  <cit> .
the utilisation of technologies such as microarray and, lately, next generation sequencing, allowed the simultaneous quantification of thousands of transcripts, shedding light into the genetics behind gene expression variation  <cit> .
comparison in arabidopsis thaliana showed that at least 46% of the genes in the genome are differentially expressed between a pair of accessions  <cit> .
variation in gene expression levels is highly heritable and specific genomic regions can be mapped, underpinning transcript abundance variation  <cit>  and underlying phenotypic diversity, causing for example differential disease susceptibility in humans  <cit>  or variation in metabolite accumulation in plants  <cit> ... previous studies in experimental segregating populations have proved the polygenic nature of gene expression variation among lines and mapped thousands of expression quantitative trait loci  in different model organisms, e.g.
: arabidopsis thaliana  <cit> , saccharomyces cerevisiae  <cit> , and caenorhabditis elegans  <cit> .
no less important is the interaction between genotypes and the environment  <cit> .
for instance, a third of the budding yeast genes showed a significant strain x condition interaction effect, demonstrating the complexity of the control of transcript abundance.
transcript level differences can originate from cis- and/or trans-regulatory changes.
cis- acting regulations represent polymorphisms in physical linkage to the gene, affecting the transcription in an allele-specific manner.
consequently, cis-eqtls are usually detected as local-eqtls, being mapped in the vicinity of the gene; whereas trans-acting regulations act distantly in a non allele-specific way and can be located anywhere in the genome, most likely detected as distant-eqtl  <cit> .
although linkage or association mapping approaches allow the identification of local- and distant-eqtls, cis and trans effects should be directly assessed by allele-specific expression assays  <cit> .
in yeast, most expression differences mapped as major trans-acting loci, where few modifiers regulate the expression of hundreds of genes  <cit> .
contrasting with these results, experiments in c. elegans revealed an opposite pattern, with an over-representation of local-eqtls compared to distant regulators  <cit> .
only few similar studies have addressed this subject in a. thaliana.
using microarray technology in recombinant inbred lines  from a cross between accessions bayreuth  and shahdara , more than  <dig>  eqtls were described using whole rosettes  <cit> .
in contrast, a parallel study using landsberg  and cape verde island  accessions as ril progenitors, detected  <dig> times less eqtls than those mapped in the bay x sha population  <cit> .
discrepancies between the two studies were also evident in the ratio between local- and distant-eqtls, likely due to differences in statistical power and significance thresholds  <cit> .
different growth conditions, stages and the lack of similar analytical strategies in these studies impairs fine comparisons between crosses, stressing the need for a systematic approach to better understand the role of natural variation in expression traits among arabidopsis accessions.
the identification of thousands of local eqtls allows testing for marks of concerted evolution in transcript abundance from genes with related function  <cit> .
several studies have addressed this issue and demonstrated the presence of signatures of selection associated with expression traits in a number of systems  <cit> .
for example, several pathways were identified in a yeast hybrid between two saccharomyces species exhibiting polygenic directional divergence in cis-regulatory variants  <cit> .
similarly, evidence for selection in expression traits involved in processes such as growth, locomotion and memory was detected in two subspecies of mus musculus  <cit> .
thus far, evidence for a directional allelic effect between arabidopsis accessions is scarce.
transcriptome analyses for  <dig> natural accessions determined an overrepresentation of differentially expressed genes within go categories, involving response to the biotic environment, among others  <cit> .
in order to identify genomic regions underlying transcript variation in a. thaliana and signatures of selection between genes with related function, we performed eqtl mapping in two well-established recombinant populations.
the analysis was essentially based on a connected-cross design between three accessions  <cit> : the reference columbia , cape verde island  and burren .
expression level variation was monitored in particular for more than  <dig>  annotated nuclear genes in  <dig> rils and eqtls were mapped using a novel r/eqtl package, accurately designed for such a study.
here, we show that the eqtl landscape is moderately conserved across populations and many of the eqtls identified are cross-specific.
moreover, we show that accounting for confounding sources in one of the crosses increases the sensitivity and the power to detect eqtls, allowing the identification of a potentially coordinated transcriptional evolution in genes with related functions.
RESULTS
eqtl analyses were performed in three displays: two on young rosette transcripts in cvi- <dig> × col- <dig>  and bur- <dig> × col- <dig>  ril sets, and one on developing seed transcripts in bay- <dig> × shahdara  ril set.
we will henceforth focus on and compare the results obtained for rosette samples, while the baysha data is briefly presented as additional file  <dig>
transcript abundance can be treated as a highly heritable phenotype affected by natural genetic diversity  <cit> .
in order to understand the genetics underlying gene expression variation among a. thaliana accessions, we utilised the catma array technology  <cit> , which has the advantage of not being sensitive to punctual polymorphisms  between accessions' transcripts thanks to the length of the probes used .
we measured expression levels for almost  <dig>  traits  in three week-old rosettes of  <dig> rils obtained from a connected-cross set between three parental accessions , including  <dig> rils from cvicol and  <dig> rils from burcol.
initially, in order to characterise expression divergence between the accessions selected in this study, the number of differentially expressed genes between pairs of accessions was estimated .
a total of  <dig>  differentially-expressed genes between col- <dig> and cvi- <dig> and  <dig>  between col- <dig> and bur- <dig>  were detected.
this cut-off value represents, at least,  <dig>  fold expression differences between the two allelic variants.
from this differentially expressed set, none of the parents showed a significant excess of up- or down-regulated genes .
about half of the genes  showing significant expression contrast between col- <dig> and bur- <dig> were also found when comparing the other pair .
interestingly,  <dig> % of this overlap was similarly up- or down-regulated in col- <dig> relative to the other two accessions, suggesting that many of the polymorphisms contributing to these differential expression patterns may arise from col- <dig> or be shared by cvi- <dig> and bur- <dig>
local and distant eqtl distributions are highly divergent between crosses
in order to detect linkage between genomic regions and transcript accumulation variation, we performed eqtl mapping implementing the novel r/eqtl package .
we identified  <dig>  significant eqtls at a fdr of 5%  from the cvicol dataset  and respectively  <dig>  and  <dig>  eqtls from the burcol dataset .
this corresponds to a 1: <dig>  ratio between sets, at most, a contrast essentially relying on eqtls of low significance.
among the set of eqtls in each cross at a fdr of 5%,  <dig>  and  <dig>  were up-regulated in col in the cvicol and burcol populations, respectively.
from the set of differentially expressed genes between the parental accessions,  <dig> % and  <dig> % of them showed a significant linkage to at least one genomic region in the cvicol and burcol crosses respectively.
when comparing the distribution of the phenotypic variance explained per eqtl in both crosses, a highly dissimilar pattern was observed .
in the cvicol set we found a higher average explained variance  with proportionally more eqtls having larger effects .
most of the traits were associated with a unique genomic region , and only a minority exhibited a significant association with two to five eqtls .
when only considering genes located in the nuclear genome,  <dig> %  and 19%  of these had at least one significant eqtl at a fdr of 5%.
altogether, these results demonstrate the complex nature of transcript variation between accessions and suggest a distinct eqtl landscape for each cross.
linkage analysis allows the examination of local and distant eqtls influencing transcript abundance  <cit> .
to detect local associations, we set a 1mb cut-off distance between the eqtl peak and the catma gst physical position.
any eqtl reaching a significance peak within 1mb from the location of the probe it controls was considered as potentially acting in cis and hence classified as a local eqtl.
in both populations we observed that the distributions of local and distant eqtls were not uniform throughout the lod scale .
distant eqtls were enriched at lower lods, compared to local eqtls, which were over-represented at high significance values.
hence, 50%  of all eqtls detected at 5% fdr mapped locally in the cvicol  population, while this number would rise to  <dig> %  if using a much more conservative fdr of  <dig> % .
moreover, compared to distant eqtls, local eqtls explained a significantly higher fraction of the phenotypic variance per trait .
these results can lead to contrasting conclusions as we vary threshold values .
the higher phenotypic variance explained by local eqtls underlies the fact that most of the phenotypic variation in expression levels is due to local associations, and trans-acting factors would have essentially minor effects on single transcripts, but overall a widespread effect over many genes.
in order to classify cross-specific and overlapping eqtls, we estimated the number of local eqtls shared between the burcol and cvicol populations.
we found that  <dig> transcripts mapped to an overlapping local eqtl interval in the populations, corresponding to  <dig> % of all local eqtls in the cvicol set .
the direction of the allelic effect for these potentially shared eqtls was consistent in  <dig> % of the cases.
moreover, within this subset, we observed a strong correlation between populations for the strength of the additive effect associated with the eqtls  and for the explained phenotypic variance .
we also compared the distribution of local eqtls along the genome between crosses.
for this, we divided the genome into 1mb windows  and estimated the fraction of transcripts with local eqtls on every bin per chromosome.
a spearman test detected a significant correlation between both distributions , suggesting major structural genomic aspects.
we determined whether dense eqtl regions were enriched for snp polymorphisms or were gene-dense regions.
a genome-wide regression analysis did not detect a significant correlation between the number of snps and eqtls within intervals for cvicol, however a marginally significant association was found for burcol .
in contrast, an expected significant correlation was found with gene density in both cases .
furthermore, we compared the top  <dig> bins  in both sets and detected  <dig> overlapping intervals with high eqtl density between crosses .
interestingly, the  <dig> intervals were not among the top polymorphic or gene-dense and therefore the overlapping cannot be solely attributed to these factors.
distant eqtl hotspots are cross-specific
we considered the set of distant eqtls mapped at 5% fdr to test for the presence of hotspots along the genome.
similarly to the local eqtl analysis, we divided the genome into 1mb windows and identified regions containing a higher number of eqtls than expected by chance, according to a permutation test .
in the cvicol population, a total of  <dig> bins contained more than the  <dig> distant eqtl peaks expected by chance, and up to  <dig> peaks were found within the strongest hotspot .
these regions included  <dig> % of all the distant eqtls detected.
moreover, on two occasions, two or more adjacent bins were above the threshold, likely highlighting the wide distribution of eqtls around a major hotspot due to a lack of precision when mapping small-effect loci.
for example, the hotspot localised on chromosome  <dig> included four bins and  <dig> % of all the distant-eqtls, the strongest one in the genome.
within this interval,  <dig> eqtls are located at the most dense bin  with an average lod of  <dig> .
a different depiction was obtained in the burcol population: we detected  <dig> bins containing more eqtls than expected by chance .
in this population the hotspots encompassed the majority  of all the distant-eqtls.
moreover, within the three major hotspots containing several significant adjacent bins on chromosome  <dig>   <dig> and  <dig> , we found a significant enrichment for distally-regulated genes of  <dig> go categories in the branch 'biological process'.
these categories included response to chitin  on chromosome  <dig> and cell-adhesion  on chromosome  <dig> .
we found little overlap between hotspots in the two populations .
from the  <dig> hotspots detected in total, only three were collinear between sets, suggesting little conservation in the major determinants of distant eqtls across the genome.
none of the hotspots co-localised with any of the major-effect developmental qtls known to segregate in a. thaliana , except in burcol where a marginally significant hotspot was found on chromosome  <dig> , which contains flc.
although  <dig> traits with a distant eqtl mapping to this region were observed, only  <dig> have been shown to interact with flc  <cit> , suggesting potentially a minor contribution from this locus.
the gene-poor centromeres rich in transposons and pseudogenes were not detected in any of the hotspot intervals.
the lack of a significant count of eqtls extended for at least 1mb from the centromere until the nearest hotspot, except in the burcol set where we detected distant eqtls-rich regions immediately downstream and upstream the centromeres on chromosome  <dig> and  <dig>  respectively .
controlling for confounding factors increases eqtl detection in the cvicol dataset
variation in transcript abundance can also be influenced by non-genetic factors, such as experimental noise and hidden factors  <cit> : greenhouse localisation, temperature variation, sampling time and other undetected experimental perturbations can affect expression levels among samples, altering the subsequent analysis.
a recent association mapping study in a. thaliana already demonstrated the greater resolution obtained when considering confounding sources between genotypes  compared with standard approaches  <cit> .
hence, it is interesting to try and control these variables in order to improve the eqtl mapping.
for this purpose, we utilised the probabilistic approach vbqtl together with r/eqtl to account for hidden confounding factors  <cit> .
we tested this approach in the cvicol dataset for  <dig>   <dig>   <dig>   <dig> and  <dig> hidden factors, finding the lowest number of eqtls for no-hidden factors and the highest for  <dig> factors .
we chose one known factor  and ten unknown factors to model the expression trait, since the lowest number of distant associations were explained away, in contrast to  <dig> and  <dig> factors.
contrarily, using the same number of factors  in the burcol dataset, we found that many of the distant associations were lost .
using  <dig> factors we mapped  <dig>  eqtls at 5% fdr  in the burcol population, while only  <dig> % of the distant associations were retained by both the standard and the vbqtl methods.
many of the distant associations exhibiting high significance values were present in both datasets; however, several distant associations with marginally significant lods were lost .
the majority of the local associations detected by the standard approach were also mapped by vbqtl .
even more, we observed a 60% increase in the number of local eqtls detected by vbqtl.
nevertheless, the large number of genetic associations explained away in the burcol set hindered the subsequent data analysis, since many distant eqtls affecting transcript levels should be considered  and not treated as confounding sources  <cit> .
we focused the utilisation of this approach on the cvicol population.
linkage analysis detected a total of  <dig>  eqtls , which is almost twice as many eqtls as with the standard method, demonstrating the sensitivity of the approach.
interestingly, the fraction of local and distant eqtls remained almost unchanged.
at the less conservative threshold  we mapped  <dig> % of the eqtls in the vicinity of the transcript, slightly less than from the analysis of the raw expression values.
moreover, the contribution to the total number of local and distant eqtls was increased by  <dig>  and  <dig>  eqtls respectively.
many of the local associations mapped by the standard approach were also mapped by vbqtl  and, contrarily to the burcol data set,  <dig> % of the distant associations were also retained .
approximately one third of the peaks that had remained only marginally significant  in the standard analysis, were now detected as significant at 5% fdr after vbqtl analysis, likely due to the negative effect of confounding sources on power in the original approach.
evidence for directional allelic effect in col- <dig> versus cvi- <dig> accessions
the independence of eqtls is crucial for the identification of selection on transcript abundance in a group of genes with related functions  <cit> .
in order to detect evidence for a directional allelic effect in the cvicol population, we used the larger set of local eqtls identified by the vbqtl approach.
we considered functional categories with an overrepresentation of eqtls and a significant bias in the number of genes with either col or cvi alleles upregulated.
to test for categories with more eqtls than expected by chance , the gene functional classification assigned by the gene ontology consortium was used  <cit> .
in particular, we focused on the branch 'biological process', since go categories here contained could give insight into the divergent adaptation events occurred in each accession.
we identified  <dig> significant overrepresented categories : "response to stress", "response to biotic and abiotic stimulus" and "transport".
from the  <dig> go sets within these categories, we only analysed those containing more than  <dig> genes .
subsequently, we estimated the number of eqtls within each set either up or down regulated in any accession and test for deviation from the genome-wide  <dig> : <dig> col:cvi eqtl ratio using a hypergeometric test  <cit> .
we found  <dig> sets showing a significant skew , all of them part of the response to stress or biotic and abiotic stimulus, including: defense, hypoxia and plant-type hypersensitive response.
interestingly, in all cases most of the members were upregulated in the col accession .
these results demonstrate the presence of a systematic directional allelic effect in sets of genes having related functions and validate the utilisation of expression traits to identify potential intraspecific diversification events.
number of eqtls per categories among go 'biological process'.
categories with a significant excess  of eqtls are indicated with 's*'.
the number of expected eqtls was estimated from a whole-genome scan considering the fraction of genes in each functional category
go terms showing a significant directional allelic effect.
the number of genes per term, eqtls observed and number of up-regulating allele from each parent is shown, as well as the p value associated with the allelic bias
discussion
we have performed quantitative genetic analyses for genome-wide expression traits in a connected-cross set between three accessions in the model organism a. thaliana.
utilising the r/eqtl package in both crosses, we observed a  <dig> % difference in the number of significant eqtls between sets, with the bur- <dig> × col- <dig> cross being the one with the highest count .
regardless of this difference, the conservative figures of eqtls obtained in both sets resemble the one previously described in rils between accessions ler- <dig> and cvi- <dig>  <cit> , but they significantly differ from a previous study involving bay- <dig> and shahdara accessions, where more than  <dig>  eqtls were detected  <cit> .
in our study, the use of a common approach for both sets allowed us to perform straightforward comparisons between them, drawing conclusions not necessarily affected by differences in statistical power  <cit> .
we observed that  <dig> % of all local eqtls detected were shared between crosses, likely due to polymorphisms solely present in col- <dig> or shared by cvi- <dig> and bur- <dig> .
an argument to support this hypothesis is the high level of shared eqtls with an additive effect in the same direction  and of similar strength .
however, the low level of overlapping eqtls demonstrates the greater number of cross-specific eqtls, possibly due to the great complexity of expression traits  <cit>  and their high potential for evolution  <cit> .
our results are consistent with a recent survey in several arabidopsis accessions that identified high allelic heterogeneity within local regulatory eqtls  <cit> , suggesting the existence of many private polymorphisms between all three accessions extending the eqtl landscape.
although cis-acting regulation among local eqtls remains to be confirmed, the co-localisation of the eqtl and the controlled gene provides a robust probabilistic way of classifying eqtls into their mode of action.
we observed that, albeit the two crosses shared a common accession, the eqtl distribution among these types is very dissimilar.
the burcol set shows three times more distant than local eqtls.
at the same fdr, the cvicol set has a very different balance, with an equivalent number of both types of eqtls  and a lower number of hotspots  encompassing  <dig> % of the trans eqtls.
nevertheless, the distribution pattern of local eqtls on each chromosome was conserved between crosses, likely due to a strong correlation with gene density and chromosome structural features.
also, we found that local eqtls were overrepresented at higher significance in all crosses , explaining a greater phenotypic variance per trait compared to distant eqtls.
our result agrees with another study in mouse, where many highly significant local-eqtls and moderately significant trans-acting associations were observed  <cit> .
furthermore, the burcol recombinant population exhibited a higher number of distant eqtls, concentrating 75% of them in  <dig> bins classified as hotspots .
one of the most interesting intervals was the hotspot at bin # <dig>  which contained the highest number of genes whose function is related to the regulation of transcription , including the strong candidate gigantea  <cit> , a gene known to be involved in diverse developmental processes  <cit> .
it is difficult to guess whether a hotspot could be the expression of a locus actually controlling the transcription of many genes, or that of a major  player that has indirect effects on many genes' expression .
the presence of pleiotropic hotspots affecting the expression of many transcripts has already been described in yeast  <cit> , where no more than  <dig> trans-eqtls explained transcript variation for  <dig>  traits.
similarly, many trans-hotspots were described in previous a. thaliana studies  <cit> .
in agreement with these studies, the hotspots detected here explained approximately 10% of the individual traits' variation, demonstrating their milder effect on transcript abundance compared to local eqtls  <cit> .
interestingly, only ~7% of the hotspots overlapped between crosses, suggesting either alternative master regulators or little conservation in complex regulatory networks  <cit> .
these conserved regions could also suggest the presence of polymorphisms in col- <dig> as major contributors for the large number of trans eqtls within these shared regions.
to increase the power in detecting eqtls and identify milder genetic variations, we have utilised the bayesian framework vbqtl, a model designed to dissect gene expression variation in order to account for confounding factors  <cit> .
the implementation of this model in the cvicol set allowed us to detect twice as many eqtls compared to the standard approach , without affecting the ratio of local versus distant eqtls.
subsequently, we used this extended set and focused the posterior analyses on local eqtls since they likely represent independent alleles, only affecting a single transcript and explaining a greater phenotypic variance, in contrast to distant eqtls that may affect many genes and show minor phenotypic effects.
the use of this new set of local eqtls allowed us to look for potential traces of selection where col alleles were overexpressed in a set of genes with related function.
assuming neutrality, no over-representation of either up-regulated or down-regulated alleles from the same accession is expected within a cluster of genes, unless a directional allelic drift has occurred, which may represent selection  <cit> .
we detected an overrepresentation of eqtls in several go functional categories and a significant skew within go terms .
the accumulation of local regulatory polymorphisms up-regulating col alleles in all cases, suggests dissimilar patterns of responses to the environment and an advantage in a particular niche  <cit>  or, else, a specific lack of cost in loosing this response.
for example, one of the clusters of genes globally up-regulated in col with respect to cvi relates to hypoxia response, which is a major determinant of submergence tolerance.
it was indeed recently shown that col and cvi contrasts for this trait, with cvi being one of the least submergence-tolerant accessions  <cit> .
interestingly, the entire set of significant categories here described are related to response to stress or biotic and abiotic response, highlighting its potential link with environmental and niche-adaptation changes.
CONCLUSIONS
our results demonstrate the importance of extending the number of mapping populations in order to understand the eqtl landscape within a species.
we found that the distribution of local and distant eqtls is moderately conserved between recombinant populations, demonstrating the complex inheritance of transcript abundance.
moreover, the identification of candidate pathways for signatures of selection is a significant step towards understanding accession diversification and their adaptation to the environment.
further studies using novel technologies, such as next generation sequencing , association mapping, or the combination of both, will help elucidating and validating transcript abundance divergence with a greater resolution within a. thaliana species.

 most traits vary quantitatively between individuals and are significantly influenced by genetic variation and its interaction with the environment. in general, transcript abundance of a gene can be considered as a quantitative trait since it differs between individuals with respect to genetic factors. the utilisation of technologies such as microarray and, lately, next generation sequencing, allowed the simultaneous quantification of thousands of transcripts, shedding light into the genetics behind gene expression variation. comparison in arabidopsis thaliana showed that at least 46% of the genes in the genome are differentially expressed between a pair of accessions. variation in gene expression levels is highly heritable and specific genomic regions can be mapped, underpinning transcript abundance variation and underlying phenotypic diversity, causing for example differential disease susceptibility in humans or variation in metabolite accumulation in plants... previous studies in experimental segregating populations have proved the polygenic nature of gene expression variation among lines and mapped thousands of expression quantitative trait loci in different model organisms, e.g.: arabidopsis thaliana, saccharomyces cerevisiae, and caenorhabditis elegans. no less important is the interaction between genotypes and the environment. for instance, a third of the budding yeast genes showed a significant strain x condition interaction effect, demonstrating the complexity of the control of transcript abundance. transcript level differences can originate from cis- and/or trans-regulatory changes. cis- acting regulations represent polymorphisms in physical linkage to the gene, affecting the transcription in an allele-specific manner. consequently, cis-eqtls are usually detected as local-eqtls, being mapped in the vicinity of the gene; whereas trans-acting regulations act distantly in a non allele-specific way and can be located anywhere in the genome, most likely detected as distant-eqtl. although linkage or association mapping approaches allow the identification of local- and distant-eqtls, cis and trans effects should be directly assessed by allele-specific expression assays. in yeast, most expression differences mapped as major trans-acting loci, where few modifiers regulate the expression of hundreds of genes. contrasting with these results, experiments in c. elegans revealed an opposite pattern, with an over-representation of local-eqtls compared to distant regulators. only few similar studies have addressed this subject in a. thaliana. using microarray technology in recombinant inbred lines from a cross between accessions bayreuth and shahdara, more than eqtls were described using whole rosettes. in contrast, a parallel study using landsberg and cape verde island accessions as ril progenitors, detected times less eqtls than those mapped in the bay x sha population. discrepancies between the two studies were also evident in the ratio between local- and distant-eqtls, likely due to differences in statistical power and significance thresholds. different growth conditions, stages and the lack of similar analytical strategies in these studies impairs fine comparisons between crosses, stressing the need for a systematic approach to better understand the role of natural variation in expression traits among arabidopsis accessions. the identification of thousands of local eqtls allows testing for marks of concerted evolution in transcript abundance from genes with related function. several studies have addressed this issue and demonstrated the presence of signatures of selection associated with expression traits in a number of systems. for example, several pathways were identified in a yeast hybrid between two saccharomyces species exhibiting polygenic directional divergence in cis-regulatory variants. similarly, evidence for selection in expression traits involved in processes such as growth, locomotion and memory was detected in two subspecies of mus musculus. thus far, evidence for a directional allelic effect between arabidopsis accessions is scarce. transcriptome analyses for natural accessions determined an overrepresentation of differentially expressed genes within go categories, involving response to the biotic environment, among others. in order to identify genomic regions underlying transcript variation in a. thaliana and signatures of selection between genes with related function, we performed eqtl mapping in two well-established recombinant populations. the analysis was essentially based on a connected-cross design between three accessions: the reference columbia, cape verde island and burren. expression level variation was monitored in particular for more than annotated nuclear genes in rils and eqtls were mapped using a novel r/eqtl package, accurately designed for such a study. here, we show that the eqtl landscape is moderately conserved across populations and many of the eqtls identified are cross-specific. moreover, we show that accounting for confounding sources in one of the crosses increases the sensitivity and the power to detect eqtls, allowing the identification of a potentially coordinated transcriptional evolution in genes with related functions.  eqtl analyses were performed in three displays: two on young rosette transcripts in cvi- × col- and bur- × col- ril sets, and one on developing seed transcripts in bay- × shahdara ril set. we will henceforth focus on and compare the results obtained for rosette samples, while the baysha data is briefly presented as additional file transcript abundance can be treated as a highly heritable phenotype affected by natural genetic diversity. in order to understand the genetics underlying gene expression variation among a. thaliana accessions, we utilised the catma array technology, which has the advantage of not being sensitive to punctual polymorphisms between accessions' transcripts thanks to the length of the probes used. we measured expression levels for almost traits in three week-old rosettes of rils obtained from a connected-cross set between three parental accessions, including rils from cvicol and rils from burcol. initially, in order to characterise expression divergence between the accessions selected in this study, the number of differentially expressed genes between pairs of accessions was estimated. a total of differentially-expressed genes between col- and cvi- and between col- and bur- were detected. this cut-off value represents, at least, fold expression differences between the two allelic variants. from this differentially expressed set, none of the parents showed a significant excess of up- or down-regulated genes. about half of the genes showing significant expression contrast between col- and bur- were also found when comparing the other pair. interestingly, % of this overlap was similarly up- or down-regulated in col- relative to the other two accessions, suggesting that many of the polymorphisms contributing to these differential expression patterns may arise from col- or be shared by cvi- and bur- local and distant eqtl distributions are highly divergent between crosses in order to detect linkage between genomic regions and transcript accumulation variation, we performed eqtl mapping implementing the novel r/eqtl package. we identified significant eqtls at a fdr of 5% from the cvicol dataset and respectively and eqtls from the burcol dataset. this corresponds to a 1: ratio between sets, at most, a contrast essentially relying on eqtls of low significance. among the set of eqtls in each cross at a fdr of 5%, and were up-regulated in col in the cvicol and burcol populations, respectively. from the set of differentially expressed genes between the parental accessions, % and % of them showed a significant linkage to at least one genomic region in the cvicol and burcol crosses respectively. when comparing the distribution of the phenotypic variance explained per eqtl in both crosses, a highly dissimilar pattern was observed. in the cvicol set we found a higher average explained variance with proportionally more eqtls having larger effects. most of the traits were associated with a unique genomic region, and only a minority exhibited a significant association with two to five eqtls. when only considering genes located in the nuclear genome, % and 19% of these had at least one significant eqtl at a fdr of 5%. altogether, these results demonstrate the complex nature of transcript variation between accessions and suggest a distinct eqtl landscape for each cross. linkage analysis allows the examination of local and distant eqtls influencing transcript abundance. to detect local associations, we set a 1mb cut-off distance between the eqtl peak and the catma gst physical position. any eqtl reaching a significance peak within 1mb from the location of the probe it controls was considered as potentially acting in cis and hence classified as a local eqtl. in both populations we observed that the distributions of local and distant eqtls were not uniform throughout the lod scale. distant eqtls were enriched at lower lods, compared to local eqtls, which were over-represented at high significance values. hence, 50% of all eqtls detected at 5% fdr mapped locally in the cvicol population, while this number would rise to % if using a much more conservative fdr of %. moreover, compared to distant eqtls, local eqtls explained a significantly higher fraction of the phenotypic variance per trait. these results can lead to contrasting conclusions as we vary threshold values. the higher phenotypic variance explained by local eqtls underlies the fact that most of the phenotypic variation in expression levels is due to local associations, and trans-acting factors would have essentially minor effects on single transcripts, but overall a widespread effect over many genes. in order to classify cross-specific and overlapping eqtls, we estimated the number of local eqtls shared between the burcol and cvicol populations. we found that transcripts mapped to an overlapping local eqtl interval in the populations, corresponding to % of all local eqtls in the cvicol set. the direction of the allelic effect for these potentially shared eqtls was consistent in % of the cases. moreover, within this subset, we observed a strong correlation between populations for the strength of the additive effect associated with the eqtls and for the explained phenotypic variance. we also compared the distribution of local eqtls along the genome between crosses. for this, we divided the genome into 1mb windows and estimated the fraction of transcripts with local eqtls on every bin per chromosome. a spearman test detected a significant correlation between both distributions, suggesting major structural genomic aspects. we determined whether dense eqtl regions were enriched for snp polymorphisms or were gene-dense regions. a genome-wide regression analysis did not detect a significant correlation between the number of snps and eqtls within intervals for cvicol, however a marginally significant association was found for burcol. in contrast, an expected significant correlation was found with gene density in both cases. furthermore, we compared the top bins in both sets and detected overlapping intervals with high eqtl density between crosses. interestingly, the intervals were not among the top polymorphic or gene-dense and therefore the overlapping cannot be solely attributed to these factors. distant eqtl hotspots are cross-specific we considered the set of distant eqtls mapped at 5% fdr to test for the presence of hotspots along the genome. similarly to the local eqtl analysis, we divided the genome into 1mb windows and identified regions containing a higher number of eqtls than expected by chance, according to a permutation test. in the cvicol population, a total of bins contained more than the distant eqtl peaks expected by chance, and up to peaks were found within the strongest hotspot. these regions included % of all the distant eqtls detected. moreover, on two occasions, two or more adjacent bins were above the threshold, likely highlighting the wide distribution of eqtls around a major hotspot due to a lack of precision when mapping small-effect loci. for example, the hotspot localised on chromosome included four bins and % of all the distant-eqtls, the strongest one in the genome. within this interval, eqtls are located at the most dense bin with an average lod of. a different depiction was obtained in the burcol population: we detected bins containing more eqtls than expected by chance. in this population the hotspots encompassed the majority of all the distant-eqtls. moreover, within the three major hotspots containing several significant adjacent bins on chromosome and, we found a significant enrichment for distally-regulated genes of go categories in the branch'biological process'. these categories included response to chitin on chromosome and cell-adhesion on chromosome. we found little overlap between hotspots in the two populations. from the hotspots detected in total, only three were collinear between sets, suggesting little conservation in the major determinants of distant eqtls across the genome. none of the hotspots co-localised with any of the major-effect developmental qtls known to segregate in a. thaliana, except in burcol where a marginally significant hotspot was found on chromosome, which contains flc. although traits with a distant eqtl mapping to this region were observed, only have been shown to interact with flc, suggesting potentially a minor contribution from this locus. the gene-poor centromeres rich in transposons and pseudogenes were not detected in any of the hotspot intervals. the lack of a significant count of eqtls extended for at least 1mb from the centromere until the nearest hotspot, except in the burcol set where we detected distant eqtls-rich regions immediately downstream and upstream the centromeres on chromosome and respectively. controlling for confounding factors increases eqtl detection in the cvicol dataset variation in transcript abundance can also be influenced by non-genetic factors, such as experimental noise and hidden factors: greenhouse localisation, temperature variation, sampling time and other undetected experimental perturbations can affect expression levels among samples, altering the subsequent analysis. a recent association mapping study in a. thaliana already demonstrated the greater resolution obtained when considering confounding sources between genotypes compared with standard approaches. hence, it is interesting to try and control these variables in order to improve the eqtl mapping. for this purpose, we utilised the probabilistic approach vbqtl together with r/eqtl to account for hidden confounding factors. we tested this approach in the cvicol dataset for and hidden factors, finding the lowest number of eqtls for no-hidden factors and the highest for factors. we chose one known factor and ten unknown factors to model the expression trait, since the lowest number of distant associations were explained away, in contrast to and factors. contrarily, using the same number of factors in the burcol dataset, we found that many of the distant associations were lost. using factors we mapped eqtls at 5% fdr in the burcol population, while only % of the distant associations were retained by both the standard and the vbqtl methods. many of the distant associations exhibiting high significance values were present in both datasets; however, several distant associations with marginally significant lods were lost. the majority of the local associations detected by the standard approach were also mapped by vbqtl. even more, we observed a 60% increase in the number of local eqtls detected by vbqtl. nevertheless, the large number of genetic associations explained away in the burcol set hindered the subsequent data analysis, since many distant eqtls affecting transcript levels should be considered and not treated as confounding sources. we focused the utilisation of this approach on the cvicol population. linkage analysis detected a total of eqtls, which is almost twice as many eqtls as with the standard method, demonstrating the sensitivity of the approach. interestingly, the fraction of local and distant eqtls remained almost unchanged. at the less conservative threshold we mapped % of the eqtls in the vicinity of the transcript, slightly less than from the analysis of the raw expression values. moreover, the contribution to the total number of local and distant eqtls was increased by and eqtls respectively. many of the local associations mapped by the standard approach were also mapped by vbqtl and, contrarily to the burcol data set, % of the distant associations were also retained. approximately one third of the peaks that had remained only marginally significant in the standard analysis, were now detected as significant at 5% fdr after vbqtl analysis, likely due to the negative effect of confounding sources on power in the original approach. evidence for directional allelic effect in col- versus cvi- accessions the independence of eqtls is crucial for the identification of selection on transcript abundance in a group of genes with related functions. in order to detect evidence for a directional allelic effect in the cvicol population, we used the larger set of local eqtls identified by the vbqtl approach. we considered functional categories with an overrepresentation of eqtls and a significant bias in the number of genes with either col or cvi alleles upregulated. to test for categories with more eqtls than expected by chance, the gene functional classification assigned by the gene ontology consortium was used. in particular, we focused on the branch'biological process', since go categories here contained could give insight into the divergent adaptation events occurred in each accession. we identified significant overrepresented categories: "response to stress", "response to biotic and abiotic stimulus" and "transport". from the go sets within these categories, we only analysed those containing more than genes. subsequently, we estimated the number of eqtls within each set either up or down regulated in any accession and test for deviation from the genome-wide: col:cvi eqtl ratio using a hypergeometric test. we found sets showing a significant skew, all of them part of the response to stress or biotic and abiotic stimulus, including: defense, hypoxia and plant-type hypersensitive response. interestingly, in all cases most of the members were upregulated in the col accession. these results demonstrate the presence of a systematic directional allelic effect in sets of genes having related functions and validate the utilisation of expression traits to identify potential intraspecific diversification events. number of eqtls per categories among go'biological process'. categories with a significant excess of eqtls are indicated with's*'. the number of expected eqtls was estimated from a whole-genome scan considering the fraction of genes in each functional category go terms showing a significant directional allelic effect. the number of genes per term, eqtls observed and number of up-regulating allele from each parent is shown, as well as the p value associated with the allelic bias discussion we have performed quantitative genetic analyses for genome-wide expression traits in a connected-cross set between three accessions in the model organism a. thaliana. utilising the r/eqtl package in both crosses, we observed a % difference in the number of significant eqtls between sets, with the bur- × col- cross being the one with the highest count. regardless of this difference, the conservative figures of eqtls obtained in both sets resemble the one previously described in rils between accessions ler- and cvi-, but they significantly differ from a previous study involving bay- and shahdara accessions, where more than eqtls were detected. in our study, the use of a common approach for both sets allowed us to perform straightforward comparisons between them, drawing conclusions not necessarily affected by differences in statistical power. we observed that % of all local eqtls detected were shared between crosses, likely due to polymorphisms solely present in col- or shared by cvi- and bur-. an argument to support this hypothesis is the high level of shared eqtls with an additive effect in the same direction and of similar strength. however, the low level of overlapping eqtls demonstrates the greater number of cross-specific eqtls, possibly due to the great complexity of expression traits and their high potential for evolution. our results are consistent with a recent survey in several arabidopsis accessions that identified high allelic heterogeneity within local regulatory eqtls, suggesting the existence of many private polymorphisms between all three accessions extending the eqtl landscape. although cis-acting regulation among local eqtls remains to be confirmed, the co-localisation of the eqtl and the controlled gene provides a robust probabilistic way of classifying eqtls into their mode of action. we observed that, albeit the two crosses shared a common accession, the eqtl distribution among these types is very dissimilar. the burcol set shows three times more distant than local eqtls. at the same fdr, the cvicol set has a very different balance, with an equivalent number of both types of eqtls and a lower number of hotspots encompassing % of the trans eqtls. nevertheless, the distribution pattern of local eqtls on each chromosome was conserved between crosses, likely due to a strong correlation with gene density and chromosome structural features. also, we found that local eqtls were overrepresented at higher significance in all crosses, explaining a greater phenotypic variance per trait compared to distant eqtls. our result agrees with another study in mouse, where many highly significant local-eqtls and moderately significant trans-acting associations were observed. furthermore, the burcol recombinant population exhibited a higher number of distant eqtls, concentrating 75% of them in bins classified as hotspots. one of the most interesting intervals was the hotspot at bin # which contained the highest number of genes whose function is related to the regulation of transcription, including the strong candidate gigantea, a gene known to be involved in diverse developmental processes. it is difficult to guess whether a hotspot could be the expression of a locus actually controlling the transcription of many genes, or that of a major player that has indirect effects on many genes' expression. the presence of pleiotropic hotspots affecting the expression of many transcripts has already been described in yeast, where no more than trans-eqtls explained transcript variation for traits. similarly, many trans-hotspots were described in previous a. thaliana studies. in agreement with these studies, the hotspots detected here explained approximately 10% of the individual traits' variation, demonstrating their milder effect on transcript abundance compared to local eqtls. interestingly, only ~7% of the hotspots overlapped between crosses, suggesting either alternative master regulators or little conservation in complex regulatory networks. these conserved regions could also suggest the presence of polymorphisms in col- as major contributors for the large number of trans eqtls within these shared regions. to increase the power in detecting eqtls and identify milder genetic variations, we have utilised the bayesian framework vbqtl, a model designed to dissect gene expression variation in order to account for confounding factors. the implementation of this model in the cvicol set allowed us to detect twice as many eqtls compared to the standard approach, without affecting the ratio of local versus distant eqtls. subsequently, we used this extended set and focused the posterior analyses on local eqtls since they likely represent independent alleles, only affecting a single transcript and explaining a greater phenotypic variance, in contrast to distant eqtls that may affect many genes and show minor phenotypic effects. the use of this new set of local eqtls allowed us to look for potential traces of selection where col alleles were overexpressed in a set of genes with related function. assuming neutrality, no over-representation of either up-regulated or down-regulated alleles from the same accession is expected within a cluster of genes, unless a directional allelic drift has occurred, which may represent selection. we detected an overrepresentation of eqtls in several go functional categories and a significant skew within go terms. the accumulation of local regulatory polymorphisms up-regulating col alleles in all cases, suggests dissimilar patterns of responses to the environment and an advantage in a particular niche or, else, a specific lack of cost in loosing this response. for example, one of the clusters of genes globally up-regulated in col with respect to cvi relates to hypoxia response, which is a major determinant of submergence tolerance. it was indeed recently shown that col and cvi contrasts for this trait, with cvi being one of the least submergence-tolerant accessions. interestingly, the entire set of significant categories here described are related to response to stress or biotic and abiotic response, highlighting its potential link with environmental and niche-adaptation changes.  our results demonstrate the importance of extending the number of mapping populations in order to understand the eqtl landscape within a species. we found that the distribution of local and distant eqtls is moderately conserved between recombinant populations, demonstrating the complex inheritance of transcript abundance. moreover, the identification of candidate pathways for signatures of selection is a significant step towards understanding accession diversification and their adaptation to the environment. further studies using novel technologies, such as next generation sequencing, association mapping, or the combination of both, will help elucidating and validating transcript abundance divergence with a greater resolution within a. thaliana species. 
here, we present the transcriptome architecture in leaf tissue of two ril sets obtained from a connected-cross design involving commonly used accessions. expression traits can vary quantitatively between individuals and have a complex inheritance. identification of the genetics underlying transcript variation can help in the understanding of phenotypic variation due to genetic factors regulating transcript abundance and shed light into divergence patterns. our results demonstrate the extensive diversity and moderately conserved eqtl landscape between crosses and validate the utilisation of expression traits to explore for candidates behind phenotypic variation among accessions. moreover, using local eqtls obtained from this approach, we detected evidence for a directional allelic effect in genes with related function, where significantly more eqtls than expected by chance were up-regulated from one of the accessions. the utilisation of the novel r/eqtl package allowed us to map and eqtls in the cvi- × col- and bur- × col- recombinant populations respectively. so far, only a limited number of studies have addressed this subject in arabidopsis, with contrasting results due to dissimilar statistical power. primary experimental data, analysis parameters, eqtl results and visualisation of lod score curves presented here are stored and accessible through the qtlstore service database http://qtlstore.versailles.inra.fr/. in agreement with previous studies, we observed a larger phenotypic variance explained by eqtls in linkage with the controlled gene, compared to distant loci.
 expression traits can vary quantitatively between individuals and have a complex inheritance. identification of the genetics underlying transcript variation can help in the understanding of phenotypic variation due to genetic factors regulating transcript abundance and shed light into divergence patterns. so far, only a limited number of studies have addressed this subject in arabidopsis, with contrasting results due to dissimilar statistical power. here, we present the transcriptome architecture in leaf tissue of two ril sets obtained from a connected-cross design involving commonly used accessions. we also present the transcriptome architecture observed in developing seeds of a third independent cross.  the utilisation of the novel r/eqtl package allowed us to map and eqtls in the cvi- × col- and bur- × col- recombinant populations respectively. in agreement with previous studies, we observed a larger phenotypic variance explained by eqtls in linkage with the controlled gene, compared to distant loci. distant eqtls hotspots were essentially not conserved between crosses, but instead, cross-specific. accounting for confounding factors using a probabilistic approach increased the mapping resolution and the number of significant associations. moreover, using local eqtls obtained from this approach, we detected evidence for a directional allelic effect in genes with related function, where significantly more eqtls than expected by chance were up-regulated from one of the accessions. primary experimental data, analysis parameters, eqtl results and visualisation of lod score curves presented here are stored and accessible through the qtlstore service database http://qtlstore.versailles.inra.fr/.  our results demonstrate the extensive diversity and moderately conserved eqtl landscape between crosses and validate the utilisation of expression traits to explore for candidates behind phenotypic variation among accessions. furthermore, this stresses the need for a wider spectrum of diversity to fully understand expression trait variation within a species. eqtlnatural variationselectionrils 