Cubillos et al. BMC Genomics 2012, 13:117 http://www.biomedcentral.com/1471-2164/13/117 hypoxia response, which is a major determinant of sub- mergence tolerance. It was indeed recently shown that Col and Cvi contrasts for this trait, with Cvi being one of the least submergence-tolerant accessions [33]. Inter- estingly, 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 under- stand the eQTL landscape within a species. We found that the distribution of local and distant eQTLs is mod- erately 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 (RNA-seq), association mapping, or the combination of both, will help elucidating and validating transcript abundance divergence with a greater resolu- tion within A. thaliana species. Methods Plant material, sample preparation and microarray hybridization We used 158 RILs from the core-population of the Cvi-0 × Col-0 set and 156 RILs from the core-population of the Bur-0 × Col-0 set [21], along with their respective par- ents. Plants were cultivated in a greenhouse in typical long day conditions (16h photoperiod) at 20°C and whole plants were collected above the roots 20 days after sow- ing, which corresponds on average to growth stage ‘1.08’ from Boyes et al. [34]. Total RNA was extracted using RNeasy Plant mini Kit (Qiagen kit #74904) according to the supplier’s instructions. For each biological sample (RIL), total RNAs were obtained by pooling RNAs from three randomized plants within a single experiment. Transcript abundance estimates were carried out using the bicolour CATMA microarrays version 5, containing 34,529 spots (including multiple internal controls) allow- ing, among others, to uniquely interrogate 26,166 nuclear genes (predicted from TAIR and EUGENE algo- rithms) in A. thaliana [22,35]. For each comparison of two RILs (random pair design), one technical replicate with fluorochrome reversal was performed (i.e. four hybridisations -2 arrays- per comparison of two RILs). Labelling of cRNAs with Cy3-dUTP or Cy5-dUTP (Per- kin-Elmer-NEN Life Science Products), hybridisation to the slides, and scanning procedures were performed as previously described [36]. Page 9 of 12 Microarray data analysis For each array, the raw data comprised the logarithm of the median feature pixel intensity at wavelengths 635 nm (red) and 532 nm (green) and no background was subtracted. An array-by-array normalisation was per- formed to remove systematic biases. First, spots consid- ered as badly formed features were excluded. Then a global intensity-dependent normalisation using the loess procedure was performed to correct the dye bias [37]. Finally, for each block, the log-ratio median calculated over the values for the entire block was subtracted from each individual log-ratio value to correct print tip effects. Differential analysis was based on the log ratios averaged on the dye-swap: the technical replicates were averaged to get one log-ratio per biological replicate and these values were used to perform a paired t-test. A trimmed variance is calculated from spots which do not display extreme variance [38]. The raw P-values were adjusted by the Bonferroni method, which controls the Family Wise Error Rate in order to keep a strong con- trol of the false positives in a multiple-comparison con- text. We considered probes as being differentially expressed after Bonferroni correction using a P < 0.05. For the eQTL analysis, normalised intensity per probe × RIL has been calculated from the normalised log-ratio by sharing the correction value between both samples [39]. To be specific, after the normalization of each array, we get a log-ratio denoted M and a mean inten- sity denoted A for each probe. First, we average M and A across the two arrays taking into account the dye switching and then calculate the normalised intensity of the two co-hybridized RILs which are equal to (2A-M)/2 and (2A + M)/2. The intensity (2A-M)/2 corresponds to the RIL labelled in green on the first array of the dye- swap. Then a between-array normalisation was per- formed to rescale the mean intensity of each slide at an arbitrary value equal to 8.5. This second round of nor- malisation makes probe × RIL intensity comparable and gives us a single expression level for each probe × RIL that can be used further for the eQTL mapping strategy. Data Deposition Microarray data was deposited at Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/, accession no. GSE28791 and at CATdb (http://urgv.evry.inra.fr/ CATdb/; Project: GNP07_RILKIT) according to the ‘Minimum Information About a Microarray Experiment’ standards. R/eqtl package and eQTL mapping eQTL mapping was performed utilising the normalised microarray data using our R package R/eqtl http://cran. r-project.org/web/packages/eqtl/, which exploits func- tions from R/qtl [40]. R/eqtl is designed specifically for