|Home | About | Journals | Submit | Contact Us | Français|
Allotetraploids carry pairs of diverged homoeologs for most genes. With the genome doubled in size, the number of putative interactions is enormous. This poses challenges on how to coordinate the two disparate genomes, and creates opportunities by enhancing the phenotypic variation. New combinations of alleles co-adapt and respond to new environmental pressures. Three stages of the allopolyploidization process - parental species divergence, hybridization, and genome duplication - have been well analyzed. The last stage of evolutionary adjustments remains mysterious.
Homoeolog-specific retention and use were analyzed in Arabidopsis suecica (As), a species derived from A. thaliana (At) and A. arenosa (Aa) in a single event 12,000 to 300,000 years ago. We used 405,466 diagnostic features on tiling microarrays to recognize At and Aa contributions to the As genome and transcriptome: 324 genes lacked Aa contributions and 614 genes lacked At contributions within As. In leaf tissues, 3,458 genes preferentially expressed At homoeologs while 4,150 favored Aa homoeologs. These patterns were validated with resequencing. Genes with preferential use of Aa homoeologs were enriched for expression functions, consistent with the dominance of Aa transcription. Heterologous networks - mixed from At and Aa transcripts - were underrepresented.
Thousands of deleted and silenced homoeologs in the genome of As were identified. Since heterologous networks may be compromised by interspecies incompatibilities, these networks evolve co-biases, expressing either only Aa or only At homoeologs. This progressive change towards predominantly pure parental networks might contribute to phenotypic variability and plasticity, and enable the species to exploit a larger range of environments.
An allotetraploid is formed when diploids from two different species, which may have diverged for millions of years, hybridize. The resulting plant, if viable, might have a competitive edge, such as broader ecological tolerance compared to its parents [1-3]. The evolutionary importance of polyploidy, of which allotetraploidy is a common form, is reflected in its prevalence in flowering plants : ancient polyploidy is apparent in all plant genomes sequenced to date and is estimated to have been involved in 15% of all plant speciation events . Furthermore, most cultivated crops have undergone polyploidization during their ancestry [5,6]. Why are polyploids so evolutionarily, ecologically, and agriculturally successful? To answer this question, one has to consider the evolutionary and genetic processes acting at different stages of polyploidization.
Allopolyploidization can be characterized by four distinct stages. Stage 1 is the divergence between parental species, with both species adapting to specific environments and adopting their own mating strategies and reproductive schedules. Directional selection can contribute to the fixation of species-specific beneficial mutations in coding and regulatory regions [7,8], while slightly deleterious mutations are introduced due to drift. In stages 2 and 3, the diverged species hybridize and increase ploidy, with the two events sometimes reversed in order . This change in ploidy enables the correct pairing at meiosis. Hybridization frequently results in phenotypic instability, widespread genomic rearrangements, epigenetic silencing, and unusual splicing [3,10-25]. Newly created polyploids often experience rapid intragenomic adjustments. Stages 2 and 3 are well-studied with artificial polyploids constructed in the laboratory [10,12-17,19,22-24] or spontaneously arising in nature [14,26].
Stage 4 is the long term evolution of homoeologous genes (that is, homologous genes from two parents joined into one polyploid genome and stably inherited). This stage occurs much slower on the evolutionary time-scale and has received considerably less attention, perhaps due to several technical limitations. Sequence analyses have historically required extensive cloning and bioinformatics. Microarrays have had to be specifically designed to distinguish between homoeologs and orthologs. Interesting patterns have been reported, but typically for a few genes [14,27-29]. Notably, the retention and expression of homoeologs is frequently biased towards one parental species. These patterns were reported on a large scale for approximately 1,400 out of 42,000 genes in cotton [30-32], and for dozens in Tragopogon . Recent studies have also discovered abundant genetic variation among independently originated or evolved accessions of Tragopogon [34-36]. What molecular evolutionary processes account for this variation among accessions? How does intraspecific variation in polyploid genomes contribute to phenotypic variation? These questions remain wide open.
Here, we focus on Arabidopsis suecica (As), a highly selfing species  found mainly in central Sweden and southern Finland . As originated 12,000 to 300,000 years ago (KYA) from a cross between a largely homozygous ovule-parent Arabidopsis thaliana (At, 2n10) and a pollen-parent Arabidopsis arenosa (Aa, 2n = 16) [39-41]. A single origin of As (2n = 26) has been established with mitochondrial, chloroplast, and nuclear DNA [39-41]. As originated south of the ice cover and spread north when the ice retreated 10,000 years ago . At is an annual, weedy, and mostly autogamous species native to Europe and central Asia but naturalized worldwide . It has undergone at least two rounds of ancient polyploidization  and is annotated with 39 thousand genes. Aa is a self-incompatible member of the Arabidopsis genus, carrying the highest level of genetic diversity among the species group . At and Aa diverged approximately 5 million years ago .
One can generate an artificial F1 allotetraploid (F1As) in the lab by performing a cross between a tetraploid At ovule-parent and a tetraploid Aa pollen donor. The resulting primary species hybrid contains two genomes from At and two from Aa. We can use this as an estimate, as the exact haplotypes that contributed to the initial hybridization event are not available, of the genomic composition and homoeolog-specific expression at the time of allopolyploid speciation [24,45,46]. Taking these patterns as reflective of the As ancestral state, we observed how evolution has shaped the As genome. As At is a selfer and Aa an outcrosser, At-originated homoeologs might have possessed more deleterious mutations due to Hill-Robertson interference . Are Aa-originated homoeologs more commonly retained? At and Aa evolved orthologous networks in which genes were finely tuned to coordinate, separately within each species. Interference of At and Aa homoeologs may cause mis-regulation within mixed As networks. This is akin to Dobzhansky-Muller incompatibilities . Do heterologous networks evolve to restore their original orthologous-like compositions? Here, we address these and other questions.
For every gene in As, we set to determine whether both At and Aa homoeologs are present in the genome and whether they are expressed evenly or in homoeolog-specific fashion . With the genome-wide Arabidopsis tiling microarray, we scanned the genomes of At, Aa, As, and F1As. We analyzed the transcriptome of As with tiling arrays and validated results with Illumina resequencing. We assembled a statistical pipeline to identify At and Aa homoeolog-originated signals, and to estimate their contribution to the As populations of DNA and RNA.
The Arabidopsis array features 3.2 million 25-base-long probes tiled throughout the complete genome at a 35-base distance. As these features are homologous to the At reference, they should, on average, exhibit a lower hybridization with Aa DNA. Probe intensities confirm this expectation. Two typical examples are shown for chromosomes 3 and 4 (Figures (Figures11 and and2;2; see Additional files 1, 2, 3,4, 5 and 6 for other examples). F1As signals are a sharp intermediate between At and Aa. As shows remarkable correspondence with F1As, with the exception of several extended regions. We hypothesize that these regions correspond to historic losses of homoeologous chromosomal regions in As.
We mapped features onto the genes and compared intensities between As and F1As; 6,790 genes exhibited differential hybridization (Wilcoxon ranked sum test, false discovery rate (FDR) <0.05). To identify large putative alterations, we scanned for clusters containing at least 30 genes with a strong unidirectional bias (at least 27 with the same bias, significant for at least 9 genes). We identified 39 clusters, encompassing 1,643 genes (Table (Table1).1). Some clusters were due to differential abundance of transposable-element-like sequences. Chr1 13.66 M, Chr1 14.00 M, Chr3 12.44 M, Chr3 13.36 M, and Chr5 11.06 M mainly consisted of copia-like, gypsy-like, or CACTA-like retrotransposons. Other regions - for instance, on Chr1 0.29 M, Chr3 0.30 M, Chr3 5.58 M, Chr3 21.60 M, and Chr3 22.98 M - appeared free from this problem (Additional file 2 includes detailed information). Interestingly, the region 1.60 M-1.78 M on chromosome 4 (Figure (Figure1)1) is coincident with the heterochromatic knob known to be hypervariable in At . The 22.98 M-23.46 M region of chromosome 3 (Figure (Figure2)2) looked like an At-homoeolog deletion. These results show that tiling arrays can be a useful tool for detecting copy number variation  and large-scale alterations in the As genome. As these analyses are based on non-normalized signals (between species), they are likely error-prone for individual genes.
To analyze the homoeolog-specific retention and expression of individual genes, we focused on 1,393,557 probes mapping to coding regions using Bowtie . Since Aa and At sequences differ at 1 out of 20 bases, some 25-base oligonucleotides designed for At are a perfect match for Aa sequences. Whenever orthologous Aa sequences mis-match to the At chip, this hybridization is weakened (hereafter termed 'diagnostic features' (DFs)). Separately for every gene, we identified a scaling factor based on probes with similar signatures of hybridization to normalize intensities between species. We then identified homoeolog-specific DFs and only retained those (405,466) robust over replicates (Figure (Figure3).3). We could only follow 24,344 genes as the fastest-evolving genes have too many DFs for normalization (Additional file 3).
We tested for deviations from an equal representation of the two homoeologs in the As genome [12,16,53]. As a reference point, we used the F1As DNA in which homoeologs are present at equal doses (Figure (Figure1).1). For each gene within the regions of putative alterations, we tested for changes in α between As and F1As, where α represents the relative contribution of Aa DF hybridization strengths in a hybrid genome. There was an upward shift in α in As compared to F1As (one-sided paired t-test, P < 2e-17), suggesting a preferential retention of homoeologs derived from the Aa parent (Figure (Figure4).4). Supporting this, more genes were called Aa-like (614) than At-like (324). This bias is significant, although moderate compared to earlier studies [30-32,34-36]. This might reflect a limited power of microarrays. For instance, we analyzed 30 genes encoded by the mitochondria organelle known to be At-derived. Only one plastid-encoded gene had enough DFs to be unambiguously classified, and was biased towards maternal At, as expected.
To identify homoeologous transcripts in As, we extracted RNA from leaf tissues and processed microarrays with the SNP-detection protocols similar to above. More than 49% of genes were called expressed, and 7,608 exhibited homoeolog-specific expression, with 3,458 and 4,150 exhibiting At-enriched and Aa-enriched DFs, respectively. Overall, we conclude that, over the 12,000 to 300,000 years, As has accumulated more deletions of At-originated homoeologs and uses the remaining At-originated homoeologs somewhat less (Table (Table2).2). Genes physically clustered together might co-express and co-evolve in transcript levels, as previously observed in flies . To test whether biases in homoeolog-specific expression were concordant between nearby genes, we calculated running averages of Δα along chromosomes (Figure (Figure5),5), and found regions with clusters of At-enriched and Aa-enriched transcription.
To validate the tiling array-based procedures above, we prepared Illumina libraries and performed RNA-sequencing of the As transcriptome. The Aa genome is not yet assembled, but we identified 52 Aa genes from GenBank and acquired an additional 50 genes from the UC Genome Center. We identified the orthologous At genes for these Aa genes and mapped the Illumina reads to both homologs. Nine genes did not contain any reads that were mapped to either homolog. For 14 genes, reads only mapped to either the Aa or the At reference. For the remaining genes, reads were aligned to both homologs and clustered as either derived from At or Aa (Figure (Figure6).6). We consider the number of uniquely mapped reads as a measure of homoeolog-specific expression. A strong correlation in Aa:At expression ratio between tiling arrays and the RNA-seq (R2 = 0.646, P < 5e-07) proves that both approaches work. This concordance is very satisfactory (Figure (Figure7)7) given that RNA samples were extracted from independently grown plants, and that microarray estimates are frequently noisy.
The summary of the Gene Ontology analysis of genes exhibiting homoeolog-specific retention and expression is shown in Tables Tables33 and and4.4. The categories 'cell communication' and 'signal transduction' were underrepresented, while 'DNA repair' and 'response to DNA damage stimulus' were overrepresented. Aa-enriched transcripts were overrepresented in the 'gene expression' category, including subprocesses involved in transcription, translation, RNA processing and gene silencing by miRNA.
Lastly, we considered homoeolog-specific expression in the context of At transcriptional networks . Of the 7,608 genes, connectedness estimates were available for 6,941 gene pairs. We tested whether bins of higher-connected gene pairs exhibited higher concordance of homoeolog-specific expression (Figure (Figure8).8). The fraction of concordant pairs was approximately 0.4 in low-connectedness bins, but increased to 0.8 for the high-connected gene pairs (R2 = 0.47, P < 0.0001). We also partitioned networks with homoeolog-specific expressions of at least two genes as co-biased for Aa (325), co-biased for At (219), or with mixed biases (302) (Table (Table5).5). The latter 'mixed' group was significantly underrepresented in comparison with random expectation (χ2 test, P < 6e-08).
In allopolyploid speciation, two genomes that have experienced long independent evolution are combined. Their genomes were shaped in different ways in response to the extrinsic environmental and intrinsic lifestyle pressures. We focused on As, a species that evolved 12 to 300 KYA from a single hybrid individual formed from an ovule of At and a pollen of Aa. Orthologous genes of At and Aa have average sequence divergence of 5% , exhibit differences in tissue-specific expression [10,24], and are located on five versus eight chromosomes. The allotetraploid hybrid initially had low fertility, if one can conclude this from the performance of artificial hybrids in the lab. This fertility can be restored through the complex interplay of genetic and epigenetic processes . Several groups have been fascinated with this rapid but complex process [10,22,24,45,46,53,56-59]. We focus on the subsequent longer-term molecular evolution, by comparing an evolved natural As with an 'unevolved' F1As hybrid.
F1As and its following generations are a model for whole-genome rearrangements and gene expression. Approximately one of ten cDNA amplified fragment length polymorphism (AFLP) bands displayed patterns that were non-additive between F1As and its parental species . One percent of bands were not detected in the parental species altogether . For AFLP fragments observed in the parents, homoeolog silencing was nearly symmetrical: 4% of At versus 5% of Aa. These patterns varied among tissues in a seemingly stochastic way. There was also some variation among accessions. In addition to AFLPs, Wang et al.  used spotted 70-mer oligonucleotide arrays to compare gene expression between At, Aa, and F1As. More than 15% of transcripts had different levels between parental species. In F1As, 5% of genes deviated in expression level from the additive mid-parent expectation, with the majority being repressed. Interestingly, 94% of these genes were more strongly expressed in the At parent, with their levels of expression in F1As resembling Aa [56,57]. In conclusion, the levels of gene expression in F1As more frequently resemble those in Aa, although homoeologs seem to have been used symmetrically and sometimes randomly. Aa-specific phenotypes, such as flower morphology, plant stature and long lifespan, are dominant in F1As (likewise, Arabidopsis lyrata phenotypes are dominant in thaliana-lyrata hybrids [56,59]). These results were confirmed and further detailed in very recent investigations [24,45,46].
We found that in As, Aa homoeologs are more frequently retained and more actively transcribed than their At counterparts. We hypothesize that these Aa-favoring biases are not random, but rather represent a signature of an evolutionary process. To explain these patterns, we propose a concept of 'homoeolog competition.' Genes are subject to detrimental mutations at approximately constant rates . Purifying selection removes these mutations with varying efficiencies depending on the gene redundancy, dominance, and other characteristics [6,21,60,61]. As some F1As homoeologs are functionally redundant, they should be progressively lost to mutations and deletions. From the initial pool of homoeologs, natural selection would preferentially maintain those with a higher contribution to fitness. In this sense, homoeologs 'compete'. Despite stoichiometric constraints to maintain stable ratios of dosage among genes , there is a well-documented shrinkage of polyploid genomes over time [6,9,12,15,18,21,25,26], as few genes are haploinsufficient .
Why would At-originated homoeologs be less valuable? Our first hypothesis is inspired by Hill and Robertson . Selfing organisms, such as At, are less capable of purging mildly deleterious mutations. This is because of severely reduced recombination in comparison to outcrossers, such as Aa [61,63,64]. This may seem paradoxical, as At maintains much less variation than Aa , which one might interpret as mutations in Aa. When selfing evolves, segregating mutations are quickly purged, as they exhibit their deleterious nature in autozygous individuals. In the short term, selfers are in fact better off . With time, however, Mullers' ratchet kicks in one slightly deleterious mutation after another, resulting in low standing variation but inferior functionality . Selfing is typical of terminal branches on phylogenetic trees, interpreted as being an evolutionary dead-end [64,65]. Thus, Aa homoeologs may contribute more to the fitness of an F1As, as they originate from an outcrossing species. In the future, we will test this hypothesis by population 'allele-specific' resequencing and applying molecular evolution tests to homoeologs separately.
Our second hypothesis involves historical factors. Suppose the southern-adapted At accession hybridized with the northern-adapted Aa accession, and that the emerging As accession spent most of the 12,000 to 300,000 years in the northern environment [37,39]. Aa-originated homoeologs would be a better fit for the environment, would be more frequently retained, and would evolve to be preferentially used . To test this hypothesis, one must sample As accessions from multiple locations, resequence their genomes and transcriptomes and identify environment-specific molecular evolution since the unique As speciation event. Our model assumes a large standing variation in the genome and transcriptome, which has been well-documented in Tragopogon [35,36]. A more direct, rather than biogeographic-type, evidence might be obtained with Gossypium . This species displays a similar strengthening of parentally skewed expression when natural allotetraploids are compared with F1 allotetraploid controls.
Thirdly, recall that the Aa transcription machinery is preferentially expressed in F1As . Homoeologs pre-adapted to function under Aa transcriptional control will then be selected for, reinforcing this initial pattern. Homoeolog-specific methylation might be at the heart of these processes [45,46]. Indirectly supporting this idea, Aa-like genes exhibited enrichment in the 'gene expression' category (with subprocesses: transcription, translation, RNA processing, and gene silencing by miRNA). Recent reports in Arabidopsis and Brassica allopolyploids indicate a high proportion of nonadditive expression for genes within these categories as well [53,67,68]. Similar results have also been shown in Senecio [69,70].
Imagine ancestral genes A1 and A2 that formed a functional dimer in the common ancestor of Aa and At 5 million years ago. These genes evolved into At1 and At2 orthologs in the At lineage, and into Aa1 and Aa2 orthologs in the Aa lineage. Within these lineages, At1 and At2 have been selected for the ability to form a dimer. Likewise, co-evolution has been taking place between Aa1 and Aa2 proteins . In F1As, along with the parental dimers At1-At2 and Aa1-Aa2, there will also be heterologous At1-Aa2 and Aa1-At2 dimers. Are these dimers likely to be functional ? Dobzhansky and Muller hypothesized that some would not be . Strongly decreased fitness of At × Aa F1 and F2 seeds, and meiotic disruptions in F1's, attest to the presence of intrinsic incompatibilities contributing to the reproductive isolation of these two species, and some genes involved have been characterized [61,62].
An allotetraploid might walk an evolutionary path to fitness restoration by preferentially co-expressing only one parental set of interacting homoeologs, with mixed networks being less common. The data confirmed our expectation that homoeologous networks in fact evolved towards pure Aa or At profiles. This type of 'D-M homoeolog conflict resolution' should be typical for polyploid ancestors and might potentially contribute to the fractionated genomes we observe today [9,72]. As we now know the identity of networks having evolved to a 'pure' parental type, our strong prediction is that the experimenter-induced heterologous state in these networks shall result in detectable reproductive losses.
When an allotetraploid is formed, the functions of homoeologs are partially redundant, and the genome is set for gene silencing and deletion. Thousands of genes affected by these processes in As were identified with tiling arrays and resequencing. These new computational approaches enable the use of widely available and economical tiling microarrays for the whole-genome analyses of species closely related to the sequenced references. In the As allotetraploid, more At-originated homoeologs are lost and silenced than Aa-originated homoeologs. We hypothesize that these Aa-favoring biases are not random, but rather represent a signature of an evolutionary process. Whenever more than one gene experiences silencing within a network, the homoeolog bias of the first event influences the likewise bias for the subsequent silencing; networks evolve towards their ancestral types. The mosaics of predominantly pure-parental networks in allotetraploids might contribute to phenotypic variability and plasticity, and enable the species to exploit a larger range of environments.
Affymetrix GeneChip® Arabidopsis Tiling 1.0R Arrays were hybridized with samples from four different sources. Genomic DNA was obtained from tetraploid At accession Ler , tetraploid Aa accession Care-1 , allotetraploid As accession Sue-1 , and an F1As produced by crossing the tetraploids At and Aa as maternal and paternal parents, respectively . cDNA was prepared from As leaf samples. All genomic DNA and cDNA samples were hybridized in three biological replicates using standard protocols.
RNA purification, cDNA synthesis and Illumina library construction was performed using the protocols of Mortazavi et al.  with the following modifications. Total RNA, mRNA, and DNA were quantified using a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA). mRNA fragmentation was performed using Fragmentation Reagent (Ambion, Austin, TX, USA) and subsequently cleaned through an RNA cleanup kit (Zymo Research, Irvine, CA, USA). Additional DNA and gel purification steps were conducted using Clean and Concentrator kits (Zymo Research). Illumina sequences are available for download at the NCBI Short Read Archive under the accession SRA025958.
The Arabidopsis Tiling Microarray is composed of over 3.2 million probe pairs tiled throughout the complete At genome. Probes are tiled at an average of 35 base pairs. Affymetrix CEL files are available for download from the public repository ArrayExpress under the accessions E-MEXP-2968 and E-MEXP-2969. To ensure that arrays within genotypes are comparable to each other, Robust Multiarray Analysis [75,76] was implemented to perform background correction. Intensities for three biological replicates were summarized with quantile normalization . In addition, intensities for the three biological replicates of As and F1As were summarized altogether with quantile normalization. Consistency and density plots may be found in the Additional files. PM probes exhibited some mismatches for the At genotype, as this array is based on a different reference; the arrays exhibited an additional lower hybridization intensity peak. PM probes from conserved exon regions were much more robust.
As expected from interspecific sequence divergence, the number of Aa higher intensity probes decreased, while the number of lower intensity probes increased. Note, however, that 'conservative features' and 'divergent features' peak at similar intensities in both species, making the analyses easier. Similar to At, Aa lower intensity probes were overrepresented in non-coding regions.
Probe intensities among three biological replicates in As were averaged and paired with the corresponding average among the three F1As replicates. For each gene, a paired Wilcoxon rank-sum test (FDR <0.05)  of all probes was used to identify genes with differential hybridization. The significance of individual genes might be misleading, but the pattern for multigene regions is robust. We scanned for windows in which at least 27 (90%) out of 30 genes exhibited unidirectional stronger or unidirectional weaker hybridization in As in comparison with F1As. We also required these differences to be significant at FDR <0.05 for at least 9 (30%) genes. Overlapping windows were collapsed to identify the entirety of these regions.
Our goal here is to select probe features enabling the comparison of At and Aa signal representation in As DNA and RNA. To enable cross-comparison of DNA and RNA, the analyses have to be made gene-by-gene, with DNA and RNA hybridization signals normalized to the same level with each gene.
First, probes representing conserved signatures between genotypes were identified and used to scale the entire gene. For every probe in a gene, its average intensity among replicates in At was compared to the average intensity in Aa. These ratios formed a unimodal distribution and the peak of this distribution was used as the scaling factor for which to normalize between genotypes for that gene. Mathematically, for probe i in the gene, the average intensity among j biological replicates in both genotypes is defined as:
where aij and tij represent the probe intensities of the jth replicate of the ith probe in Aa and At, respectively. Defining Xi as:
The scaling factor, xmax is defined as:
The value for xmax was estimated using the mlv function in R, which calculates the kernel density and searches for x that maximizes that estimated density function. From hereon, we replace all aij values with rescaled values represented by product(xmax,aij). We disregarded genes whose f(x) failed the Shapiro-Wilks normality test. This normalization method is similar to one recently outlined by Robinson and Oshlack , where a scaling parameter is used to normalize between two samples.
Second, we identified single feature polymorphisms or DFs between At and Aa using a Welch t-test of log2-transformed values, followed by controlling FDR to be smaller than 0.05. These approaches enabled us to analyze homoeolog-specific retention in 24,344 out of approximately 39,000 At genes.
If an As gene retained both parental homoeologs, we should observe an equal mix of At and Aa signals. A linear model was used to determine whether As has probe intensities within a gene contributed by i) both parents (mixed), ii) parental At only (At-like), or iii) parental Aa only (Aa-like). For a gene with n DFs, the vector of intensities, S = [S1, S2,..., Sn], may be contributed by corresponding Aa- and At-specific signals, such that S = α1•A + β1•T and the contribution of Aa, α1, can be estimated using a simple linear regression. Specifically:
where i = 1,2,..,n, j = 1,2,3 for the three biological replicates, and Ai and Ti are the mean intensities in Aa and At, respectively. εij are error terms that are independent random variables from a normal distribution with a mean 0 and variance σ2. The strength of our experimental design is in F1As, in which a null model holds true for genomic DNA. For F1As, this expectation is:
To detect deviations from the null, we tested whether α1 is significantly different from α2. Under the null hypothesis that α1 = α2, and assuming α + β = 1:
follows an F distribution with 1 and 6n - 1 degrees of freedom. This assumption of α + β = 1 can be made since the contributions of Aa and At are weighted. The bias was labeled as Aa-like if α1 > α2 and as At-like if α1 < α2. To account for multiple testing issues arising from thousands of genes tested, Benjamini-Hochberg's FDR was employed to adjust the significance level at 0.05 .
As with all linear regression models, we assume that the error terms follow a normal distribution. We investigated this by applying a Shapiro-Wilks test on each gene to ensure that they were normal. We removed over 7,000 genes that failed these tests. We found little discrepancy for the results of the analyses when α1 was defined as the At contribution. We also determined significance by performing a permutation test for each gene and found little discrepancy with the F distribution shown above.
Since we are estimating the relative contribution of Aa rather than the absolute, the expression level of every gene in the As transcriptome was normalized to identical hybridization levels with its corresponding genomic DNA. This was done using probes representing conserved signatures, identified as previously described. We then analyzed the homoeolog-specific expression with the same linear model approach as above, using DFs identified between RNA and DNA, and α found in As DNA as the null reference point. When these intensities of DFs are biased in one direction, we can determine homoeolog-specific expression. Furthermore, for each gene, α was estimated by regressing over all DFs in the set, minimizing spurious effects of individual probes. Forty-nine percent of genes were expressed. Distributions of intensities for conserved features in As DNA and RNA prior to and after gene-wise normalization are shown in the Additional files. The homoeolog-specific expression was assayed in 18,876 genes.
Pair-ended 72-base Illumina reads were aligned and mapped allowing up to 10 mismatches using bwa  to 102 Aa transcript sequences and their orthologous At sequences. A pairwise global alignment identified SNPs and short insertion/deletion variants between orthologous Aa and At gene pairs. Reads that mapped to either of the two orthologs were scanned for these variants to ensure that they were clustered with the appropriate ortholog (Figure (Figure6).6). The number of reads mapped to each ortholog was normalized to FPK (fragments per kilobase of exon) to account for slightly variable sequence length between orthologs. This analysis and its results are summarized in Figures Figures66 and and7,7, and in Additional file 5.
Note that although extant accessions of Aa, At, and F1As were used, As was formed 12 to 300 KYA, perhaps from different accessions. DFs and Illumina resequencing may potentially result in misleading conclusions. Nevertheless, 5 million years of sequence divergence between Aa and At compares favorably with the smaller amount of standing sequence variation and with the unaccounted extra divergence since As formation. From the above resequencing data, we estimated the divergence of the Aa homoeolog within As from the homologous gene in Aa. Likewise, we estimated the divergence of the At homoeolog within As from the homologous gene in At. Consistent with high sequence variation in Aa , the divergence from parental homologs is larger in Aa, as sequence variation in natural At is very limited . This would result in fewer Aa-like calls, and lower biases detected in this manuscript. Note that, as expected from , stronger expressed genes appear more conserved and exhibit lesser Aa and At divergences (Figure (Figure99).
Aa: Arabidopsis arenosa; AFLP: amplified fragment length polymorphism; As: Arabidopsis suecica; At: Arabidopsis thaliana; DF: diagnostic feature; F1As: F1 artificial allotetraploid; FDR: false discovery rate; KYA: thousand years ago; SNP: single-nucleotide polymorphism.
The authors declare that they have no competing interests.
PLC performed the computational and statistical analysis of the data, carried out the molecular resequencing, and drafted the manuscript. BPD performed sequence extraction and microarray experiments, and drafted the manuscript. MM performed sequence extraction and microarray experiments. LC participated in the design of the study. SVN conceived the study, participated in its design and coordination, and drafted the manuscript. All authors read and approved the final manuscript.
Differential hybridization between As and F1As. Excel table showing 6,790 genes with differential hybridization between As and F1As (Wilcoxon ranked sum test, FDR <0.05).
Differentially hybridized clusters between As and F1As. Excel table showing 1,643 genes found within differentially hybridized clusters between As and F1As. Clusters contain at least 30 genes with a strong unidirectional bias, where at least 27 genes have the same bias, and significant for at least 9 genes.
Outline of genes included and analyzed. Excel table outlining the number of genes discarded and included at each step in the analysis.
Homoeolog-specific retention in As DNA. Excel table showing 938 genes with homoeolog-specific retention in As DNA.
Comparison of homoeolog-specific expression estimated from At tiling microarray and Illumina resequencing. Excel table showing comparison of expression for 102 genes using both At tiling microarrays and Illumina resequencing.
Summary of probe hybridization intensities between At, Aa, As, and F1As. Probe hybridization intensities are shown for various regions throughout the genome (Figures S1 to S12). Density plots are shown for probe hybridization of DNA for PM and MM probes (Figures S13 to S16). A density plot is shown for conserved probes in As DNA and As RNA before and after gene-level normalization.
BPD and LC were supported by grant DBI0733857 from NSF Plant Genome Research Program. The authors are grateful to Joseph Fass, Meric Lieberman and Victor Missirian at the UC Genome Center for providing of A. arenosa sequences. The authors would also like to thank the anonymous reviewers for their comments and suggestions during the review of the manuscript.