BACKGROUND
genetic and epigenetic variation form the basis for local adaptation and speciation processes, and are becoming increasingly accessible through advances in genomic and bioinformatic tools.
the advent of microarray and ultra-high throughput sequencing  technologies have thus brought about a renewed interest in evolutionary questions, with a prospect for gaining novel insights at the whole-genome level.
these opportunities have spurred genome-wide surveys of single nucleotide polymorphisms   <cit>  and methylation polymorphisms in many organisms including plants, for example in multiple accessions of the genetic model organism arabidopsis thaliana and in closely related species .
in attempts to identify causative genetic changes in plant adaptations, classical linkage analysis and genome-wide association studies  have successfully mapped traits governing the performance under local environmental conditions to snps at specific loci  <cit> .
structural variation in the form of gene copy number variation  polymorphism is an influential component of natural genetic diversity that markedly contributes to phenotypic variation  <cit> .
however, cnv has been addressed in noticeably fewer studies because of technical difficulties in its comprehensive and reliable assessment  <cit> .
short cnvs consisting of insertions or deletions below  <dig> kb in size can be readily detected based on uhts technologies.
however, the identification of cnvs comprising from  <dig> kb up to one or multiple genes has generally remained challenging.
genome-wide analyses in human and other mammalian model organisms revealed cnvs to be much more abundant than previously known, e.g.
affecting 10% of the mouse genome and 12% of the human genome .
cnvs have been implicated in human disease etiology, and evidence for adaptive cnvs is also emerging  <cit> .
in comparison to mammalian genomes, gene duplications and deletions especially from whole genome duplications appear to be even more abundant in plant genomes  <cit> .
single-gene and segmental duplications as well as whole-genome duplications have been hypothesized to propel adaptive evolution and speciation.
in plants, this view is supported by recent reports on cultivar-specific boron tolerance in barley  <cit> , aluminium tolerance in maize  <cit>  and species-wide heavy metal tolerance in the wild plant arabidopsis halleri  <cit> , all supporting the role of gene copy number expansion in plant adaptation to abiotic stress.
population genomic data, for example from arabidopsis thaliana and zea mays, have identified an unexpectedly high abundance of cnvs  <cit> , generating interest in the contribution of structural mutations to genome plasticity.
ten percent of maize genes were found to exhibit copy number polymorphisms, and an experimental evolution study in a. thaliana reported de novo structural mutations resulting in  <dig> copy number variant genes after only  <dig> generations  <cit> .
although between-species genome comparisons have remained difficult to date, the few existing studies have supported the hypothesis that gene copy number expansions, and especially those involving tandem duplications  <cit> , might underlie plant adaptations to environmental stress  <cit> .
given that novel functions are much more likely to be generated by adaptive specialization of one of several pre-existing copies of a duplicated gene than by an entirely novel gene  <cit> , such comparative studies are key to understanding the patterns of genomic polymorphisms associated with adaptation and speciation.
the availability of a well-annotated genome sequence and a wealth of knowledge on gene functions for arabidopsis thaliana, as well as for several closely related species that have diverged over short evolutionary timespans, render arabidopsis a suitable model genus to study adaptation and speciation processes  <cit> .
one of its species is arabidopsis halleri — a wild outcrossing, zn and cd hyperaccumulating and hypertolerant species that is naturally found on both highly metal-contaminated and non-contaminated soils   <cit> .
its genome is expected to be about 25% larger than that of its non-accumulating, non-tolerant sister species a. lyrata, and about 65% larger than the size of the genome of a. thaliana, from which both species diverged between  <dig>   <cit>  and  <dig> million years ago  <cit> .
the striking phenotypic contrast between a. halleri and the two other closely related species, despite a high sequence similarity within coding regions  <cit> , provides an exceptional opportunity to elucidate molecular evolutionary patterns reflecting the influence of natural soil characteristics on adaptation and speciation.
fig.
<dig> comparison of the metal hyperaccumulator species arabidopsis halleri to the closely related non-hyperaccumulator species a. lyrata and a. thaliana.
a representative photograph is shown for each species, together with the estimated evolutionary distances separating them, given as the divergence times from a common ancestor  <cit> .
listed below are key phenotypic and genomic characteristics.
mya, million years ago
previous cross-species transcriptomics studies identified a number of differentially expressed candidate genes for the metal hyperaccumulation/hypertolerance trait of a. halleri .
among these, heavy metal atpase  <dig> , which encodes a metal pump that acts as an exporter of zn2+ and cd2+ from specific cell types, was shown to be necessary both for the hyperaccumulation of zn and for the full extent of zn and cd hypertolerance  <cit> .
strongly increased hma <dig> transcript levels in a. halleri were attributed to a lineage-specific tandem triplication combined with cis-regulatory mutations  <cit> .
an analysis of sequence polymorphism in the genomic region of hma <dig> gene copy number expansion demonstrated strong positive selection, as well as selection for enhanced hma <dig> gene product dosage  <cit> .
another candidate gene, nicotianamine synthase  <dig> , encodes an enzyme that catalyses the biosynthesis of the low-molecular-weight metal chelator nicotianamine from s-adenosyl methionine, and was shown to contribute to zn hyperaccumulation  <cit> .
in addition to hma <dig>  several other transition metal homeostasis candidate genes of a. halleri were demonstrated to be copy number expanded through the dna gel  blot technique  <cit> .
the objective of the work presented here was to identify genes exhibiting copy number expansion in a. halleri at a genome-wide scale, in relation to the known species-specific extreme traits.
we conducted a survey of gene copy number divergence  in a. halleri relative to a. thaliana by employing array-comparative genomic hybridization  in a cross-species manner using the ath <dig> microarray designed for a. thaliana.
in order to test whether the identified cnds are species- and thus possibly trait-specific, our analysis included a. lyrata as a third species, which is a non-tolerant non-hyperaccumulator like a. thaliana, but shares with a. halleri an equal phylogenetic distance from a. thaliana.
we devised a novel routine for evaluating cross-species array-cgh data, which is based on the quantification and subsequent global correction of the effects of sequence divergence on hybridization signal intensities.
our procedure operates without loss of probe information, which is crucial for retaining statistical power for cnv estimation further downstream.
our predictions of genic cnds were validated against a small set of genes with known copy number in a. halleri  <cit>  and against a set of genes predicted to be copy number expanded or reduced according to the a. lyrata reference genome sequence  <cit> .
gene copy number expansions in a. halleri, but not in a. lyrata, were found to be significantly enriched for metal homeostasis functions.
conversely, biotic stress functions were significantly enriched among genes exhibiting copy number reduction in a. halleri, but not in a. lyrata.
these results suggest that between-species divergence in gene copy numbers reflects adaptive evolution of metal hyperaccumulation, a species-specific trait of a. halleri that has been proposed to provide an elemental defence against biotic stress  <cit> .
RESULTS
metal hyperaccumulation and hypertolerance in a. halleri have previously been attributed to the constitutively enhanced expression of a number of metal homeostasis genes, several of which were additionally shown to be expanded in genomic copy number through dna gel blots  <cit> , bac sequencing  <cit>  or other methods  <cit> .
here, the technique of cross-species array-cgh was employed for a genome-wide assessment of between-species divergence in gene copy number.
fragmented and labelled nuclear genomic dna samples from a. thaliana and from the two closely related heterologous species a. halleri and a. lyrata, were hybridized in duplicate to affymetrix ath <dig> genechip®; microarrays   <cit> .
the challenge in cross-species hybridizations of genomic dna is that the target gene sequences of the heterologous species inevitably contain mismatches relative to the probe sequences of the reference species on the microarray, thus reducing the efficiency of hybridization and resulting in lowered signals  <cit> .
we employed a novel approach to correct this bias by using a signal adjustment strategy, which - unlike previous methods  <cit>  — accounts for sequence mismatches through a global adjustment of cross-species hybridization signal intensities.
in brief, we implemented a two-step normalization scheme .
the first step was a conventional within-species normalization, which was applied to raw signal intensities from each pair of two replicate microarray hybridizations of the same target species.
the second step was an adjustment of normalized signal intensities through the calculation and application of a species-specific global scaling factor for compensating the effects of sequence divergence from a. thaliana.
fig.
<dig> overview of the data analysis workflow.
flowchart summarizing our two-step normalization approach for the processing of cross-species genomic hybridization data, consisting of within-species normalization and global scaling of signals through species-specific signal correction factors, followed by the final prediction of copy number divergent genes.
grey arrows and backgrounds mark the auxiliary steps taken for the determination of species-specific global scaling factors with the aid of reference gene datasets
to implement our strategy, we began by establishing a representative subset of probes for which curated sequence data was available from a. halleri, termed reference dataset .
a similar reference dataset was also generated at random from the available reference genome sequence of a. lyrata.
for each heterologous target species, the signal correction factor was calculated from the statistical distribution of the occurrence of mismatches and the ensuing effect on hybridization signal intensity as measured in the respective reference dataset.
subsequently, the normalized and corrected cross-species hybridization data were analysed for differential signals between species in order to identify putative copy number divergent genes.
finally, a comparison between copy number alterations in a. lyrata and a. halleri enabled us to identify species-specific copy number alterations.
consequences of inter-species sequence divergence for mismatch occurrence between probe sequences and heterologous target sequences
for the adjustment of microarray signals in cross-species array-cgh, we generated one reference dataset of representative, curated sequence data for each of the two heterologous target species.
the a. halleri reference dataset comprised  <dig> genes, yielding  <dig> matching probe sequences on the microarray .
because of the lack of a reference genome, these data corresponded to previously obtained sequences from a. halleri ssp.
halleri   <cit> .
the a. lyrata reference dataset comprised  <dig> genes with  <dig> matching probe sequences on the microarray, obtained from the published reference genome  <cit>  .
the number and positions of mismatches between each heterologous target sequence and the corresponding microarray probe sequence was determined .
for a. halleri and a. lyrata, respectively,  <dig> and 35% of all probe sequences were fully conserved across species,  <dig> and 29% contained only a single mismatch with respect to the  <dig> nucleotide-long probe sequence, and  <dig> and 30% of sequences contained between  <dig> and  <dig> mismatches compared to the corresponding probe sequence .
fig.
<dig> frequency distribution of mismatch occurrence between microarray probe sequences and heterologous target gene sequences.
shown is the percentage of a. thaliana probes on the ath <dig> array that display no mismatches up to  <dig> mismatches  when hybridized to non-a.
thaliana genomic dna from either a. halleri  or a. lyrata .
the expected frequency distribution  is shown by the grey line, and was calculated based on the average coding sequence identity  within transcribed regions between a. halleri and a. thaliana
we computed the expected distribution of the total number of mismatches of a given nucleotide sequence with respect to the probe sequence on the microarray .
the expected distribution  of mismatches in cross-species hybridization of a. halleri gdna closely matched our observations for the reference dataset .
both the observed and expected probabilities of the occurrence of more than  <dig> mismatches by comparison to the probe sequence were negligible .
finally, the observed mismatch distributions for a. halleri and a. lyrata were highly similar to each other , thus confirming similar levels of sequence divergence from a. thaliana  and allowing us to use the expected distribution of mismatches calculated for a. halleri also in a. lyrata hybridizations.
both the a. halleri and a. lyrata lineages are thought to have diverged from the common ancestor with a. thaliana at the same point in the past .
our observations support the theory that a correlation exists between the levels of sequence divergence and actual phylogenetic distances between species, as was estimated, for example, based on cross-species array-cgh data  <cit> .
quantification of effects of sequence mismatches on signal intensities in cross-species microarray hybridization
sequence mismatches are known to be the single most confounding factor biasing the signals of cross-species array hybridizations.
a previous study has estimated sequence mismatches to account for at least 40% of the average noise in microarray hybridization  <cit> , and several studies confirm mismatches as the primary cause of failure of conventional normalization techniques in cross-species microarray data analysis  <cit> .
as a result of sequence divergence, sequence mismatches are expected to reduce the hybridization efficiency of genomic dna from a. halleri and a. lyrata to the a. thaliana probe sequences on the ath <dig> microarray, resulting in lowered overall hybridization signal intensity.
after background correction of raw data , we examined the influence of the total number and positions of mismatches on the normalized hybridization signal intensities using the probe signal intensities from our a. halleri and a. lyrata reference datasets.
as expected, hybridization signal intensity decreased with increasing number of mismatches in a probe.
the largest decrease in signal intensity by  <dig> and 40% in a. halleri and a. lyrata, respectively, was observed for a single mismatch .
additional mismatches had only small effects, with a total of four mismatches resulting in a further reduction of probe signal intensity by  <dig> and 7% in a. halleri and a. lyrata, respectively.
there were only minor differences between a. halleri and a. lyrata in the dependence of signal intensity on the number of mismatches.
note that the published a. lyrata reference genome, which was used here to determine the number of mismatches between a. lyrata reference dataset target sequences and the corresponding ath <dig> probe sequences, was from the north american subspecies lyrata, whereas our hybridization experiments were conducted with the european ssp.
petrea  <cit> .
a stark reduction in diversity has been reported in a. lyrata ssp.
lyrata by comparison to a. lyrata ssp.
petrea, and several studies  report differentiation between the two sub-species that have been isolated from each other for between  <dig> and  <dig> thousand years  <cit> .
fig.
<dig> dependence of hybridization signal intensity on number and position of mismatches with respect to the probe sequence on the ath <dig> array.
a values are arithmetic means  of background-corrected raw probe signal intensity ratios for non-a.
thaliana gdna relative to a. thaliana gdna hybridizations, shown as a function of the total number of mismatches of the heterologous target sequence compared to the corresponding a. thaliana 25-mer probe sequence.
b independence of hybridization signal intensity from the position of a single mismatch with respect to the probe sequence.
values are arithmetic means  of background-corrected raw probe signal intensity ratios for non-a.
thaliana gdna relative to a. thaliana gdna hybridizations, shown as a function of mismatch position in the heterologous target sequence compared to the corresponding a. thaliana probe sequence.
black circles represent the representative a. halleri reference dataset; white diamonds represent the representative a. lyrata reference dataset
surprisingly, we observed a noisy profile of signal intensity over different positions of a single mismatch along the probe sequence instead .
the expected sharp drop in signal intensity when a single mismatch is positioned in the centre  of a probe sequence, as proposed by affymetrix for so-called mismatch  probes  <cit> , was not detected here.
this finding is in agreement with a number of previous studies  <cit> , which have pointed out that experimental data do not conform to this postulate and that, in fact, for some probes signal intensity was even found to be higher for mm probes than for perfectly matching  probes  <cit> .
consequently, position-based effects on hybridization signal intensity are hard to construct, and accordingly, the most popular normalization methods no longer take the information from mm probes into account.
therefore, for our between-species normalization strategy, we did not consider the influence of sequence mismatch position.
we estimated the incremental signal correction factor s
k for a probe with k mismatches as the average of the ratio of normalized hybridization signal intensity of a. thaliana to the respective signal intensity of the heterologous species.
for each probe containing  <dig> to  <dig> mismatches, incremental signal correction factors were weighted by their probability of occurrence , followed by the calculation of the arithmetic mean to yield species-specific global scaling factors .
these global signal correction factors of  <dig>  for hybridizations of a. halleri gdna and  <dig>  for hybridizations of a. lyrata gdna were employed to scale the hybridization signal intensities of the respective cross-species microarray hybridizations.
cross-species normalization and validation of copy number divergent genes
the median raw signal intensities for the heterologous species a. halleri and a. lyrata were lower than those for the ath <dig> target model species a. thaliana, namely by  <dig> and 36%, respectively .
after applying conventional within-species vsn normalizations, median normalized signal intensities were more uniform across replicates within each species .
however, the differences between species were large, with median signal intensities for a. halleri and a. lyrata which were  <dig> and 49% lower, respectively, than for a. thaliana.
upon the subsequent between-species scaling of vsn-normalized signal intensities from the a. halleri and a. lyrata hybridizations employing species-specific global scaling factors , median signal intensities became more similar across species and remained 17% lower for a. halleri and 11% lower for a. lyrata, respectively, than for a. thaliana .
following normalization and scaling of hybridization signals,  <dig>  and  <dig> copy number expanded genes  were identified in a. halleri and a. lyrata , respectively .
furthermore,  <dig> genes predicted to be copy number reduced  were identified in a. halleri and  <dig> in a. lyrata .
overall,  <dig> cnes and  <dig> cnrs are shared between a. halleri and a. lyrata compared to the reference species a. thaliana .
thus, based on a.
thaliana as a reference, an about 3-fold larger number of copy number divergent genes was detected in a. halleri than in a. lyrata, whereas nucleotide sequence divergence from a. thaliana was similar in both heterologous species within transcribed regions .
the observed difference between a. halleri and a. lyrata was not merely a spurious result caused by a higher level of polymorphism between the two replicate a. halleri samples than between the a. lyrata replicates.
this was confirmed by performing all data processing steps with two additional pairs of a. halleri replicates, each consisting of either one of the two single a. halleri hybridizations and one a. halleri replicate generated in silico by simulating between-replicate variation as observed in a. lyrata, respectively .
consequently, our results suggest that the rate of either acquisition or maintenance of gene copy number changes in the genome can differ between closely-related lineages or species.
this is in stark contrast to the general stability of base substitution rates normalized to genome size and generation  <cit> .
fig.
<dig> distribution of signal intensities before and after normalization and scaling.
boxplot of  background-corrected raw hybridization signal intensities,  normalized signal intensities after vsn normalization of the replicate arrays of each species, respectively, and  signal intensities after the application of species-specific global scaling factors to the normalized data.
boxes show median, and upper and lower quartiles, of l
o
g
<dig> probe signal intensities for each gdna hybridization.
upper and lower horizontal bars mark all values lying within  <dig>  times the inter-quartile range.
replicate hybridizations are denoted  <dig> and  <dig> and grouped by species
to evaluate the reliability of our predicted cnds, we compared our results to genes of known copy number status.
for a. halleri, we used a set of  <dig> genes .
the evaluation of microarray-based predictions of cross-species gene cnd against known copy number status  indicated  <dig> % specificity,  <dig> % precision and  <dig> % sensitivity of our cross-species array-cgh based cnd estimation.
for a. lyrata, a complete reference genome sequence is available  <cit> .
this provides an opportunity for more extensive data validation by comparing our predictions of gene cnds with predictions based on the reference genome sequence.
orthology predictions retrieved from ensembl plants indicated that, relative to all a. thaliana genes represented on the ath <dig> genechip,  <dig>  orthologous genes of a. lyrata are copy number expanded, whereas  <dig>  genes are reduced in copy number or deleted .
from this set of genes, we further chose the subset of highly conserved multi-copy genes of ≥ 95% average sequence identity between transcribed regions of all paralogs with respect to their a. thaliana ortholog.
our final reference dataset of highly conserved copy number expanded genes in a. lyrata, as predicted based on ensembl plants, contained  <dig> genes.
by comparison, our method yielded  <dig> % specificity,  <dig> % precision and  <dig> % sensitivity .
these scores suggest that our procedure is conservative, making predictions with high reliability, and sacrificing sensitivity for a higher specificity.
a direct comparative assessment of the performance of our array-cgh based method for a. halleri and a. lyrata is not possible because of the differing qualities of experimental evidence underlying the two datasets available for validation.
darby et al.
2011d
machado et al.
2010e
a
a. lyrata genes sharing ≥ 95% sequence identity with their closest a. thaliana homologue  <cit>  are termed highly conserved
bheaders of half-columns refer to total number of cnes predicted by ensembl plants  alone or additionally by a. lyrata genome analysis , respectively, as given in parentheses here.
shown are commonalities with these two groups of genes  or data referring to these two groups of genes
ctrue positive is a cne detected based on array-cgh that was previously predicted to be a cne by ensembl plants  <cit>  alone, or additionally by a. lyrata genome  <cit>  analysis
d <cit>
e <cit>
we compared the performance of our array-cgh based approach with that of the two previous studies that also aimed at estimating cnd using the array-cgh technique  <cit> .
we reproduced the normalization and scaling strategies of machado and renn  and darby et al.
as described  <cit> , with few small modifications necessary to apply these methods to our array-cgh platform .
the method of darby et al.
resulted in the prediction of a  <dig> -fold elevated number of gene copy number expansions.
out of the two previously published methods, maximum sensitivity, specificity and precision of the detection of copy number expansion among highly conserved genes were  <dig> ,  <dig>  and  <dig> %, respectively, all inferior to our method .
even for the genes that are not highly conserved but predicted to be copy number expanded concordantly by both ensembl plants and the a. lyrata genome project, our method reports higher sensitivity, specificity and precision –  <dig>   <dig>  and  <dig> % respectively than previous studies  <cit>  –  <dig> ,  <dig>  and  <dig> % .
specificity and precision of our method were also superior concerning copy number reductions or gene deletions .
functional analysis of copy number divergent genes of a. halleri
after identifying the sets of genes exhibiting copy number divergence by comparison to the reference species a. thaliana in either of the two heterologous species according to array-cgh, we evaluated these for any enrichment of functional categories using the mapman ontology  <cit> .
copy number expanded genes of a. halleri showed a statistically significant enrichment for transition metal homeostasis-related gene functions , and for mitochondrial electron transport/atp synthesis-related functions , relative to all genes represented on the array .
among copy number reduced genes of a. halleri, there was a significant enrichment for biotic stress-related functions .
in contrast, none of these functional categories was detected to be significantly enriched among genes divergent in copy number in a. lyrata by comparison to a. thaliana.
in the light of the two extreme traits specific to a. halleri — namely zn/cd hyperaccumulation and associated hypertolerance — the genome-wide overrepresentation of metal homeostasis genes among copy number expanded genes is remarkable .
indeed, a high occurrence of copy number expansion was reported among arabidopsis halleri metal homeostasis candidate genes identified through the presence of elevated transcript levels in a. halleri compared to a. thaliana  <cit> .
moreover, copy number expansion is known to contribute to high hma <dig> transcript levels in a. halleri, which in turn are necessary for both metal hyperaccumulation and the full extent of metal hypertolerance  <cit> .
hma <dig> gene copy number expansion is not limited to a. halleri, but also found in the zn/cd hyperaccumulator species noccaea caerulescens, similarly associated with strongly elevated transcript levels  <cit> .
a. halleri mtp <dig> is another copy number-expanded candidate gene, for which several lines of evidence suggest an involvement in zn hypertolerance  <cit> .
it was not known to date whether these findings on individual candidate genes pertain at the genome-wide level, but this is now supported by array-cgh data presented here.
our data additionally confirm the previous finding of zip <dig> copy number expansion  <cit> .
in contrast, our array-cgh analysis did not detect hma <dig> as copy number expanded in a. halleri, although this is well established.
one of the transition metal homeostasis candidate genes newly identified to be copy number expanded in a. halleri is nas <dig>  which was demonstrated to be highly expressed in roots of a. halleri  <cit>  and to contribute to zn hyperaccumulation  <cit> .
array-cgh also predicts ahhma <dig> to be copy number expanded.
this candidate gene was reported as highly expressed in a. halleri, and an ahhma <dig> cdna confers zn and cd tolerance upon heterologous expression in yeast  <cit> .
finally, atpcr <dig> and atpcr <dig> have been implicated in the export of zn and cd, respectively, from cells  <cit>  and appear to be copy number expanded in a. halleri.
indeed, transcript levels were found to be higher in a. halleri than in a. thaliana in a previous microarray hybridization study  <cit> , but this was not explicitly reported because of an ambiguous assignment of the probe set to this pair of highly similar a. thaliana genes.
sams <dig> encodes an enzyme that catalyses the biosynthesis of the substrate for nicotianamine synthase and was previously identified to be more highly expressed in a. halleri than in a. thaliana; and it is now identified as copy number expanded.
expanding from these findings, the enrichment of mitochondrial electron transport functions among the genes copy number expanded in a. halleri was surprising.
a majority of  <dig> copy number expansions, half of which are shared with a. lyrata, affect complex i acting as nadh ubiquinone oxidoreductase in oxidative phosphorylation .
a more parsimonious explanation for the common expansions could be that these genes have been deleted in a. thaliana relative to the common ancestor.
the two genes privately copy number expanded in a. halleri have predicted roles in atp synthesis and cytochrome c biogenesis, respectively, and one is annotated as a mitochondrial dicarboxylate carrier.
fig.
<dig> significantly enriched functional classes among  copy number expanded and  copy number reduced genes of a. halleri.
bars/brackets represent the proportion/number of genes in each functional class among all copy number expanded or reduced genes of a. halleri  or a. lyrata , and among all nuclear genes represented on the array  for comparison.
asterisks mark statistically significant enrichment by comparison to the set of nuclear genes represented on the ath <dig> microarray .
the analysis was based on mapman  <cit>  functional categories .
shown are only top-level categories of the ontological hierarchy
a. halleri vs. a. thaliana
a. lyrata vs. a. thaliana
l
o
g
<dig> fcb
l
o
g
<dig> fcb
mfm1
pcr1; pcr2
<dig>
<dig>
chx25
261448_atc
itp2
<dig>
<dig>
chx2
267304_atd
zip6
fer4
266718_atc,d
mtp1
mt2a
fer3
cch
sams2
as8
253658_atc,d
hma3
hcc2
mt2b
248048_atc
nas2
<dig>
<dig>
ysl4
zip3
ptac17
ti1
250944_atd
hipp6
hipp43
nas-like
fp3
copt <dig>  emb1513
259871_atc
256930_atd
oasa2
253413_atc
zip9
pdf <dig>
sams3/ mat4*
mtp7
260913_atd
sams1
chx16
257365_x_atc,d
pdf <dig> b
249255_atc
chx18
263838_atc,d
mat3**
hipp27
cax11
252694_atc
pdf <dig>
zip11
mt3
258987_atd
hcc1
a. halleri vs. a. thaliana
a. lyrata vs. a. thaliana
l
o
g
<dig> fcb
l
o
g
<dig> fcb
<dig>
<dig>
b12-1
<dig>
<dig>
dic2
<dig>
<dig>
<dig>
<dig>
b <dig>
<dig>
<dig>
atpq
249627_atd
emb1467
ccme
hcc2
phb6; phb1
b14
a. halleri vs. a. thaliana
a. lyrata vs. a. thaliana
l
o
g
<dig> fcb
l
o
g
<dig> fcb
– <dig>
– <dig>
– <dig>
– <dig>
rlp40; rlp42
ttr1; wrky16
rlp34; rlp53
transition metal homeostasis.
mitochondrial electron transport/atp synthesis.
biotic stress mapman functional classes
ashort names of genes previously reported to be copy-number expanded in a. halleri based on other experimental approaches are given in bold type
b
l
o
g
<dig> fold-change  is underlined for genes that are concordantly copy number expanded/reduced in both a. halleri and a. lyrata
chigher transcript abundance in roots of a. halleri vs. a. thaliana
dhigher transcript abundance in shoots of a. halleri vs. a. thaliana
∞synthesis of cytochrome c oxidase;
$zinc-regulated transporter, iron-regulated transporter-like; *at3g <dig> was annotated as sams <dig> in ref.
talke et al.
2006; **at2g <dig> was annotated as sams <dig> in talke et al.
2006; &nadh-ubiquinone oxidoreductase, mitochondrial respiratory chain complex i;#coiled coil  or toll-interleukin receptor -nucleotide binding site -leucine rich repeat ;+nucleotide-binding  motif-containing domain shared by human apaf- <dig>  certain plant r gene products, and caenorhabditis elegans ced- <dig>
among the genes predicted to be copy number reduced according to our array-cgh results, we observed an overrepresentation of biotic defence functions .
this is an interesting observation with respect to the proposed ecological role of metal hyperaccumulation in plants .
almost all biotic stress-related genes predicted to be copy number reduced in a. halleri encode members of large protein families typically involved in plant innate immunity and designated as disease resistance  genes, such as the predominating tir-nbs-lrr receptor kinases .
infecting pathogens generate characteristic molecular patterns that can be specifically recognized by cognate r gene products, which subsequently trigger a localized cell death response that is essential for plant disease resistance.
the enrichment of r gene-related biotic stress functions among genes reduced in copy number in the hyperaccumulator a. halleri supports the elemental defence hypothesis as well as the trade-off hypothesis for the evolutionary role of elemental defence.
accordingly, elemental defences through metal hyperaccumulation allow for the loss of r genes, thus alleviating the fitness costs associated with r gene expression.
at lower levels of the ontological hierarchy, post-translational modification functions, and genes encoding glycine-rich and crinkly-like proteins were over-represented among the copy number expanded genes of a. halleri, whereas dc <dig> domain-containing, as well as protein degradation/e3-scf-f-box and cysteine protease functions were significantly enriched among the genes predicted to be copy number reduced.
furthermore, genes predicted to be copy number reduced in both a. halleri and a. lyrata according to our array-cgh, showed a significant enrichment of the mapman functional category or bin “dna.synthesis/chromatin structure.retrotransposon/transposase” .
contrary to this observation, it is known that transposable element gene families have undergone an expansion in the a. lyrata genome relative to a. thaliana  <cit> .
contributions of segmental duplications or deletions to copy number divergence
next we tested whether large segmental duplications or deletions make major contributions to cnd from a. thaliana in a. halleri or a. lyrata.
we detected only few large segmental duplications  and deletions  comprising more than  <dig> kb , by using a sliding window approach and the set parameters .
a total of six predicted segmental deletion regions, located in positions corresponding to a. thaliana chromosomes  <dig> to  <dig>  were shared between a. halleri and a. lyrata.
most detected segmental deletions were located in or around regions corresponding to a. thaliana centromers.
these segmental deletions were composed almost entirely of transposable element genes  <cit>  and are thus likely to include regions of low sequence conservation rather than segmental deletions .
the sequences of centromeric regions are known to be highly divergent in both nucleotide composition and length even between closely related species  <cit> .
none of the mapman functional classes of top level ontological hierarchy that were significantly enriched among copy number divergent genes  was found similarly enriched among the genes contained in segmental duplications or deletions.
overall, cnds, and especially cnes, appeared to be dispersed evenly in the genomes of both a. halleri and a. lyrata .
a survey of highly conserved duplications in the human genome of  <dig> to 98% nucleotide sequence identity and >  <dig> kb in length revealed a similar pattern.
only one-third of duplicated genes were found in clusters or segments larger than  <dig> kb, whereas the remaining two-thirds were found dispersed in the euchromatic, i.e.
gene-rich, regions all over the genome  <cit> .
our results suggest that copy number divergence between either a. halleri or a. lyrata and a. thaliana is the result of a large number of small-scale events and not of a small number of events affecting large genomic regions.
fig.
<dig> genomic distribution of copy number alterations in a. halleri and a. lyrata, relative to the a. thaliana genome.
the five chromosomes of a. thaliana are represented in different colours, with centromeric regions shaded in grey and radial axis labels indicating genomic position in mbp.
round symbols mark the positions of genes found to be copy number expanded  and reduced , respectively, in a. halleri  and a. lyrata , with radial positions scaled to the l
o
g
<dig> signal ratio of the heterologous species vs. a. thaliana.
line symbols mark the positions of 20-kb genomic regions, within which a minimum of  <dig> kb region is found to be copy number expanded  or reduced , respectively.
black line symbols mark positions of at least  <dig> consecutive concordant red or blue lines, corresponding to 40-kb genomic regions, in which a minimum of  <dig> kb is found to be copy number expanded or reduced.
the centre of the radial plot shows a zoomed  <dig> -mbp region , with red circles representing copy number expanded genes in a. halleri  and a. lyrata , labelled by short gene names or agi codes.
genes  and intergenic regions  are marked; tickmarks are spaced by  <dig> kb
discussion
cnv detection pipelines employing arrays are still generally considered to be more accurate than sequencing-based algorithms , particularly for the detection of large duplicated segments or duplications of high sequence conservation  <cit> .
de novo assembly of non-model genomes often collapses paralogous gene copies into a single locus.
while the traditional shotgun assembly resulted in a pronounced under-detection of highly identical gene cnv regions of sizes >  <dig> kb  <cit> , de novo assembly of short reads missed 99%  <cit>  of all known sequence duplications  <cit> .
by contrast, the higher the similarity between paralogs, the higher is the probability of their detection by hybridization-based methods.
this makes array-cgh methods particularly useful for detecting recent gene duplications or gene copy number divergence in emergent species.
recent and thus almost identical gene duplicates are expected to underlie the distinctive traits of species after short divergence times of  <dig> to  <dig> mya, such as human and apes, or a. halleri and a. thaliana, as in our study  <cit> .
cross-species array-cgh methods thus fill a gap in the identification of gene cnvs and should be considered an important complementary approach to other cnv detection methods and for genome annotation.
previously, array-cgh techniques were demonstrated to be quick, reliable and cost-effective in the detection of cnvs within the same species or lineage.
although whole genome tiling arrays are most widely used for this purpose, cross-species comparative genomic studies have also used gene expression arrays to compare the coding portions of genomes, for example the identification of genomic islands of speciation in three diverging populations of anopheles gambiae  <cit>  and the identification of human-specific gene duplications and contractions at a genome-wide level  <cit> .
the use of expression arrays for determining inter-species genic copy number divergence via array-cgh was generally demonstrated to give reliable and reproducible results  <cit> .
to date, affymetrix ath <dig> microarrays have been used in cross-species genomic dna hybridizations for selecting a subset of probes for use in subsequent cross-species transcriptomics studies, but not with the aim to determine cnvs  <cit> .
conversely, combined with the global scaling approach developed here, array-cgh data could be used for adjusting hybridization signals in cross-species transcriptomics in order to improve accuracy without the loss of probe information.
by comparison to previously published studies, besides being the only one in plants, our study decidedly benefits from the choice of species that are closely related to the reference species.
our method yields erroneous results when the probes on ath <dig> gene chip have been sourced from genomic regions prone to high sequence divergence such as transposable elements.
we found an enrichment of retrotransposon/transposase functional category in the cnrs, contrary to the known expansion of te families in a. lyrata relative to a. thaliana.
we calculated an average of  <dig> % sequence divergence between the ath <dig> probe sequences corresponding to a. thaliana transposable element genes and the corresponding sequences of the a. lyrata reference genome, identified as the best blast hits.
we thus propose that our result can be explained by the high degree of sequence divergence of transposable element gene families from a. thaliana in both heterologous species, as they diverged from a. thaliana at approximately the same time.
indeed, other sequence based orthology methods like ensembl-compara also report the te genes as missing orthologs in a. lyrata  <cit> .
prior to the publication of the a. lyrata genome sequence, several population genomic studies predicted a reduction in te copy number in comparison to a. thaliana  <cit> .
similarly, a well-known cne gene in a. halleri - hma <dig> was not detected by our methodology.
on the ath <dig> genechip, probe sequences for hma <dig> are located within the 3-´region of the hma <dig> coding sequence, which encodes the c-terminal cytoplasmic regulatory domain of the protein and is highly divergent in a. halleri .
this unusually high nucleotide sequence divergence explains our false negative result for hma <dig> copy number expansion.
additionally, a significant enrichment of genes “not assigned” were found among copy number reduced genes of both a. halleri and a. lyrata.
among our predictions of cnrs common to both heterologous target species, we would generally expect false positives corresponding to genes that are highly divergent from their a. thaliana orthologues and that belong to fast-evolving groups of genes.
the sensitivity of our method for copy number reductions  also appears to be lower than other methods but this is because employing the procedures of both previous studies resulted in the prediction of a substantial, more than 26-fold excess of gene copy number reductions or deletions  when compared to the high-confidence set of predictions shared by both ensembl plants  <cit>  and the a. lyrata reference genome project  <cit> .
thus, specifically accounting for inter-species sequence divergence in the analysis of cross-species array-cgh data enhances accuracy and prevents the under-reporting of true gene orthologs.
we also note that our superior results in comparison to previous methods might reflect the fact that each method is most suited to the array platform and scientific goal for which it was developed.
for example, single base substitutions are likely to be less deleterious in hybridizations to 500-bp probes on spotted cdna arrays  <cit>  than to the much shorter 25-nt probes of the ath <dig> array.
different from the ath <dig> array, tiling arrays do not require any summarization of signals from single probes within a probeset  <cit> .
extrapolating from our results that several genes known to be highly-expressed in a. halleri relative to a. thaliana are also copy number expanded, there might be a general association of gene copy number expansion with enhanced transcript levels.
indeed, one of the ways in which gene copy number divergence can confer an adaptive advantage is by effecting changes in gene product dosage, with a high probability of generating a physiological effect.
in support of this, drastic phenotypes have been demonstrated to result from modification of transcript levels resulting from alterations in gene copy number .
regarding the enrichment of mitochondrial electron transport functions in cne genes of a. halleri, a function of mitochondrial proteins in metal hypertolerance or hyperaccumulation has not been reported to date, but complex iii of the mitochondrial electron transport chain was proposed to be a major molecular target for cd toxicity in plants  <cit> .
alternatively, our observation may relate to metabolic alterations allowing enhanced dicarboxylate accumulation.
malate and citrate have been identified as quantitatively important metal chelators in hyperaccumulator plants  <cit> , thought to operate mostly inside vacuoles and in the apoplastic xylem.
an important role of these anions at storage sites might be supplying a charge balance rather than chelation.
finally, these copy number expansions could reflect a need for enhanced mitochondrial respiration to energize processes involved in metal hyperaccumulation or hypertolerance.
future experiments will have to address these alternative hypotheses on the possible roles of copy number expanded genes associated with mitochondrial electron transport and atp synthesis in metal hypertolerance or hyperaccumulation or other traits of a. halleri.
the biotic defense functions were found enriched in cnr genes of a. halleri.
the elemental defence hypothesis postulates that the extraordinarily high concentrations of metals accumulated in leaves of metal hyperaccumulator plants act as a defence against herbivory or pathogen attack  <cit> .
furthermore, the trade-off hypothesis postulates that metal hyperaccumulation-based defences allow plants to reduce the production and thus fitness costs of secondary metabolite or protein-based organic defences, incurred through their direct metabolic expense or indirect metabolic expense in scarce nutrients such as sulfur and nitrogen  <cit> .
a quantitative comparison of metabolic expenses of organic and elemental defences, however, has not been possible to date.
indeed, secondary metabolism-related functions were not found to be enriched among copy number reduced genes of a. halleri in this study.
thus, our results do not provide support for the commonly proposed metabolic cost trade-off hypothesis regarding elemental defences in hyperaccumulators.
a fitness cost of up to 9% was reported for the expression of a single r gene, rpm <dig>  in a. thaliana  <cit> .
recent work has identified plant r genes as an important group of genes underlying instances of hybrid necrosis  <cit> .
a. halleri is an obligate outcrosser, for which post-zygotic hybrid incompatibility is expected to have more severe consequences than for its largely selfing relative a. thaliana.
in metal hyperaccumulator plants, fitness costs of r gene expression could additionally involve the recently proposed inadvertent activation of r protein-mediated defence signalling by internal heavy metals  <cit> .
the loss of r genes in a. halleri can be interpreted to suggest that for a plant capable of metal hyperaccumulation, threats from biotic stress are reduced on soils permitting the hyperaccumulation of heavy metals .
our results thus suggest that cnv resulting from structural mutations are important for plant evolution occurring upon interaction with both biotic stress and abiotic environmental factors.
genetic alterations similar to those reported here for a. halleri might thus underlie other reports of niche adaptation and consequent speciation driven by herbivore pressure alone, or in combination with soil edaphic factors  <cit> .
CONCLUSIONS
in this work, we generated conservative estimates of genome-wide between-species gene copy number divergence in two close relatives of the model plant a. thaliana, based on cross-species microarray hybridization data.
we devised a novel normalization strategy, which requires less existing sequence data, yields more comprehensive results and performs better than previously developed methods.
functions in transition metal homeostasis were found to be most highly enriched among copy number expanded genes in the metal hyperaccumulator a. halleri, but not in the non-accumulating sister-species a. lyrata, despite a shared overall genome structure and an equivalent evolutionary divergence from a. thaliana.
in combination with the enrichment of biotic stress functions among copy number reduced genes of a. halleri also observed here, our finding lends support to the elemental defence hypothesis for the evolution of the metal hyperaccumulation trait  <cit> .
our results highlight the genome-wide importance of gene copy number alterations in adaptive evolution and suggest that genome scans for copy number divergence can identify functional networks that have been targets of natural selection.
we propose that our findings are applicable in ecotoxicology for identifying the types and targets of environmental change-mediated stress in suitable indicator organisms.
finally, our study has identified novel candidate genes for the future improvement of the molecular physiological understanding of metal hyperaccumulation and associated hypertolerance in plants.

 genetic and epigenetic variation form the basis for local adaptation and speciation processes, and are becoming increasingly accessible through advances in genomic and bioinformatic tools. the advent of microarray and ultra-high throughput sequencing technologies have thus brought about a renewed interest in evolutionary questions, with a prospect for gaining novel insights at the whole-genome level. these opportunities have spurred genome-wide surveys of single nucleotide polymorphisms and methylation polymorphisms in many organisms including plants, for example in multiple accessions of the genetic model organism arabidopsis thaliana and in closely related species. in attempts to identify causative genetic changes in plant adaptations, classical linkage analysis and genome-wide association studies have successfully mapped traits governing the performance under local environmental conditions to snps at specific loci. structural variation in the form of gene copy number variation polymorphism is an influential component of natural genetic diversity that markedly contributes to phenotypic variation. however, cnv has been addressed in noticeably fewer studies because of technical difficulties in its comprehensive and reliable assessment. short cnvs consisting of insertions or deletions below kb in size can be readily detected based on uhts technologies. however, the identification of cnvs comprising from kb up to one or multiple genes has generally remained challenging. genome-wide analyses in human and other mammalian model organisms revealed cnvs to be much more abundant than previously known, e.g. affecting 10% of the mouse genome and 12% of the human genome. cnvs have been implicated in human disease etiology, and evidence for adaptive cnvs is also emerging. in comparison to mammalian genomes, gene duplications and deletions especially from whole genome duplications appear to be even more abundant in plant genomes. single-gene and segmental duplications as well as whole-genome duplications have been hypothesized to propel adaptive evolution and speciation. in plants, this view is supported by recent reports on cultivar-specific boron tolerance in barley, aluminium tolerance in maize and species-wide heavy metal tolerance in the wild plant arabidopsis halleri, all supporting the role of gene copy number expansion in plant adaptation to abiotic stress. population genomic data, for example from arabidopsis thaliana and zea mays, have identified an unexpectedly high abundance of cnvs, generating interest in the contribution of structural mutations to genome plasticity. ten percent of maize genes were found to exhibit copy number polymorphisms, and an experimental evolution study in a. thaliana reported de novo structural mutations resulting in copy number variant genes after only generations. although between-species genome comparisons have remained difficult to date, the few existing studies have supported the hypothesis that gene copy number expansions, and especially those involving tandem duplications, might underlie plant adaptations to environmental stress. given that novel functions are much more likely to be generated by adaptive specialization of one of several pre-existing copies of a duplicated gene than by an entirely novel gene, such comparative studies are key to understanding the patterns of genomic polymorphisms associated with adaptation and speciation. the availability of a well-annotated genome sequence and a wealth of knowledge on gene functions for arabidopsis thaliana, as well as for several closely related species that have diverged over short evolutionary timespans, render arabidopsis a suitable model genus to study adaptation and speciation processes. one of its species is arabidopsis halleri — a wild outcrossing, zn and cd hyperaccumulating and hypertolerant species that is naturally found on both highly metal-contaminated and non-contaminated soils. its genome is expected to be about 25% larger than that of its non-accumulating, non-tolerant sister species a. lyrata, and about 65% larger than the size of the genome of a. thaliana, from which both species diverged between and million years ago. the striking phenotypic contrast between a. halleri and the two other closely related species, despite a high sequence similarity within coding regions, provides an exceptional opportunity to elucidate molecular evolutionary patterns reflecting the influence of natural soil characteristics on adaptation and speciation. fig. comparison of the metal hyperaccumulator species arabidopsis halleri to the closely related non-hyperaccumulator species a. lyrata and a. thaliana. a representative photograph is shown for each species, together with the estimated evolutionary distances separating them, given as the divergence times from a common ancestor. listed below are key phenotypic and genomic characteristics. mya, million years ago previous cross-species transcriptomics studies identified a number of differentially expressed candidate genes for the metal hyperaccumulation/hypertolerance trait of a. halleri. among these, heavy metal atpase, which encodes a metal pump that acts as an exporter of zn2+ and cd2+ from specific cell types, was shown to be necessary both for the hyperaccumulation of zn and for the full extent of zn and cd hypertolerance. strongly increased hma transcript levels in a. halleri were attributed to a lineage-specific tandem triplication combined with cis-regulatory mutations. an analysis of sequence polymorphism in the genomic region of hma gene copy number expansion demonstrated strong positive selection, as well as selection for enhanced hma gene product dosage. another candidate gene, nicotianamine synthase, encodes an enzyme that catalyses the biosynthesis of the low-molecular-weight metal chelator nicotianamine from s-adenosyl methionine, and was shown to contribute to zn hyperaccumulation. in addition to hma several other transition metal homeostasis candidate genes of a. halleri were demonstrated to be copy number expanded through the dna gel blot technique. the objective of the work presented here was to identify genes exhibiting copy number expansion in a. halleri at a genome-wide scale, in relation to the known species-specific extreme traits. we conducted a survey of gene copy number divergence in a. halleri relative to a. thaliana by employing array-comparative genomic hybridization in a cross-species manner using the ath microarray designed for a. thaliana. in order to test whether the identified cnds are species- and thus possibly trait-specific, our analysis included a. lyrata as a third species, which is a non-tolerant non-hyperaccumulator like a. thaliana, but shares with a. halleri an equal phylogenetic distance from a. thaliana. we devised a novel routine for evaluating cross-species array-cgh data, which is based on the quantification and subsequent global correction of the effects of sequence divergence on hybridization signal intensities. our procedure operates without loss of probe information, which is crucial for retaining statistical power for cnv estimation further downstream. our predictions of genic cnds were validated against a small set of genes with known copy number in a. halleri and against a set of genes predicted to be copy number expanded or reduced according to the a. lyrata reference genome sequence. gene copy number expansions in a. halleri, but not in a. lyrata, were found to be significantly enriched for metal homeostasis functions. conversely, biotic stress functions were significantly enriched among genes exhibiting copy number reduction in a. halleri, but not in a. lyrata. these results suggest that between-species divergence in gene copy numbers reflects adaptive evolution of metal hyperaccumulation, a species-specific trait of a. halleri that has been proposed to provide an elemental defence against biotic stress.  metal hyperaccumulation and hypertolerance in a. halleri have previously been attributed to the constitutively enhanced expression of a number of metal homeostasis genes, several of which were additionally shown to be expanded in genomic copy number through dna gel blots, bac sequencing or other methods. here, the technique of cross-species array-cgh was employed for a genome-wide assessment of between-species divergence in gene copy number. fragmented and labelled nuclear genomic dna samples from a. thaliana and from the two closely related heterologous species a. halleri and a. lyrata, were hybridized in duplicate to affymetrix ath genechip®; microarrays. the challenge in cross-species hybridizations of genomic dna is that the target gene sequences of the heterologous species inevitably contain mismatches relative to the probe sequences of the reference species on the microarray, thus reducing the efficiency of hybridization and resulting in lowered signals. we employed a novel approach to correct this bias by using a signal adjustment strategy, which - unlike previous methods — accounts for sequence mismatches through a global adjustment of cross-species hybridization signal intensities. in brief, we implemented a two-step normalization scheme. the first step was a conventional within-species normalization, which was applied to raw signal intensities from each pair of two replicate microarray hybridizations of the same target species. the second step was an adjustment of normalized signal intensities through the calculation and application of a species-specific global scaling factor for compensating the effects of sequence divergence from a. thaliana. fig. overview of the data analysis workflow. flowchart summarizing our two-step normalization approach for the processing of cross-species genomic hybridization data, consisting of within-species normalization and global scaling of signals through species-specific signal correction factors, followed by the final prediction of copy number divergent genes. grey arrows and backgrounds mark the auxiliary steps taken for the determination of species-specific global scaling factors with the aid of reference gene datasets to implement our strategy, we began by establishing a representative subset of probes for which curated sequence data was available from a. halleri, termed reference dataset. a similar reference dataset was also generated at random from the available reference genome sequence of a. lyrata. for each heterologous target species, the signal correction factor was calculated from the statistical distribution of the occurrence of mismatches and the ensuing effect on hybridization signal intensity as measured in the respective reference dataset. subsequently, the normalized and corrected cross-species hybridization data were analysed for differential signals between species in order to identify putative copy number divergent genes. finally, a comparison between copy number alterations in a. lyrata and a. halleri enabled us to identify species-specific copy number alterations. consequences of inter-species sequence divergence for mismatch occurrence between probe sequences and heterologous target sequences for the adjustment of microarray signals in cross-species array-cgh, we generated one reference dataset of representative, curated sequence data for each of the two heterologous target species. the a. halleri reference dataset comprised genes, yielding matching probe sequences on the microarray. because of the lack of a reference genome, these data corresponded to previously obtained sequences from a. halleri ssp. halleri. the a. lyrata reference dataset comprised genes with matching probe sequences on the microarray, obtained from the published reference genome. the number and positions of mismatches between each heterologous target sequence and the corresponding microarray probe sequence was determined. for a. halleri and a. lyrata, respectively, and 35% of all probe sequences were fully conserved across species, and 29% contained only a single mismatch with respect to the nucleotide-long probe sequence, and and 30% of sequences contained between and mismatches compared to the corresponding probe sequence. fig. frequency distribution of mismatch occurrence between microarray probe sequences and heterologous target gene sequences. shown is the percentage of a. thaliana probes on the ath array that display no mismatches up to mismatches when hybridized to non-a. thaliana genomic dna from either a. halleri or a. lyrata. the expected frequency distribution is shown by the grey line, and was calculated based on the average coding sequence identity within transcribed regions between a. halleri and a. thaliana we computed the expected distribution of the total number of mismatches of a given nucleotide sequence with respect to the probe sequence on the microarray. the expected distribution of mismatches in cross-species hybridization of a. halleri gdna closely matched our observations for the reference dataset. both the observed and expected probabilities of the occurrence of more than mismatches by comparison to the probe sequence were negligible. finally, the observed mismatch distributions for a. halleri and a. lyrata were highly similar to each other, thus confirming similar levels of sequence divergence from a. thaliana and allowing us to use the expected distribution of mismatches calculated for a. halleri also in a. lyrata hybridizations. both the a. halleri and a. lyrata lineages are thought to have diverged from the common ancestor with a. thaliana at the same point in the past. our observations support the theory that a correlation exists between the levels of sequence divergence and actual phylogenetic distances between species, as was estimated, for example, based on cross-species array-cgh data. quantification of effects of sequence mismatches on signal intensities in cross-species microarray hybridization sequence mismatches are known to be the single most confounding factor biasing the signals of cross-species array hybridizations. a previous study has estimated sequence mismatches to account for at least 40% of the average noise in microarray hybridization, and several studies confirm mismatches as the primary cause of failure of conventional normalization techniques in cross-species microarray data analysis. as a result of sequence divergence, sequence mismatches are expected to reduce the hybridization efficiency of genomic dna from a. halleri and a. lyrata to the a. thaliana probe sequences on the ath microarray, resulting in lowered overall hybridization signal intensity. after background correction of raw data, we examined the influence of the total number and positions of mismatches on the normalized hybridization signal intensities using the probe signal intensities from our a. halleri and a. lyrata reference datasets. as expected, hybridization signal intensity decreased with increasing number of mismatches in a probe. the largest decrease in signal intensity by and 40% in a. halleri and a. lyrata, respectively, was observed for a single mismatch. additional mismatches had only small effects, with a total of four mismatches resulting in a further reduction of probe signal intensity by and 7% in a. halleri and a. lyrata, respectively. there were only minor differences between a. halleri and a. lyrata in the dependence of signal intensity on the number of mismatches. note that the published a. lyrata reference genome, which was used here to determine the number of mismatches between a. lyrata reference dataset target sequences and the corresponding ath probe sequences, was from the north american subspecies lyrata, whereas our hybridization experiments were conducted with the european ssp. petrea. a stark reduction in diversity has been reported in a. lyrata ssp. lyrata by comparison to a. lyrata ssp. petrea, and several studies report differentiation between the two sub-species that have been isolated from each other for between and thousand years. fig. dependence of hybridization signal intensity on number and position of mismatches with respect to the probe sequence on the ath array. a values are arithmetic means of background-corrected raw probe signal intensity ratios for non-a. thaliana gdna relative to a. thaliana gdna hybridizations, shown as a function of the total number of mismatches of the heterologous target sequence compared to the corresponding a. thaliana 25-mer probe sequence. b independence of hybridization signal intensity from the position of a single mismatch with respect to the probe sequence. values are arithmetic means of background-corrected raw probe signal intensity ratios for non-a. thaliana gdna relative to a. thaliana gdna hybridizations, shown as a function of mismatch position in the heterologous target sequence compared to the corresponding a. thaliana probe sequence. black circles represent the representative a. halleri reference dataset; white diamonds represent the representative a. lyrata reference dataset surprisingly, we observed a noisy profile of signal intensity over different positions of a single mismatch along the probe sequence instead. the expected sharp drop in signal intensity when a single mismatch is positioned in the centre of a probe sequence, as proposed by affymetrix for so-called mismatch probes, was not detected here. this finding is in agreement with a number of previous studies, which have pointed out that experimental data do not conform to this postulate and that, in fact, for some probes signal intensity was even found to be higher for mm probes than for perfectly matching probes. consequently, position-based effects on hybridization signal intensity are hard to construct, and accordingly, the most popular normalization methods no longer take the information from mm probes into account. therefore, for our between-species normalization strategy, we did not consider the influence of sequence mismatch position. we estimated the incremental signal correction factor s k for a probe with k mismatches as the average of the ratio of normalized hybridization signal intensity of a. thaliana to the respective signal intensity of the heterologous species. for each probe containing to mismatches, incremental signal correction factors were weighted by their probability of occurrence, followed by the calculation of the arithmetic mean to yield species-specific global scaling factors. these global signal correction factors of for hybridizations of a. halleri gdna and for hybridizations of a. lyrata gdna were employed to scale the hybridization signal intensities of the respective cross-species microarray hybridizations. cross-species normalization and validation of copy number divergent genes the median raw signal intensities for the heterologous species a. halleri and a. lyrata were lower than those for the ath target model species a. thaliana, namely by and 36%, respectively. after applying conventional within-species vsn normalizations, median normalized signal intensities were more uniform across replicates within each species. however, the differences between species were large, with median signal intensities for a. halleri and a. lyrata which were and 49% lower, respectively, than for a. thaliana. upon the subsequent between-species scaling of vsn-normalized signal intensities from the a. halleri and a. lyrata hybridizations employing species-specific global scaling factors, median signal intensities became more similar across species and remained 17% lower for a. halleri and 11% lower for a. lyrata, respectively, than for a. thaliana. following normalization and scaling of hybridization signals, and copy number expanded genes were identified in a. halleri and a. lyrata, respectively. furthermore, genes predicted to be copy number reduced were identified in a. halleri and in a. lyrata. overall, cnes and cnrs are shared between a. halleri and a. lyrata compared to the reference species a. thaliana. thus, based on a. thaliana as a reference, an about 3-fold larger number of copy number divergent genes was detected in a. halleri than in a. lyrata, whereas nucleotide sequence divergence from a. thaliana was similar in both heterologous species within transcribed regions. the observed difference between a. halleri and a. lyrata was not merely a spurious result caused by a higher level of polymorphism between the two replicate a. halleri samples than between the a. lyrata replicates. this was confirmed by performing all data processing steps with two additional pairs of a. halleri replicates, each consisting of either one of the two single a. halleri hybridizations and one a. halleri replicate generated in silico by simulating between-replicate variation as observed in a. lyrata, respectively. consequently, our results suggest that the rate of either acquisition or maintenance of gene copy number changes in the genome can differ between closely-related lineages or species. this is in stark contrast to the general stability of base substitution rates normalized to genome size and generation. fig. distribution of signal intensities before and after normalization and scaling. boxplot of background-corrected raw hybridization signal intensities, normalized signal intensities after vsn normalization of the replicate arrays of each species, respectively, and signal intensities after the application of species-specific global scaling factors to the normalized data. boxes show median, and upper and lower quartiles, of l o g probe signal intensities for each gdna hybridization. upper and lower horizontal bars mark all values lying within times the inter-quartile range. replicate hybridizations are denoted and and grouped by species to evaluate the reliability of our predicted cnds, we compared our results to genes of known copy number status. for a. halleri, we used a set of genes. the evaluation of microarray-based predictions of cross-species gene cnd against known copy number status indicated % specificity, % precision and % sensitivity of our cross-species array-cgh based cnd estimation. for a. lyrata, a complete reference genome sequence is available. this provides an opportunity for more extensive data validation by comparing our predictions of gene cnds with predictions based on the reference genome sequence. orthology predictions retrieved from ensembl plants indicated that, relative to all a. thaliana genes represented on the ath genechip, orthologous genes of a. lyrata are copy number expanded, whereas genes are reduced in copy number or deleted. from this set of genes, we further chose the subset of highly conserved multi-copy genes of ≥ 95% average sequence identity between transcribed regions of all paralogs with respect to their a. thaliana ortholog. our final reference dataset of highly conserved copy number expanded genes in a. lyrata, as predicted based on ensembl plants, contained genes. by comparison, our method yielded % specificity, % precision and % sensitivity. these scores suggest that our procedure is conservative, making predictions with high reliability, and sacrificing sensitivity for a higher specificity. a direct comparative assessment of the performance of our array-cgh based method for a. halleri and a. lyrata is not possible because of the differing qualities of experimental evidence underlying the two datasets available for validation. darby et al. 2011d machado et al. 2010e a a. lyrata genes sharing ≥ 95% sequence identity with their closest a. thaliana homologue are termed highly conserved bheaders of half-columns refer to total number of cnes predicted by ensembl plants alone or additionally by a. lyrata genome analysis, respectively, as given in parentheses here. shown are commonalities with these two groups of genes or data referring to these two groups of genes ctrue positive is a cne detected based on array-cgh that was previously predicted to be a cne by ensembl plants alone, or additionally by a. lyrata genome analysis d e we compared the performance of our array-cgh based approach with that of the two previous studies that also aimed at estimating cnd using the array-cgh technique. we reproduced the normalization and scaling strategies of machado and renn and darby et al. as described, with few small modifications necessary to apply these methods to our array-cgh platform. the method of darby et al. resulted in the prediction of a -fold elevated number of gene copy number expansions. out of the two previously published methods, maximum sensitivity, specificity and precision of the detection of copy number expansion among highly conserved genes were, and %, respectively, all inferior to our method. even for the genes that are not highly conserved but predicted to be copy number expanded concordantly by both ensembl plants and the a. lyrata genome project, our method reports higher sensitivity, specificity and precision – and % respectively than previous studies –, and %. specificity and precision of our method were also superior concerning copy number reductions or gene deletions. functional analysis of copy number divergent genes of a. halleri after identifying the sets of genes exhibiting copy number divergence by comparison to the reference species a. thaliana in either of the two heterologous species according to array-cgh, we evaluated these for any enrichment of functional categories using the mapman ontology. copy number expanded genes of a. halleri showed a statistically significant enrichment for transition metal homeostasis-related gene functions, and for mitochondrial electron transport/atp synthesis-related functions, relative to all genes represented on the array. among copy number reduced genes of a. halleri, there was a significant enrichment for biotic stress-related functions. in contrast, none of these functional categories was detected to be significantly enriched among genes divergent in copy number in a. lyrata by comparison to a. thaliana. in the light of the two extreme traits specific to a. halleri — namely zn/cd hyperaccumulation and associated hypertolerance — the genome-wide overrepresentation of metal homeostasis genes among copy number expanded genes is remarkable. indeed, a high occurrence of copy number expansion was reported among arabidopsis halleri metal homeostasis candidate genes identified through the presence of elevated transcript levels in a. halleri compared to a. thaliana. moreover, copy number expansion is known to contribute to high hma transcript levels in a. halleri, which in turn are necessary for both metal hyperaccumulation and the full extent of metal hypertolerance. hma gene copy number expansion is not limited to a. halleri, but also found in the zn/cd hyperaccumulator species noccaea caerulescens, similarly associated with strongly elevated transcript levels. a. halleri mtp is another copy number-expanded candidate gene, for which several lines of evidence suggest an involvement in zn hypertolerance. it was not known to date whether these findings on individual candidate genes pertain at the genome-wide level, but this is now supported by array-cgh data presented here. our data additionally confirm the previous finding of zip copy number expansion. in contrast, our array-cgh analysis did not detect hma as copy number expanded in a. halleri, although this is well established. one of the transition metal homeostasis candidate genes newly identified to be copy number expanded in a. halleri is nas which was demonstrated to be highly expressed in roots of a. halleri and to contribute to zn hyperaccumulation. array-cgh also predicts ahhma to be copy number expanded. this candidate gene was reported as highly expressed in a. halleri, and an ahhma cdna confers zn and cd tolerance upon heterologous expression in yeast. finally, atpcr and atpcr have been implicated in the export of zn and cd, respectively, from cells and appear to be copy number expanded in a. halleri. indeed, transcript levels were found to be higher in a. halleri than in a. thaliana in a previous microarray hybridization study, but this was not explicitly reported because of an ambiguous assignment of the probe set to this pair of highly similar a. thaliana genes. sams encodes an enzyme that catalyses the biosynthesis of the substrate for nicotianamine synthase and was previously identified to be more highly expressed in a. halleri than in a. thaliana; and it is now identified as copy number expanded. expanding from these findings, the enrichment of mitochondrial electron transport functions among the genes copy number expanded in a. halleri was surprising. a majority of copy number expansions, half of which are shared with a. lyrata, affect complex i acting as nadh ubiquinone oxidoreductase in oxidative phosphorylation. a more parsimonious explanation for the common expansions could be that these genes have been deleted in a. thaliana relative to the common ancestor. the two genes privately copy number expanded in a. halleri have predicted roles in atp synthesis and cytochrome c biogenesis, respectively, and one is annotated as a mitochondrial dicarboxylate carrier. fig. significantly enriched functional classes among copy number expanded and copy number reduced genes of a. halleri. bars/brackets represent the proportion/number of genes in each functional class among all copy number expanded or reduced genes of a. halleri or a. lyrata, and among all nuclear genes represented on the array for comparison. asterisks mark statistically significant enrichment by comparison to the set of nuclear genes represented on the ath microarray. the analysis was based on mapman functional categories. shown are only top-level categories of the ontological hierarchy a. halleri vs. a. thaliana a. lyrata vs. a. thaliana l o g fcb l o g fcb mfm1 pcr1; pcr2 chx25 261448_atc itp2 chx2 267304_atd zip6 fer4 266718_atc,d mtp1 mt2a fer3 cch sams2 as8 253658_atc,d hma3 hcc2 mt2b 248048_atc nas2 ysl4 zip3 ptac17 ti1 250944_atd hipp6 hipp43 nas-like fp3 copt emb1513 259871_atc 256930_atd oasa2 253413_atc zip9 pdf sams3/ mat4* mtp7 260913_atd sams1 chx16 257365_x_atc,d pdf b 249255_atc chx18 263838_atc,d mat3** hipp27 cax11 252694_atc pdf zip11 mt3 258987_atd hcc1 a. halleri vs. a. thaliana a. lyrata vs. a. thaliana l o g fcb l o g fcb b12-1 dic2 b atpq 249627_atd emb1467 ccme hcc2 phb6; phb1 b14 a. halleri vs. a. thaliana a. lyrata vs. a. thaliana l o g fcb l o g fcb – – – – rlp40; rlp42 ttr1; wrky16 rlp34; rlp53 transition metal homeostasis. mitochondrial electron transport/atp synthesis. biotic stress mapman functional classes ashort names of genes previously reported to be copy-number expanded in a. halleri based on other experimental approaches are given in bold type b l o g fold-change is underlined for genes that are concordantly copy number expanded/reduced in both a. halleri and a. lyrata chigher transcript abundance in roots of a. halleri vs. a. thaliana dhigher transcript abundance in shoots of a. halleri vs. a. thaliana ∞synthesis of cytochrome c oxidase; $zinc-regulated transporter, iron-regulated transporter-like; *at3g was annotated as sams in ref. talke et al. 2006; **at2g was annotated as sams in talke et al. 2006; &nadh-ubiquinone oxidoreductase, mitochondrial respiratory chain complex i;#coiled coil or toll-interleukin receptor -nucleotide binding site -leucine rich repeat;+nucleotide-binding motif-containing domain shared by human apaf- certain plant r gene products, and caenorhabditis elegans ced- among the genes predicted to be copy number reduced according to our array-cgh results, we observed an overrepresentation of biotic defence functions. this is an interesting observation with respect to the proposed ecological role of metal hyperaccumulation in plants. almost all biotic stress-related genes predicted to be copy number reduced in a. halleri encode members of large protein families typically involved in plant innate immunity and designated as disease resistance genes, such as the predominating tir-nbs-lrr receptor kinases. infecting pathogens generate characteristic molecular patterns that can be specifically recognized by cognate r gene products, which subsequently trigger a localized cell death response that is essential for plant disease resistance. the enrichment of r gene-related biotic stress functions among genes reduced in copy number in the hyperaccumulator a. halleri supports the elemental defence hypothesis as well as the trade-off hypothesis for the evolutionary role of elemental defence. accordingly, elemental defences through metal hyperaccumulation allow for the loss of r genes, thus alleviating the fitness costs associated with r gene expression. at lower levels of the ontological hierarchy, post-translational modification functions, and genes encoding glycine-rich and crinkly-like proteins were over-represented among the copy number expanded genes of a. halleri, whereas dc domain-containing, as well as protein degradation/e3-scf-f-box and cysteine protease functions were significantly enriched among the genes predicted to be copy number reduced. furthermore, genes predicted to be copy number reduced in both a. halleri and a. lyrata according to our array-cgh, showed a significant enrichment of the mapman functional category or bin “dna.synthesis/chromatin structure.retrotransposon/transposase”. contrary to this observation, it is known that transposable element gene families have undergone an expansion in the a. lyrata genome relative to a. thaliana. contributions of segmental duplications or deletions to copy number divergence next we tested whether large segmental duplications or deletions make major contributions to cnd from a. thaliana in a. halleri or a. lyrata. we detected only few large segmental duplications and deletions comprising more than kb, by using a sliding window approach and the set parameters. a total of six predicted segmental deletion regions, located in positions corresponding to a. thaliana chromosomes to were shared between a. halleri and a. lyrata. most detected segmental deletions were located in or around regions corresponding to a. thaliana centromers. these segmental deletions were composed almost entirely of transposable element genes and are thus likely to include regions of low sequence conservation rather than segmental deletions. the sequences of centromeric regions are known to be highly divergent in both nucleotide composition and length even between closely related species. none of the mapman functional classes of top level ontological hierarchy that were significantly enriched among copy number divergent genes was found similarly enriched among the genes contained in segmental duplications or deletions. overall, cnds, and especially cnes, appeared to be dispersed evenly in the genomes of both a. halleri and a. lyrata. a survey of highly conserved duplications in the human genome of to 98% nucleotide sequence identity and > kb in length revealed a similar pattern. only one-third of duplicated genes were found in clusters or segments larger than kb, whereas the remaining two-thirds were found dispersed in the euchromatic, i.e. gene-rich, regions all over the genome. our results suggest that copy number divergence between either a. halleri or a. lyrata and a. thaliana is the result of a large number of small-scale events and not of a small number of events affecting large genomic regions. fig. genomic distribution of copy number alterations in a. halleri and a. lyrata, relative to the a. thaliana genome. the five chromosomes of a. thaliana are represented in different colours, with centromeric regions shaded in grey and radial axis labels indicating genomic position in mbp. round symbols mark the positions of genes found to be copy number expanded and reduced, respectively, in a. halleri and a. lyrata, with radial positions scaled to the l o g signal ratio of the heterologous species vs. a. thaliana. line symbols mark the positions of 20-kb genomic regions, within which a minimum of kb region is found to be copy number expanded or reduced, respectively. black line symbols mark positions of at least consecutive concordant red or blue lines, corresponding to 40-kb genomic regions, in which a minimum of kb is found to be copy number expanded or reduced. the centre of the radial plot shows a zoomed -mbp region, with red circles representing copy number expanded genes in a. halleri and a. lyrata, labelled by short gene names or agi codes. genes and intergenic regions are marked; tickmarks are spaced by kb discussion cnv detection pipelines employing arrays are still generally considered to be more accurate than sequencing-based algorithms, particularly for the detection of large duplicated segments or duplications of high sequence conservation. de novo assembly of non-model genomes often collapses paralogous gene copies into a single locus. while the traditional shotgun assembly resulted in a pronounced under-detection of highly identical gene cnv regions of sizes > kb, de novo assembly of short reads missed 99% of all known sequence duplications. by contrast, the higher the similarity between paralogs, the higher is the probability of their detection by hybridization-based methods. this makes array-cgh methods particularly useful for detecting recent gene duplications or gene copy number divergence in emergent species. recent and thus almost identical gene duplicates are expected to underlie the distinctive traits of species after short divergence times of to mya, such as human and apes, or a. halleri and a. thaliana, as in our study. cross-species array-cgh methods thus fill a gap in the identification of gene cnvs and should be considered an important complementary approach to other cnv detection methods and for genome annotation. previously, array-cgh techniques were demonstrated to be quick, reliable and cost-effective in the detection of cnvs within the same species or lineage. although whole genome tiling arrays are most widely used for this purpose, cross-species comparative genomic studies have also used gene expression arrays to compare the coding portions of genomes, for example the identification of genomic islands of speciation in three diverging populations of anopheles gambiae and the identification of human-specific gene duplications and contractions at a genome-wide level. the use of expression arrays for determining inter-species genic copy number divergence via array-cgh was generally demonstrated to give reliable and reproducible results. to date, affymetrix ath microarrays have been used in cross-species genomic dna hybridizations for selecting a subset of probes for use in subsequent cross-species transcriptomics studies, but not with the aim to determine cnvs. conversely, combined with the global scaling approach developed here, array-cgh data could be used for adjusting hybridization signals in cross-species transcriptomics in order to improve accuracy without the loss of probe information. by comparison to previously published studies, besides being the only one in plants, our study decidedly benefits from the choice of species that are closely related to the reference species. our method yields erroneous results when the probes on ath gene chip have been sourced from genomic regions prone to high sequence divergence such as transposable elements. we found an enrichment of retrotransposon/transposase functional category in the cnrs, contrary to the known expansion of te families in a. lyrata relative to a. thaliana. we calculated an average of % sequence divergence between the ath probe sequences corresponding to a. thaliana transposable element genes and the corresponding sequences of the a. lyrata reference genome, identified as the best blast hits. we thus propose that our result can be explained by the high degree of sequence divergence of transposable element gene families from a. thaliana in both heterologous species, as they diverged from a. thaliana at approximately the same time. indeed, other sequence based orthology methods like ensembl-compara also report the te genes as missing orthologs in a. lyrata. prior to the publication of the a. lyrata genome sequence, several population genomic studies predicted a reduction in te copy number in comparison to a. thaliana. similarly, a well-known cne gene in a. halleri - hma was not detected by our methodology. on the ath genechip, probe sequences for hma are located within the 3-´region of the hma coding sequence, which encodes the c-terminal cytoplasmic regulatory domain of the protein and is highly divergent in a. halleri. this unusually high nucleotide sequence divergence explains our false negative result for hma copy number expansion. additionally, a significant enrichment of genes “not assigned” were found among copy number reduced genes of both a. halleri and a. lyrata. among our predictions of cnrs common to both heterologous target species, we would generally expect false positives corresponding to genes that are highly divergent from their a. thaliana orthologues and that belong to fast-evolving groups of genes. the sensitivity of our method for copy number reductions also appears to be lower than other methods but this is because employing the procedures of both previous studies resulted in the prediction of a substantial, more than 26-fold excess of gene copy number reductions or deletions when compared to the high-confidence set of predictions shared by both ensembl plants and the a. lyrata reference genome project. thus, specifically accounting for inter-species sequence divergence in the analysis of cross-species array-cgh data enhances accuracy and prevents the under-reporting of true gene orthologs. we also note that our superior results in comparison to previous methods might reflect the fact that each method is most suited to the array platform and scientific goal for which it was developed. for example, single base substitutions are likely to be less deleterious in hybridizations to 500-bp probes on spotted cdna arrays than to the much shorter 25-nt probes of the ath array. different from the ath array, tiling arrays do not require any summarization of signals from single probes within a probeset. extrapolating from our results that several genes known to be highly-expressed in a. halleri relative to a. thaliana are also copy number expanded, there might be a general association of gene copy number expansion with enhanced transcript levels. indeed, one of the ways in which gene copy number divergence can confer an adaptive advantage is by effecting changes in gene product dosage, with a high probability of generating a physiological effect. in support of this, drastic phenotypes have been demonstrated to result from modification of transcript levels resulting from alterations in gene copy number. regarding the enrichment of mitochondrial electron transport functions in cne genes of a. halleri, a function of mitochondrial proteins in metal hypertolerance or hyperaccumulation has not been reported to date, but complex iii of the mitochondrial electron transport chain was proposed to be a major molecular target for cd toxicity in plants. alternatively, our observation may relate to metabolic alterations allowing enhanced dicarboxylate accumulation. malate and citrate have been identified as quantitatively important metal chelators in hyperaccumulator plants, thought to operate mostly inside vacuoles and in the apoplastic xylem. an important role of these anions at storage sites might be supplying a charge balance rather than chelation. finally, these copy number expansions could reflect a need for enhanced mitochondrial respiration to energize processes involved in metal hyperaccumulation or hypertolerance. future experiments will have to address these alternative hypotheses on the possible roles of copy number expanded genes associated with mitochondrial electron transport and atp synthesis in metal hypertolerance or hyperaccumulation or other traits of a. halleri. the biotic defense functions were found enriched in cnr genes of a. halleri. the elemental defence hypothesis postulates that the extraordinarily high concentrations of metals accumulated in leaves of metal hyperaccumulator plants act as a defence against herbivory or pathogen attack. furthermore, the trade-off hypothesis postulates that metal hyperaccumulation-based defences allow plants to reduce the production and thus fitness costs of secondary metabolite or protein-based organic defences, incurred through their direct metabolic expense or indirect metabolic expense in scarce nutrients such as sulfur and nitrogen. a quantitative comparison of metabolic expenses of organic and elemental defences, however, has not been possible to date. indeed, secondary metabolism-related functions were not found to be enriched among copy number reduced genes of a. halleri in this study. thus, our results do not provide support for the commonly proposed metabolic cost trade-off hypothesis regarding elemental defences in hyperaccumulators. a fitness cost of up to 9% was reported for the expression of a single r gene, rpm in a. thaliana. recent work has identified plant r genes as an important group of genes underlying instances of hybrid necrosis. a. halleri is an obligate outcrosser, for which post-zygotic hybrid incompatibility is expected to have more severe consequences than for its largely selfing relative a. thaliana. in metal hyperaccumulator plants, fitness costs of r gene expression could additionally involve the recently proposed inadvertent activation of r protein-mediated defence signalling by internal heavy metals. the loss of r genes in a. halleri can be interpreted to suggest that for a plant capable of metal hyperaccumulation, threats from biotic stress are reduced on soils permitting the hyperaccumulation of heavy metals. our results thus suggest that cnv resulting from structural mutations are important for plant evolution occurring upon interaction with both biotic stress and abiotic environmental factors. genetic alterations similar to those reported here for a. halleri might thus underlie other reports of niche adaptation and consequent speciation driven by herbivore pressure alone, or in combination with soil edaphic factors.  in this work, we generated conservative estimates of genome-wide between-species gene copy number divergence in two close relatives of the model plant a. thaliana, based on cross-species microarray hybridization data. we devised a novel normalization strategy, which requires less existing sequence data, yields more comprehensive results and performs better than previously developed methods. functions in transition metal homeostasis were found to be most highly enriched among copy number expanded genes in the metal hyperaccumulator a. halleri, but not in the non-accumulating sister-species a. lyrata, despite a shared overall genome structure and an equivalent evolutionary divergence from a. thaliana. in combination with the enrichment of biotic stress functions among copy number reduced genes of a. halleri also observed here, our finding lends support to the elemental defence hypothesis for the evolution of the metal hyperaccumulation trait. our results highlight the genome-wide importance of gene copy number alterations in adaptive evolution and suggest that genome scans for copy number divergence can identify functional networks that have been targets of natural selection. we propose that our findings are applicable in ecotoxicology for identifying the types and targets of environmental change-mediated stress in suitable indicator organisms. finally, our study has identified novel candidate genes for the future improvement of the molecular physiological understanding of metal hyperaccumulation and associated hypertolerance in plants. 
for a genome-wide assessment of the role of copy number divergence in lineage-specific environmental adaptation, we conducted cross-species array comparative genome hybridizations of three plant species and developed a global signal scaling procedure to adjust for sequence divergence. gene copy number divergence between species is a form of genetic polymorphism that contributes significantly to both genome size and phenotypic variation. in a. halleri, transition metal homeostasis functions are enriched twofold among the genes detected as copy number expanded. our results provide genome-wide support for a link between evolutionary adaptation and cnd in a. halleri as shown previously for heavy metal atpase moreover our results support the hypothesis that elemental defences, which result from the hyperaccumulation of toxic metals, allow the reduction of classical defences against biotic stress as a trade-off. in plants, copy number expansions of single genes were implicated in cultivar- or species-specific tolerance of high levels of soil boron, aluminium or calamine-type heavy metals, respectively. moreover, biotic stress functions including mostly disease resistance gene-related genes are enriched twofold among genes detected as copy number reduced, when compared to the abundance of these functions among all genes.
 gene copy number divergence between species is a form of genetic polymorphism that contributes significantly to both genome size and phenotypic variation. in plants, copy number expansions of single genes were implicated in cultivar- or species-specific tolerance of high levels of soil boron, aluminium or calamine-type heavy metals, respectively. arabidopsis halleri is a zinc- and cadmium-hyperaccumulating extremophile species capable of growing on heavy-metal contaminated, toxic soils. in contrast, its non-accumulating sister species a. lyrata and the closely related reference model species a. thaliana exhibit merely basal metal tolerance.  for a genome-wide assessment of the role of copy number divergence in lineage-specific environmental adaptation, we conducted cross-species array comparative genome hybridizations of three plant species and developed a global signal scaling procedure to adjust for sequence divergence. in a. halleri, transition metal homeostasis functions are enriched twofold among the genes detected as copy number expanded. moreover, biotic stress functions including mostly disease resistance gene-related genes are enriched twofold among genes detected as copy number reduced, when compared to the abundance of these functions among all genes.  our results provide genome-wide support for a link between evolutionary adaptation and cnd in a. halleri as shown previously for heavy metal atpase moreover our results support the hypothesis that elemental defences, which result from the hyperaccumulation of toxic metals, allow the reduction of classical defences against biotic stress as a trade-off. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords cross-speciesarray-cghmetal hyperaccumulationcnvarabidopsis halleritoll-interleukin receptor-nucleotide binding site-leucine rich repeat protein familyresistance genes 15th international conference on bioinformatics incob 2016queenstown, singapore 21- september http://incob apbionet.org/issue-copyright-statement© the author 2016 