LETTERS

Whole-genome resequencing reveals loci under selection during chicken domestication

Carl-Johan Rubin 1 *, Michael C. Zody 1,2 *, Jonas Eriksson 1 , Jennifer R. S. Meadows 1 , Ellen Sherwood 3 , Matthew T. Webster 1 , Lin Jiang 1 , Max Ingman 4 , Ted Sharpe 2 , Sojeong Ka [START_REF] Crawford | Poultry Breeding and Genetics[END_REF] , Finn Hallbo ¨o ¨k5 , Francois Besnier [START_REF]A genetic variation map for chicken with 2.8 million single-nucleotide polymorphisms[END_REF] , O ¨rjan Carlborg 6 , Bertrand Bed'hom [START_REF] Havenstein | Growth, livability, and feed conversion of 1957 versus 2001 broilers when fed representative 1957 and 2001 broiler diets[END_REF] , Miche `le Tixier-Boichard 7 , Per Jensen 8 , Paul Siegel 9 , Kerstin Lindblad-Toh 1,2 & Leif Andersson [START_REF] Darwin | On the Origins of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life[END_REF][START_REF]A genetic variation map for chicken with 2.8 million single-nucleotide polymorphisms[END_REF] Domestic animals are excellent models for genetic studies of phenotypic evolution [START_REF] Darwin | On the Origins of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life[END_REF][START_REF] Darwin | The Variation of Animals and Plants under Domestication[END_REF][START_REF] Andersson | Genetic dissection of phenotypic diversity in farm animals[END_REF] . They have evolved genetic adaptations to a new environment, the farm, and have been subjected to strong human-driven selection leading to remarkable phenotypic changes in morphology, physiology and behaviour. Identifying the genetic changes underlying these developments provides new insight into general mechanisms by which genetic variation shapes phenotypic diversity. Here we describe the use of massively parallel sequencing to identify selective sweeps of favourable alleles and candidate mutations that have had a prominent role in the domestication of chickens (Gallus gallus domesticus) and their subsequent specialization into broiler (meat-producing) and layer (egg-producing) chickens. We have generated 44.5-fold coverage of the chicken genome using pools of genomic DNA representing eight different populations of domestic chickens as well as red jungle fowl (Gallus gallus), the major wild ancestor [START_REF] Eriksson | Identification of the yellow skin gene reveals a hybrid origin of the domestic chicken[END_REF] . We report more than 7,000,000 single nucleotide polymorphisms, almost 1,300 deletions and a number of putative selective sweeps. One of the most striking selective sweeps found in all domestic chickens occurred at the locus for thyroid stimulating hormone receptor (TSHR), which has a pivotal role in metabolic regulation and photoperiod control of reproduction in vertebrates. Several of the selective sweeps detected in broilers overlapped genes associated with growth, appetite and metabolic regulation. We found little evidence that selection for loss-of-function mutations had a prominent role in chicken domestication, but we detected two deletions in coding sequences that we suggest are functionally important. This study has direct application to animal breeding and enhances the importance of the domestic chicken as a model organism for biomedical research.

For most of their history, domestic chicken populations have been bred for two purposes, egg laying and meat production [START_REF] Crawford | Poultry Breeding and Genetics[END_REF] . The effective chicken population size must have been huge in the past, before specialized commercial populations were established during the twentieth century, as a large proportion of farms kept a group of chickens interconnected with other groups by trade between regions, countries and continents. This is consistent with the extensive sequence diversity present in domestic chicken (,5 single nucleotide polymorphisms (SNPs) per kilobase (kb) in pairwise comparisons) [START_REF]A genetic variation map for chicken with 2.8 million single-nucleotide polymorphisms[END_REF] . During the twentieth century, specialized layer and broiler breeds were established to circumvent the inherent conflict in selecting for both growth traits (meat production) and reproductive traits (egg production) in the same bird. This approach, combined with the implementation of modern breeding methods based on quantitative genetics theory, has been extremely successful in improving productivity [START_REF] Havenstein | Growth, livability, and feed conversion of 1957 versus 2001 broilers when fed representative 1957 and 2001 broiler diets[END_REF] .

Our experimental design involved the resequencing of birds representing eight populations of domestic chickens as well as populations of red jungle fowl (Fig. 1). The aim was to identify the most common allele at the majority of polymorphic sites in the genome and to identify selective sweeps shared by populations selected for the same trait. We used the Applied Biosystems SOLiD technology to generate short reads (35 base pairs (bp)) to a sequence depth of four-to fivefold coverage from pools of chickens, each representing different populations, which is an effective approach for finding common alleles favoured by positive selection. The sequencing of DNA pools from individual lines is also effective for detecting selective sweeps, because heterozygosity can be calculated in sliding windows from sequences drawn from a pool of haplotypes. We sampled four different layer populations, four broiler populations and a pool of red jungle fowl representing two different zoo populations (Supplementary Table 1). Finally, we resequenced the particular red jungle fowl female (from the partly inbred UCD 001 line) previously used to generate the classical Sanger-sequenced reference chicken genome [START_REF] Chicken | Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution[END_REF] . The inclusion of the reference bird provided an excellent quality check of the SOLiD data.

We aligned the 35-bp reads to the chicken reference sequence. The uniquely placed reads covered 92% of the 1,043 megabases (Mb) in the current genome assembly (Supplementary Table 1), and these were used for all further analyses of sequence variation. The 90 Mb not covered by any read from any line were mostly repetitive sequences but to some extent probably represent sequences underrepresented after the emulsion-PCR step. We identified SNPs for each population individually and also for several pools of populations. Three independent reads of the same non-reference nucleotide were required to declare a position polymorphic. This criterion was chosen to minimize the number of false-positives. Bioinformatic analysis showed that the great majority of SNPs (.95%) reported in this study represent true SNPs (Supplementary Information), and this conclusion was supported by experimental verification of more than 300 SNPs (Supplementary Table 2). Furthermore, we eliminated ,40,000 putative SNPs because we did not find any support for the reference allele, and assume that these represent rare sequence errors in the current assembly (Supplementary Fig. 1). In total, 7,453,845 SNP loci of which we are confident remained after this filtering and were used in the subsequent analyses. Supplementary Table 1 shows the sequence coverage and the number of non-reference alleles found in each population. We used the allele frequencies of these SNPs to construct a genetic distance tree summarizing the genetic relationships among populations (Supplementary Information and Supplementary Fig. 2).

To detect putative selective sweeps, we searched the genome for regions with high degrees of fixation. The major challenge in such an analysis is to distinguish true sweeps from fixation due to genetic drift. In general, there is a high degree of genetic diversity in the chicken genome. However, some of the populations included were maintained with a limited effective population size, which could lead to genetic drift and fixation (Supplementary Table 1). To lessen the problem of fixation due to genetic drift within populations, we conducted sweep analyses with pooled sequence data from three groups of populations selected for similar purposes: all domestic lines, the two commercial broiler lines (CB-1 and CB-2) and three layer populations (WL-A, WL-B and RIR). We calculated the pooled heterozygosity, H p , in sliding 40-kb windows along the autosomes from sequence reads corresponding to the most and least frequently observed alleles at all SNP positions (Methods).

The distributions of observed H p values and the Z transformations of H p , ZH p , are plotted in Fig. 2a for the three comparisons. Owing to the complex and partly unknown demography of these populations, it is difficult to define a strict threshold that distinguishes true sweeps from regions of homozygosity caused by drift. We propose that the best way to evaluate the many putative sweeps is to cross-reference these results with others such as verification of the sweep in additional chicken populations and data on co-localization with major quantitative trait loci (QTLs) and differentially expressed genes. We focus our description of putative sweeps on those reaching a ZH p score of 26 or less, as these are in the extreme lower end of the distribution (Fig. 2a).

Only ,0.1% of the windows (n 5 58) had a ZH p score of less than or equal to 26 in the all-domestic comparison, and the corresponding fractions were ,0.05% for the broilers and ,0.01% for the layers (Fig. 2c and Supplementary Table 3). We used the yellow skin allele of BCDO2 locus as a proof of principle showing that this approach could reveal an established sweep 4 ; all domestic birds in this study were assumed to be homozygous for the yellow skin allele. We observed an H p value of 0.09 (ZH p 5 28.2) over BCDO2 (Fig. 2b), and the region perfectly overlapped the previously defined sweep [START_REF] Eriksson | Identification of the yellow skin gene reveals a hybrid origin of the domestic chicken[END_REF] . The region did not show complete fixation, as a few domestic birds carried a wild-type haplotype.

The number of putative sweeps reaching the ZH p -score threshold of 26 was highest for the all-domestic comparison: 58 windows representing 21 loci passed this threshold. Further genetic and functional studies of these loci are well justified because many of them, as well as some that did not reach the significance threshold, may have contributed significantly to chicken domestication. Three putative sweeps had ZH p scores that were more extreme (that is, lower) than that of yellow skin/BCDO2. One was located on chromosome 1 in a non-coding region upstream of SEMA3A, which encodes semaphorin 3A, an axon guidance molecule with an essential role in brain development [START_REF] Tran | Semaphorin regulation of cellular morphology[END_REF] . The second occurred in a non-coding region 160 kb upstream of the gene for V-set and transmembranedomain-containing protein 2A (VSTM2A), which is a predicted target-SNARE gene on chromosome 2, and the third occurred at the locus encoding thyroid stimulating hormone receptor (TSHR) on chromosome 5. We decided to further evaluate TSHR because it had the lowest ZH p score (29.2) and because of the well-established biological significance of TSHR signalling for metabolic regulation and reproduction [START_REF] Dobney | Genetics and animal domestication: new windows on an elusive process[END_REF][START_REF] Yoshimura | Light-induced hormone conversion of T 4 to T 3 regulates photoperiodic response of gonads in birds[END_REF][START_REF] Nakao | Thyrotrophin in the pars tuberalis triggers photoperiodic response[END_REF][START_REF] Hanon | Ancestral TSH mechanism signals summer in a photoperiodic mammal[END_REF] . The sweep region showed almost complete fixation over a 40-kb region (Fig. 3). The TSHR region also approached significance for the individual broiler and layer comparisons, with ZH p scores of 4.7 and 4.9, respectively. Both groups showed complete homozygosity at the TSHR locus but failed to reach the threshold of 26 because the sweep regions did not cover a full 40-kb window. In the all-domestic comparison, the consistent sweep signal from all eight populations resulted in an extremely low ZH p score even though the region of complete homozygosity across all populations did not cover an entire 40-kb window.

We decided to screen eight SNPs from the TSHR region using birds representing 36 populations with geographic origins ranging from Iceland to China (Supplementary Table 4). Every domestic chicken tested, representing commercial as well as local populations, carried at least one copy of the sweep haplotype; seven birds were heterozygous for the haplotype and 264 were homozygous. This is remarkable given the extensive genetic diversity present in the domestic chicken for most parts of the genome. Thus, TSHR may be a domestication locus in chicken, that is, a locus where essentially all individuals of a domesticated species carry a mutant allele. We observed the sweep haplotype at an intermediate frequency in red jungle fowl representing zoo populations (Supplementary Table 4), and believe that the most likely explanation for this is that many zoo populations have a history of some hybridization with domestic chicken.

We next searched for candidate mutations that may have been the target for the TSHR sweep and identified a non-conservative aminoacid substitution, namely glycine to arginine at residue 558 (Fig. 3c). This residue is located at the border between transmembrane region and the following extracellular loop (Fig. 3d). Glycine at this position is conserved among all known vertebrate TSHR sequences and it is almost completely conserved among other members of this family of glycoprotein-hormone receptors (see the Glycoprotein-hormone Receptors Information System; http://gris.ulb.ac.be/). A bioinformatic analysis using DASher [START_REF] Messina | DASher: a stand-alone protein sequence client for DAS, the Distributed Annotation System[END_REF] indicated that the glycine-to-arginine substitution pushes this residue outwards from the membrane and may therefore influence ligand interaction. Thus, this missense mutation is a candidate causal mutation for the TSHR sweep. Currently, we cannot formally exclude the possibility that the identified missense substitution is a slightly deleterious mutation that has hitch-hiked on a linked causative mutation. However, careful examination, including Sanger-based resequencing of the 40-kb sweep region to close any sequence gaps, did not reveal any strong suspects (GenBank accession numbers GU323554-GU323556).

A possibility is that the domestic TSHR allele confers an advantage as regards increased metabolic activity and growth. However, this is unlikely, because the TSHR locus does not coincide with any of the 13 growth loci detected in our previous QTL study of an intercross between White Leghorn and red jungle fowl [START_REF] Kerje | The twofold difference in adult size between the red junglefowl and white leghorn chickens is largely explained by a limited number of QTLs[END_REF] , where the two TSHR alleles are segregating. It is now well established that TSHR signalling between the pars tuberalis, of the pituitary gland, and ependymal cells in the hypothalamus regulates photoperiod control of reproduction in birds and mammals [START_REF] Yoshimura | Light-induced hormone conversion of T 4 to T 3 regulates photoperiodic response of gonads in birds[END_REF][START_REF] Nakao | Thyrotrophin in the pars tuberalis triggers photoperiodic response[END_REF][START_REF] Hanon | Ancestral TSH mechanism signals summer in a photoperiodic mammal[END_REF] . The TSHR sweep may therefore be related to a classical feature of domestic animals, namely the absence of the strict regulation of seasonal reproduction found in natural populations [START_REF] Price | Animal Domestication and Behavior[END_REF] . We plan to follow up the discovery of the TSHR sweep with functional receptor studies as well as careful phenotypic characterization of our White Leghorn/red jungle fowl intercrosses.

The screen for selective sweeps in populations of commercial broilers revealed many loci that make sense in relation to selection for muscle growth (Fig. 2c). A region on chromosome 1 with strong support contains both the gene for insulin-like growth factor 1 (IGF1), which is a candidate gene for growth, as well as the gene for pro-melanin-concentrating hormone (PMCH), which has important roles in appetite and metabolic regulation [START_REF] Shimada | Mice lacking melanin-concentrating hormone are hypophagic and lean[END_REF] . IGF1 also shows a strong indication of a selective sweep in the layer comparison (ZH p , 25.6; Supplementary Table 3) but not in the all-domestic comparison, as different haplotypes are close to fixation in broilers and layers. A functionally related candidate gene affecting growth traits and with a central role for insulin signalling, insulin receptor (INSR), also lies within a sweep region.

The sweep at the TBC1D1 (TBC1 (tre-2/USP6, BUB2, cdc16) domain family, member 1) locus is particularly interesting, because it is located at the major QTL explaining differences in growth between broilers and layers in three independent studies [START_REF] Sewalem | Mapping of quantitative trait loci for body weight at three, six, and nine weeks of age in a broiler layer cross[END_REF][START_REF] Zhou | Genome-wide linkage analysis to identify chromosomal regions affecting phenotypic traits in the chicken. I. Growth and average daily gain[END_REF][START_REF] Ambo | Quantitative trait loci for performance traits in a broiler 3 layer cross[END_REF] . Strong further support for the TBC1D1 sweep comes from the fact that the sweep haplotype is also fixed in the high and low growth lines (Fig. 1). This implies that the sweep took place during the early development of broiler chickens, as the high and low growth lines were established in 1957 and have been kept as closed populations since then. TBC1D1 has been associated with susceptibility to obesity in humans [START_REF] Stone | TBC1D1 is a candidate for a severe obesity geneand evidence for a gene/ gene interaction in obesity predisposition[END_REF] , and a loss-of-function mutation in this gene causes leanness in mice [START_REF] Chadt | Tbc1d1 mutation in lean mouse strain confers leanness and protects from diet-induced obesity[END_REF] . It has recently been shown that insulin-stimulated phosphorylation of TBC1D1 is required for the translocation of the main glucose transporter, GLUT4 (also known as SLC2A4), to the cell surface of mouse C2C12 myoblasts [START_REF] Peck | Insulin-stimulated phosphorylation of the Rab GTPase-activating protein TBC1D1 regulates GLUT4 translocation[END_REF] . This study implies that billions of broiler chickens around the world carry a mutant TBC1D1 haplotype and could be used for in-depth functional studies.

It has been proposed that loss of function may be an important factor in rapid evolution, such as occurs during domestication [START_REF] Olson | When less is more: gene loss as an engine of evolutionary change[END_REF] . Deleterious mutations may also accumulate as a result of relaxed selection or hitch-hiking during selective sweeps [START_REF] Lu | The accumulation of deleterious mutations in rice genomes: a hypothesis on the cost of domestication[END_REF] . We used the resequencing data to look for two specific loss-of-function mutations:

LETTERS

stop-codon mutations and deletions (Supplementary Information and Supplementary Tables 5 and6). We identified almost 1,300 deletions that were fixed or close to fixation in at least one population. Only 16 of these were longer than 6.4 kb and none was longer than 67 kb (Supplementary Fig. 3). We found little evidence that selection for loss-of-function mutations has had a prominent role in chicken domestication. However, we detected seven deletions in coding sequences (Supplementary Table 7), one of which has previously been described as functional and another of which shows strong evidence of function. The deletion in the growth hormone receptor (GHR) gene has previously been reported to be a causative mutation for sex-linked dwarfism [START_REF] Agarwal | Dysfunctional growth hormone receptor in a strain of sex-linked dwarf chicken: evidence for a mutation in the intracellular domain[END_REF] and has been used in some commercial broiler lines to reduce growth and feed consumption in parental lines. We also found a novel deletion that removes all but the first exon of the gene SH3RF2 (SH3 domain containing ring finger 2) (Fig. 4a). The deletion is fixed in the high growth line and occurs at a low frequency in the low growth line and in CB-1. SH3RF2 lies within a QTL region for body weight detected in a cross between the high and low growth lines [START_REF] Jacobsson | Many QTLs with minor additive effects are associated with a large difference in growth between two selection lines in chickens[END_REF] . We genotyped 400 birds from generation F 8 of this intercross, which provides a much better mapping resolution than generation F 2 (ref. 28). The analysis revealed a highly significant association between the presence of the deletion and increased growth (P , 0.001; Fig. 4b). Body weight at 70 days for Del/Del birds (600 g) was 20% greater than for WT/WT birds (500 g). We could take advantage of the fact that the deletion was also present at a low frequency in the low growth line. If this is the causative mutation, it should be associated with the same QTL effect in generation F 2 whether it was inherited from an F 0 chicken from the high growth line or the low growth line. This provides high resolution in QTL mapping because the lines were crossed after generation 41 of the selection experiment, meaning that they are separated by 82 generations of meioses, disrupting the associations between the QTL and linked markers. The estimated effects of the deletion haplotype originating from the high growth line were indistinguishable from those of it originating from the low growth line (Fig. 4c). This shows that the deletion or a mutation in the very near vicinity must be the causative mutation for the QTL. Expression analysis using hypothalamus messenger RNA revealed SH3RF2 expression in the low growth line but not in the high growth line, which is expected because the latter is fixed for the deletion (Fig. 4d). The result is of interest as it is well established that chickens from the high growth line have a genetic defect in hypothalamic appetite regulation [START_REF] Dunnington | Long-term divergent selection for eight-week body weight in white Plymouth rock chickens[END_REF] .

The results strongly suggest that the deletion, eliminating most of the coding sequence of SH3RF2, is the causative mutation for a QTL contributing to the remarkable difference in growth between the high and low growth selection lines. The nature of the mutation, that is, loss of function of a gene that is well conserved among birds and mammals, suggests that the deletion is unlikely to be just a neutral linked marker. SH3RF2 encodes SH3 domain containing ring finger 2 and is expressed in brain and muscle. The high and low growth lines now provide the first animal model with which the functional significance of this gene may be explored.

The present study casts light on the genetic basis of domestication, but also has implications both for the use of chicken as a model organism for biomedical research and for the application of genomics to practical chicken breeding (Supplementary Information). Chicken was selected for this study because the fairly small genome size (,1 Gb) made it easier to achieve sufficient sequence coverage to obtain a reasonable power to detect sweeps, deletions and high-quality SNPs. However, the successful outcome of this approach suggests that it should be applied to other domesticated species as well as to natural populations, where it may reveal the genetic basis for rapid evolutionary adaptations.

METHODS SUMMARY

SOLiD sequencing. We prepared libraries from pools of genomic DNA (except for the UCD 001 reference individual) and sequenced them using SOLiD (version 2 was used for all runs except two WL-A runs, for which we used version 1) according to standard manufacturer protocols. SNP detection and analysis. We called SNPs using the Corona Lite pipeline from Life Technologies. To call a SNP, at least three reads with different start sites supporting the non-reference allele had to be present. We called SNPs on each individual line and also by bioinformatically combining the data from multiple pools to improve sensitivity to rare alleles. All reported SNPs are unique sites from the union of all runs. Selective-sweep analysis. We used allele counts at SNP positions to identify signatures of selection in sliding 40-kb windows, for pools of sequence data. For each pool and SNP, we determined the numbers of reads corresponding to the most and least abundant allele (n MAJ and n MIN ). For each window in each breed pool, we calculated a pooled heterozygosity score: H p 5 2Sn MAJ Sn MIN / (Sn MAJ 1 Sn MIN ) [START_REF] Darwin | The Variation of Animals and Plants under Domestication[END_REF] , where Sn MAJ and Sn MIN are the sums of n MAJ and, respectively, n MIN for all SNPs in the window. Individual H p values were then Z-transformed as follows: ZH p 5 (H p 2 mH p )/sH p . Detection of deletions and stop codons. We first scanned the genome for putative deletions in each line. A putative deletion was defined as 100 bp or greater with no read coverage. We then tested all of these regions against coverage by the reference bird to remove those that were uncovered because of sequencing or alignment bias.

Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

METHODS

Animals. We included genomic DNA samples from the following chicken populations. Four layer lines were sampled: eleven males from White Leghorn line 13 (WL-A) developed at the Swedish University of Agricultural Sciences [START_REF] Liljedahl | Scandinavian selection and crossbreeding experiment with laying hens. I. Background and general outline[END_REF] and previously used for QTL mapping experiments [START_REF] Kerje | The twofold difference in adult size between the red junglefowl and white leghorn chickens is largely explained by a limited number of QTLs[END_REF] ; eight males from a commercial White Leghorn line (WL-B); ten males from the obese strain (OS) established in 1955 from a White Leghorn line as a model for autoimmune thyroiditis [START_REF] Cole | Hereditary hypothyroidism in the domestic fowl[END_REF] ; and eight males from a commercial Rhode Island Red population (RIR; Hubbard ISA) collected as part of the AvianDiv project [START_REF] Hillel | Biodiversity of 52 chicken populations assessed by microsatellite typing of DNA pools[END_REF] . Four broiler lines were sampled: ten males from a commercial broiler line (Ross 308, denoted CB-1 here), ten females from a second commercial line (CB-2), maintained by Hubbard ISA and collected by the AvianDiv project [START_REF] Hillel | Biodiversity of 52 chicken populations assessed by microsatellite typing of DNA pools[END_REF] ; and seven males and four females each from the high and low growth selection lines, both established from White Plymouth Rock chickens in 1957 [START_REF] Dunnington | Long-term divergent selection for eight-week body weight in white Plymouth rock chickens[END_REF] and used in previous QTL mapping experiments [START_REF] Jacobsson | Many QTLs with minor additive effects are associated with a large difference in growth between two selection lines in chickens[END_REF][START_REF] Carlborg | Epistasis and the release of genetic variation during long-term selection[END_REF] . Eight males representing two different zoo populations of red jungle fowl (RJ) were pooled. Finally, the single RJ female from the partly inbred UCD 001 line used to generate the chicken genome sequence 8 was also included. SOLiD sequencing. Libraries were prepared from pools of birds (except the UCD 001 reference individual) and sequenced using SOLiD, version 2 (except for two WL-A runs, for which we used version 1), according to standard manufacturer protocols. Reads were aligned in colour space to the reference assembly (version 2.1, Washington University) using the MAPREADS program allowing for up to three mismatches (including 'valid adjacent' changes as a single mismatch) and no indels. Only reads aligning uniquely in the genome were retained. SNP detection and analysis. SNPs were called using the Corona Lite pipeline from Life Technologies. To call a SNP, at least three reads with different start sites supporting the non-reference allele had to be present. We called SNPs on each individual line and also by bioinformatically combining the data from multiple pools to improve sensitivity to rare alleles. All reported SNPs are unique sites from the union of all runs. Following determination of all unique variant sites, we in silico genotyped all variant sites in all lines using custom PERL scripts, to compare all uniquely aligned reads in colour space to both the reference and variant alleles to count reference and variant alleles observed in each line.

Genetic differentiation was measured between lines with pairwise fixation index (F ST ), using an estimator introduced in ref. 34. This statistic was calculated across all SNPs identified in any line within the assembled chromosome sequences (that is, SNPs assigned to the '_random' and 'Un' sequences were excluded). Then a phylogenetic tree was constructed from these estimates using the neighbour-joining method implemented in PHYLIP [START_REF] Felsenstein | PHYLIP -Phylogeny inference package (version 3.2)[END_REF] . SNP validation. A VeraCode GoldenGate assay (Illumina) was designed targeting 384 SNPs that were called in the White Leghorn/red jungle fowl comparison. All SNPs were located in putative selective sweep regions identified in an initial screen for chromosomal regions showing a high degree of homozygosity. A total of 271 domestic chickens representing 36 populations and 51 red jungle fowl from a number of zoo populations were genotyped using the standard protocol provided by Illumina. The GoldenGate assays were read using a BeadXpress Reader (Illumina) and data were analysed on the GENOMESTUDIO V2009.1 software (Illumina). To estimate the false-negative rate of our SNP-calling pipeline, the sequence from the yellow skin allele 4 (GenBank EU334162) determined by classical Sanger-sequencing was used. By comparing the sequence against the reference genome sequence, a total of 175 SNP positions were detected. Those SNP positions were cross-referenced with the SNP positions called in the alldomestic comparison. Estimation of total SNP errors. All calls made at known SNP sites on the Z chromosome (ascertained in lines other than the reference bird) were reviewed and the number of valid base calls of each type made at that position in the reference bird was counted. For all sites with two or more calls, we assumed that any minority calls were errors (rather than using the reference base, to avoid counting errors at miscalled reference bases). As expected given that the Z chromosome in the reference female is hemizygous, we predominantly saw positions with two or three reads. We took the single-read, perposition error rate to be that seen at three reads of coverage (0.00644). We note that this includes both sequence detection errors and artefacts introduced by alignment errors and is approximately ten times worse than the manufacturer's specification of rates of accuracy after dibase correction (99.94%); independent assessment of single-base accuracy from invalid colour changes suggests that raw sequencing accuracy for these runs was very close to the manufacturer's specification. Because errors should be random, it was assumed that one-third of them would represent each incorrect base, and we thus computed for each depth of coverage in the genome the probability of seeing three or more of any given erroneous base. These probabilities were all multiplied by the total number of bases seen at each depth and summed over all depths. Selective-sweep analysis. Allele counts at identified SNP positions were used to identify signatures of selection in sliding 40-kb windows with a step size of 20 kb for three pools of sequence data (breed pools): all eight domestic lines, the two commercial broiler lines (CB-1 and CB-2) and three layer populations (WL-A, WL-B and RIR). At each detected SNP position, we counted the numbers of reads corresponding to the most and least frequently observed allele (n MAJ and n MIN , respectively) in each breed pool. This was done even if a particular SNP had not been originally detected by reads from the particular breed pool where the selective-sweep scan was performed. To prevent windows containing very few SNPs from adding spurious fixation signals, we omitted seven windows where only 1-10 SNPs had been detected and analysed 46,906 windows along chicken chromosomes 1-28 represented in the version 2.1 draft assembly. For each window, we calculated the pooled heterozygosity, H p , values using the formula H p 5 2Sn MAJ Sn MIN /(Sn MAJ 1 Sn MIN ) [START_REF] Darwin | The Variation of Animals and Plants under Domestication[END_REF] , where Sn MAJ and Sn MIN are breedpool-specific sums of n MAJ and, respectively, n MIN counted at all SNPs in the window. The distributions of H p values resembled normal distributions (Fig. 2a), and we therefore transformed the H p values into Z scores: ZH p 5 (H p 2 mH p )/ sH p . The Z transformation allowed us to place the three breed pools into the same framework because ZH p values indicate the number of standard deviations by which H p deviates from the mean. We applied a threshold of ZH p # 26 for putative selective sweeps because windows below this threshold represent the extreme lower end of the distribution (Fig. 2a). Calculation of identity scores. We calculated identity scores, IS, to visualize haplotype sharing in pairwise comparisons at the TSHR locus (Fig. 3b). For each identified SNP, we determined the fraction of reads that corresponded to the reference genome allele, F, in each chicken line. The IS values of individual SNPs were then calculated as IS 5 1 2 (jF line1 2 F line2 j), with SNPs assessed only if at least one read was obtained in each line. The IS value for a window was the mean of all SNP IS values observed in the window for a specific comparison. QTL analysis. The phenotypic effect of the SH3RF2 deletion was estimated in two ways. First, the genotype-phenotype map for the three genotypes (Del/Del, Del/WT and WT/WT) was estimated without accounting for the line origin of the alleles and accounting for the observed allele frequencies using the R package NOIA [START_REF] ´lvarez-Castro | A unified model for functional and statistical epistasis and its application in quantitative trait loci analysis[END_REF][START_REF] Le Rouzic | Estimation of genetic effects and genotypephenotype maps[END_REF] , version 0.92. Second, the effect of the deletion was estimated independently within the low growth and high growth lines by computing the IBD (relationship) matrix following ref. 38, under the assumption of three subsets of alleles in generation F 2 (Del originating in F 0 birds from the high growth or low growth lines and wild type (WT) always originating in F 0 birds from the low growth line). Detection of deletions and stop codons. We first scanned the genome for putative deletions in each line. A putative deletion was defined as 100 bp or greater with no read coverage. We then tested all of these regions against coverage by the reference bird to remove those that were uncovered because of sequencing or alignment bias. We computed the binomial p value of an N:0 split of reference and line reads given the chromosomal total reference:line ratio. To correct for testing ,1,000,000 putative deletions in all lines, we required a P value of ,10 28 to consider a region deleted. However, if a putative deleted region reached this significance threshold in one line, a P value of ,10 25 was considered sufficient to call the presence of an overlapping deletion in other lines. Then we merged by single linkage all overlapping deletion events in multiple lines, on the assumption that such events were likely to be identical. Finally, we merged runs of consecutive deletions that occurred in the same line(s) if the sum of deleted bases was at least 90% of the total span length; manual inspection revealed that such events were usually incorrectly split because of sequences that attracted spurious read placement, or possibly because the line was not 100% fixed for the event. Two such events overlapped genes and were successfully validated as single deletions, supporting this 'bridging'. Because of the absence of W sequences in male birds and the large amount of unplaced W sequences, the W chromosome and all random and unplaced chromosomes were excluded from this analysis. To prevent regions of unusually high coverage in the reference bird from causing falsepositives, we ignored reads in any line overlapping any region that had coverage in any line exceeding that expected given all of our sequencing (Poisson P # 0.05, Bonferroni corrected for sites observed across all lines).

To detect stop codons, the Ensembl gene set was filtered to include only models that made a complete protein product of $100 amino acids. We then computed the impact of the alternate allele of every detected SNP on its in-frame codon for each transcript overlapping it. If either the reference or the alternate codon was a stop, it was retained for manual review. During manual review, we discarded any putative stop variant if the transcript failed to show homology to any known protein or if the variant codon lay outside a region of conservation in the alignment and was not supported by chicken EST evidence. doi:10.1038/nature08832

Validation of loss-of-function mutations. Putative deletions or introduced stop codons were amplified using approximately 20 ng of genomic DNA from individuals taken from the appropriate pool in 10-ml reactions containing 10 mM Tris-HCL (pH 8.3), 50 mM KCl, 2.5 mM MgCl 2 , 0.2 mM dNTP, 1 mM each of forward and reverse primer and 0.5 U AmpliTaq Gold (Applied Biosystems). Validation of products was performed by direct sequencing using BigDye Terminator v3.1 chemistry (Applied Biosystems) and DNA sequences were aligned and edited using CODONCODE ALIGNER (version 3.0.2, CodonCode). The true break points of validated deletions that affected coding sequencing were used to design multiplex PCR assays. These reactions were amplified as above, but with 1 mM of each multiplex primer (m1, m2 and m3) and visualized on 2% agarose gels. Multiplex assays were also applied to a diversity panel (Supplementary Table 9) to assess the prevalence of the mutation in a global setting. All primers were designed in PRIMER3 39 (version 0.4.0) and are given in Supplementary Table 9, and new sequences have been deposited in GenBank (accession numbers GU323512-GU323553). SH3RF2 expression study. We dissected hypothalamus from four-day-old birds from both the high growth and low growth lines (three male and seven female chicks from each line, for a total of 20 samples). RNA was isolated with TRIzol (Invitrogen) and purified with RNeasy Mini Kit (Qiagen). The quality of the total RNA was checked with the Agilent 2100 Bioanalyser (Agilent Technologies) and the RNA samples were frozen at 270 uC until used. Complementary DNA synthesis from total RNA and two-step quantitative PCR were performed as described previously [START_REF] Ka | Proviral integrations and expression of endogenous Avian leucosis virus during long term selection for high and low body weight in two chicken lines[END_REF] . Briefly, quantitative PCR with reverse transcription was performed using the iQ SYBR Green Supermix (Bio-Rad Laboratories) in combination with the MyiQ Single-Colour Real-Time PCR Detection System (Bio-Rad Laboratories), according to the manufacturer's protocol. Primers were designed with PRIMER EXPRESS software (version 1.5, Applied Biosystems) and are listed in Supplementary Table 9. The collected data were analysed using the 2 {DDCT method 41 in comparison with the average of the four-day-old highgrowth-line male expression level, and all results were normalized to the chicken b-actin gene.

Figure 1 |

 1 Figure 1 | Chicken lines resequenced. Schematic tree based on breeding history (branch lengths not to scale). The minimum time of divergence, when known, is indicated. Pools of genomic DNA from the following populations were resequenced: red jungle fowl (n 5 8), commercial broiler 1 (CB-1; n 5 10), commercial broiler 2 (CB-2; n 5 10), high growth line (n 5 11), low growth line (n 5 11), Rhode Island Red (RIR; n 5 8), obese strain (n 5 10), White Leghorn line 13 (WL-A; n 5 11) and a commercial White Leghorn line (WL-B; n 5 8).The single red jungle fowl female from the partly inbred UCD 001 line used to generate the reference genome sequence 8 was also included.

Figure 2 |

 2 Figure 2 | Selective-sweep analysis of the chicken genome. a, Distributions of pooled heterozygosity values, H p , and corresponding Z transformations, ZH p , for all 40-kb windows (n 5 46,906). Bins of H p and ZH p are presented along the x axes. m, mean; s, standard deviation. b, Selective sweep at the yellow skin/BCDO2 locus 4 . Genes in the vicinity of the sweep are shown below the graph. * marks the sweep interval defined using sequence data from all pools of domestic chickens. ** marks the 23-kb interval previously associated with yellow skin 4 . c, The negative end of the ZH p distribution presented along chromosomes 1-28. The horizontal dashed lines indicate the threshold at ZH p 5 26. Genes residing within 20 kb of a window with ZH p # 26 are indicated by their gene names. The location of TSHR is also indicated in the plots for commercial broilers and layers, even though ZH p did not pass the threshold of 26. Ensembl genes are denoted ENS:XXXXX, where XXXXX are the last five digits of the Ensembl chicken gene annotation (ENSGALG:11 digits). For windows with no gene overlap, the presence of PhastCons evolutionary conserved regions (ECRs) and expressed sequence tags (ESTs) is indicated with a single asterisk. ** indicates a gene candidate close to a window with no genes overlapping 620 kb.

Figure 3 |

 3 Figure 3 | Details of TSHR mutation and sweep region. a, Pooled heterozygosities, H p , in the alldomestic data set, plotted for 40-kb windows on chromosome 5. The most highly fixed window (grey box) overlaps the TSHR gene. b, Degree of haplotype sharing in pairwise comparisons among populations. SNP frequencies in individual lines were used to calculate identity scores in 40-kb windows. Greyscale/purple boxes (left) indicate the comparison performed on that row according to the shade/colour assigned to each chicken line (left versus right column). Heat-map colours (right) indicate identity scores. LR, layers; BR, broilers; RJ, red jungle fowl pool; Ref, reference genome haplotype; Dom, domestic chicken; OS, obese strain; High/Low, high/low growth line. c, A TSHR missense mutation Gly558Arg is highly frequent in domestic chicken. Dashes indicate positions identical to the red jungle fowl sequence. Polar and positively charged amino acids are green and blue, respectively. d, Schematic structure of chicken TSHR based on the model provided by the Glycoprotein-hormone Receptors Information System (http://gris.ulb.ac.be/). In the magnification, the amino acids are coloured according to their properties: blue, basic; red, acidic; green, polar; grey, hydrophobic. The amino-acid substitution is indicated by the red arrowhead (amino acid 558 in UniProt Q3V5M1).

Received 22

 22 September 2009; accepted 8 January 2010. Published online 10 March 2010.

Macmillan Publishers Limited. All rights reserved ©2010

Macmillan Publishers Limited. All rights reserved

Acknowledgements Thanks are due to J. Altimiras, J. Dodgson, F. Imsland and S. Kerje for providing genomic DNA samples. We would also like to thank E. Sonnhammar for help with the bioinformatic analysis of the TSHR sequence, the Uppsala Genome Center for SOLiD sequencing and the Uppsala Multidisciplinary Center for Advanced Computational Science for help with data storage and analysis. This work was supported by grants from the Swedish Foundation for Strategic Research, the Knut and Alice Wallenberg Foundation, the Swedish National Infrastructure for Computing and The Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning. K.L.-T. and O ¨.C. are supported by European Young Investigator awards.

Author Contributions C.-J.R. was responsible for the selective-sweep analysis; M.C.Z. was responsible for detection of deletions and stop-codon mutations, and contributed to the selective-sweep analysis and SNP detection; J.E. was responsible for SNP verification and experimental validation of selective sweeps, J.R.S.M. was responsible for experimental validation of deletions and stop-codon mutations; E.S. and T.S. contributed to bioinformatic analysis; M.T.W. was responsible for the phylogenetic analysis; L.J. and M.I. contributed to SNP detection; S.K. and F.H. contributed to the experimental analysis of SH3RF2; F.B. and O ¨.C. performed the QTL analysis of SH3RF2; B.B., M.T.B., P.J. and P.S. provided biological samples of critical importance; K.L.-T. contributed to the experimental design and analysis; and L.A. designed and coordinated the study. L.A. wrote the paper together with all co-authors, who contributed sections and edited drafts to yield the final version of the manuscript.


The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla * Transfer RNA (tRNA) values were computed on exons. { Micro RNAs (miRNAs) are members of known conserved miRNA families. { Eudicotyledons are represented by P. trichocarpa and A. thaliana. 1 Magnoliophyta (most flowering plants) are represented by P. trichocarpa, A. thaliana and O. sativa.

The French-Italian Public Consortium for Grapevine Genome Characterization*

The analysis of the first plant genomes provided unexpected evidence for genome duplication events in species that had previously been considered as true diploids on the basis of their genetics 1-3 . These polyploidization events may have had important consequences in plant evolution, in particular for species radiation and adaptation and for the modulation of functional capacities [START_REF] Urgv | Evry Genomique Ve ´ge ´tale, 2 rue Gaston Cre ´mieux[END_REF][5][START_REF] Jaillon | Genome duplication in the teleost fish Tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF][7][START_REF] Maere | Modeling gene and genome duplications in eukaryotes[END_REF][START_REF] Blanc | Functional divergence of duplicated genes formed by polyploidy during Arabidopsis evolution[END_REF][START_REF] Seoighe | Genome duplication led to highly selective expansion of the Arabidopsis thaliana proteome[END_REF] . Here we report a high-quality draft of the genome sequence of grapevine (Vitis vinifera) obtained from a highly homozygous genotype. The draft sequence of the grapevine genome is the fourth one produced so far for flowering plants, the second for a woody species and the first for a fruit crop (cultivated for both fruit and beverage). Grapevine was selected because of its important place in the cultural heritage of humanity beginning during the Neolithic period 11 . Several large expansions of gene families with roles in aromatic features are observed. The grapevine genome has not undergone recent genome duplication, thus enabling the discovery of ancestral traits and features of the genetic organization of flowering plants. This analysis reveals the contribution of three ancestral genomes to the grapevine haploid content. This ancestral arrangement is common to many dicotyledonous plants but is absent from the genome of rice, which is a monocotyledon. Furthermore, we explain the chronology of previously described whole-genome duplication events in the evolution of flowering plants.

All grapevine varieties are highly heterozygous; preliminary data showed that there was as much as 13% sequence divergence between alleles, which would hinder reliable contig assembly when a wholegenome shotgun strategy was used for sequencing. Our consortium therefore selected the grapevine PN40024 genotype for sequencing. This line, originally derived from Pinot Noir, has been bred close to full homozygosity (estimated at about 93%) by successive selfings, permitting a high-quality whole-genome shotgun assembly.

A total of 6.2 million end-reads were produced by our consortium, representing an 8.4-fold coverage of the genome. Within the assembly, performed with Arachne 12 , 316 supercontigs represent putative allelic haplotypes that constitute 11.6 million bases (Mb). These values are in good fit with the 7% residual heterozygosity of PN40024 assessed by using genetic markers. When considering only one of the haplotypes in each heterozygous region, the assembly (Table 1a) consists of 19,577 contigs (N 50 5 65.9 kilobases (kb), where N 50 corresponds to the size of the shorter supercontig or contig in a subset representing half of the assembly size) and 3,514 supercontigs (N 50 5 2.07 Mb) totalling 487 Mb. This value is close to the 475 Mb previously reported for the grapevine genome size [START_REF] Lodhi | A molecular marker based linkage map of Vitis[END_REF] .

Using a set of 409 molecular markers from the reference grapevine map 14 , 69% of the assembled 487 Mb, arranged into 45 ultracontigs and 51 single supercontigs, were anchored along the 19 linkage groups. Thirty-seven ultracontigs and 22 single supercontigs were oriented, representing 61% of the genome assembly (Supplementary Tables 2 and3). This assembly has been annotated by using a combination of evidence. The major features of the genome annotation are presented in Table 1b. The 8.4-fold draft sequence of the grapevine genome contains a set of 30,434 protein-coding genes (an average of 372 codons and 5 exons per gene). This value is considerably lower than the 45,555 protein-coding genes reported for the poplar (Populus trichocarpa) genome, which has a similar size, at 485 Mb (ref. 1), and even lower than the 37,544 protein-coding genes identified in the 389 Mb of the rice genome 2 .

Three different approaches revealed that 41.4% (average value) of the grapevine genome is composed of repetitive/transposable elements (TEs), a slightly higher proportion than that identified in the rice genome, which has a somewhat smaller size 2 . The distribution of repeats and TEs along the chromosomes is quite uneven (see below). All classes and superfamilies of TEs are represented in the grapevine genome, with a large prevalence of class I elements over class II and helitrons (rolling-circle transposons) (Supplementary Table 7). An analysis of the distribution of the repetitive elements in the different fractions of the grapevine genome based on the current annotation shows that introns are quite rich in repeats and TEs (data not shown). In addition, 12.4% of the intron sequence contains transposons as determined using our set of manually annotated elements, most of which (75%) correspond to LINE (long interspersed element) retrotransposons, which therefore seem to have contributed specifically to the intron size observed in grapevine (Supplementary Table 8).

In eukaryotes with large genomes, the coding and repeated elements are distributed over the chromosomes and may be more or less interlaced, hence defining gene-poor and gene-rich regions. It has previously been noticed that the distribution of the genes along the chromosomes of rice and Arabidopsis thaliana is fairly homogeneous 2,3 . In contrast, we observe large regions that alternate between high and low gene density in V. vinifera (Supplementary Figs 2 and3). As expected, the density of TEs reflects a pattern substantially complementary to gene density. We observe a similar characteristic in the genome sequence of poplar, therefore indicating a dynamic for the invasion of TEs that is shared with the grapevine (Supplementary Fig. 3).

A striking feature of the grapevine proteome lies in the existence of large families related to wine characteristics, which have a higher gene copy number than in the other sequenced plants. Stilbene synthases (STSs) drive the synthesis of resveratrol, the grapevine phytoalexin that has been associated with the health benefits associated with moderate consumption of red wine 15,[START_REF] Baur | Therapeutic potential of resveratrol: the in vivo evidence[END_REF] . The family of genes encoding STSs has a noticeable expansion: 43 genes have been identified. Of these, 20 have previously been shown to be expressed after infection by Plasmopara viticola, thus confirming that they are likely to be functional. The terpene synthases (TPSs) drive the synthesis of terpenoids; these secondary metabolites are major components of resins, essential oils and aromas (their relative abundance is directly correlated with the aromatic features of wines [START_REF] Mateo | Monoterpenes in grape juice and wines[END_REF] ) and are involved in plant-environment interactions. In comparison with the 30-40 genes of this family in Arabidopsis, rice and poplar, the grapevine TPS family is more than twice as large, with 89 functional genes and 27 pseudogenes. Classification based on known plant homologues reveals that the subclass of putative monoterpene synthases represents only 15% of the Arabidopsis TPS family [START_REF] Aubourg | Genomic analysis of the terpenoid synthase (AtTPS) gene family of Arabidopsis thaliana[END_REF] whereas this subclass represents 40% of the grapevine TPS family. This result suggests a high diversification of grapevine monoterpene synthases that specifically produce C 10 terpenoids present in aroma (such as geraniol, linalool, cineole and a-terpineol). Furthermore, the grapevine genome annotation has also revealed genes encoding homologues to the two forms of geranyl diphosphate synthases (GPPSs), the enzymes that produce the substrate for monoterpene synthases: both the homodimeric GPPS and the heterodimeric form are present; the latter is present only in plants such as Mentha piperita and Clarkia breweri, which produce large quantities of monoterpenes [START_REF] Tholl | Formation of monoterpenes in Antirrhinum majus and Clarkia breweri flowers involves heterodimeric geranyl diphosphate synthases[END_REF] . Most of the STS and TPS genes occur as 20 clusters, including up to 33 paralogous genes located in a 680-kb stretch.

Because global duplication events seem to be a frequent event in plant evolution [START_REF] Adams | Polyploidy and genome evolution in plants[END_REF] , we searched the genome of V. vinifera for paralogous regions by using protein sequence similarity. Paralogous regions are defined as chromosome fragments in which homologous genes are present in clusters. Statistical analysis [START_REF] Simillion | The hidden duplication past of Arabidopsis thaliana[END_REF] of these clusters reveals that 94.5% have high probability of being paralogous (P , 10 24 ; Supplementary Table 11). Most Vitis gene regions have two different paralogous regions, which we have grouped together as triplets (Supplementary Fig. 5; coverage details in Supplementary Table 10). We conclude that the present-day grapevine haploid genome originated from the contribution of three ancestral genomes. It is yet to be demonstrated whether this content came from a true hexaploidization event or through successive genome duplications. The resulting plant had a diploid content that corresponds to the three full diploid contents of the three ancestors; it may therefore be described as a 'palaeo-hexaploid' organism. A number of rearrangements have affected the original three complements after the formation of the palaeo-hexaploid state. However, the gene order has been sufficiently conserved to permit the alignment of most regions with their two siblings.

We explored the time of formation of the palaeo-hexaploid arrangement by comparing grapevine gene regions with those of other completely sequenced plant genomes. If the palaeo-hexaploid complement is present in another species, it should result in a onefor-one pairing of gene regions between the two species considered. In contrast, if another species's genome evolved before palaeohexaploid formation, it should result in a one-to-three relationship between the other species and the grapevine genome. The available genome sequences were those of poplar 1 , Arabidopsis 3 and rice (Oryza sativa 2 ), of which poplar is considered to be most closely related to grapevine. All clusters constructed between the orthologues in the three comparisons have P , 10 24 (Table 1c). When the gene order in poplar is compared with that in grapevine, there are two clear distributions. First, the grapevine regions align with two poplar segments, as would be expected from a recent whole-genome duplication (WGD) in the poplar lineage 1 . Second, each of the three grapevine regions that form a homologous triplet recognizes different pairs of poplar segments (Fig. 1a and Supplementary Fig. 6). This shows that the palaeo-hexaploidy observed in grapevine was already present in its common ancestor with poplar.

Poplar belongs to the Eurosid I clade. The sister clade to Eurosid I is that of Eurosid II, which contains the model species Arabidopsis. Its gene order was compared with that in the grapevine genome. Two distributions appear: first, most grapevine regions correspond to four Arabidopsis segments (Supplementary Fig. 7); second, each component of a triplicated group in grapevine recognizes four different regions in Arabidopsis (Fig. 1b). This shows that the grapevine palaeo-hexaploidy was present in the common ancestor to Arabidopsis and grapevine, and therefore that it is a trait common to all Eurosids. This is confirmed by the homology level distribution between paralogues of the grapevine, indicating a lower conservation than between Vitis/Arabidopsis orthologues (Supplementary Fig. 4). The Eurosid group contains many economically important flowering plants such as legumes, cotton and Brassicaceae. Our present results establish these species as having a palaeo-hexaploid common ancestor. The grapevine/Arabidopsis comparison also reveals that the Arabidopsis lineage underwent two WGDs after its separation from the Eurosid I clade [START_REF] Simillion | The hidden duplication past of Arabidopsis thaliana[END_REF][START_REF] Bowers | Unravelling angiosperm genome evolution by phylogenetic analysis of chromosomal duplication events[END_REF][START_REF] Vision | The origins of genomic duplications in Arabidopsis[END_REF][START_REF] Blanc | A recent polyploidy superimposed on older large-scale duplications in the Arabidopsis genome[END_REF] . This contradicts some models based on more indirect evidence that placed the most ancient of these two duplications at the base of the Eurosid group, or even earlier [START_REF] Urgv | Evry Genomique Ve ´ge ´tale, 2 rue Gaston Cre ´mieux[END_REF][START_REF] Adams | Polyploidy and genome evolution in plants[END_REF][START_REF] Simillion | The hidden duplication past of Arabidopsis thaliana[END_REF][START_REF] Bowers | Unravelling angiosperm genome evolution by phylogenetic analysis of chromosomal duplication events[END_REF] . Some studies had also suggested a possible third duplication event in the distant past of the Arabidopsis lineage, potentially at the base of the angiosperm radiation. The controversy about this third event is now resolved by the Vitis genome comparisons: this event corresponds to the palaeo-hexaploidy formation that remains evident in the grapevine genome but has been difficult to characterize in Arabidopsis and poplar because of the more recent WGDs. In particular, the Arabidopsis genome lineage has undergone many rearrangements and chromosome fusions such that the ancestral gene order is particularly difficult to deduce from this species (Fig. 2).

Grapevines, like Arabidopsis and poplar, are dicotyledonous plants that diverged from monocotyledons about 130-240 Myr ago [START_REF] Wolfe | Date of the monocot-dicot divergence estimated from chloroplast DNA sequence data[END_REF][START_REF] Crane | The origin and early diversification of angiosperms[END_REF] .

Because rice is a monocotyledon, we assessed the presence or absence of palaeo-hexaploidy in its genome sequence. The observed pattern is the opposite of that seen for Arabidopsis and poplar: constituents of a grapevine triplet are generally orthologous to the same group of rice regions (Fig. 1c and Supplementary Fig. 11). Because rice and grapevine are phylogenetically distant, it is more difficult to detect relations of orthology across the two whole genomes: rearrangements, duplication and gene loss have affected the gene orders differently in the two lineages (Supplementary Fig. 10). Even with this limitation, we observed numerous cases of one-to-three relationships between rice and grapevine (Supplementary Figs 8, 9 and 11); 23% of orthologous blocks include the paralogous regions that originate from the grapevine palaeo-hexaploidy. For Arabidopsis, this number is as low as 1.4% (this difference is significant at 5%: x 2 5 8.9; Supplementary Table 12), despite the fact that the Arabidopsis genome has suffered many gene losses since its two WGDs. These gene losses would be expected to obscure the orthologous relations with the grapevine genome, but they are clearly insufficient to explain the high number of one-to-three relationships observed in the rice-grapevine comparison. The most probable explanation for this excess is that the rice ancestor did not exhibit the palaeo-hexaploidy observed in the grapevine, poplar and Arabidopsis. These findings are summarized in Fig. 3: the triplicated arrangement is apparent after the separation of the monocotyledons and dicotyledons and before the spread of the Eurosid clade. Future genome sequencing projects for other clades of dicotyledons, such as Solanaceae or basal eudicots, will help in situating the triplication event more precisely, and eventually in establishing its precise nature (hexaploidization or genome duplications at distant times).

Public access to the grapevine genome sequence will help in the identification of genes underlying the agricultural characteristics of this species, including domestication traits. A selective amplification of genes belonging to the metabolic pathways of terpenes and tannins has occurred in the grapevine genome, in contrast with other plant genomes. This suggests that it may become possible to trace the diversity of wine flavours down to the genome level. Grapevine is also a crop that is highly susceptible to a large diversity of pathogens including powdery mildew, oidium and Pierce disease. Other Vitis species such as V. riparia or V. cinerea, which are known to be resistant to several of these pathogens, are interfertile with V. vinifera and can be used for the introduction of resistance traits by advanced backcrosses [START_REF] Eshed | An introgression line population of Lycopersicon pennellii in the cultivated tomato enables the identification and fine mapping of yieldassociated QTL[END_REF] or by gene transfer. Access to the Vitis sequence and the exploitation of synteny will speed up this process of introgression of pathogen resistance traits. As a consequence of this, it is hoped that it will also prompt a strong decrease in pesticide use.

The high quality of the assembly, due mainly to the highly homozygous nature of the PN40024 line, enables the discovery of three ancestral genomes constituting the diploid content of grapevine. The Greek historian Thucydides wrote that Mediterranean people began to emerge from ignorance when they learnt to cultivate olives and grapes. This first characterization of the grapevine genome, with its indication of a palaeo-hexaploid ancestral genome for many dicotyledonous plants, addresses fundamental questions related to the origin and importance of this event in the history of flowering plants. Future work may help in correlating the differential fates of the three gene complements with phenotypic traits of dicotyledonous species.

METHODS SUMMARY

Gene annotation. Protein-coding genes were predicted by combining ab initio models, V. vinifera complementary DNA alignments, and alignments of proteins and genomic DNA from other species. The integration of the data was performed with GAZE [START_REF] Howe | GAZE: a generic framework for the integration of gene-prediction data by dynamic programming[END_REF] . Details are given in Supplementary Information. Paralogous and orthologous gene sets. Statistical testing of homologous regions was performed as described in ref. 21.

Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

METHODS

Genome sequencing. The V. vinifera PN40024 genome was sequenced with the use of a whole-genome shotgun strategy. All data were generated by paired-end sequencing of cloned inserts using Sanger technology on ABI3730xl sequencers. Supplementary Table 2 gives the number of reads obtained per library. Genome assembly and chromosome anchoring. All reads were assembled with Arachne [START_REF] Jaffe | Whole-genome sequence assembly for mammalian genomes: Arachne 2[END_REF] . We obtained 20,784 contigs that were linked into 3,830 supercontigs of more than 2 kb. The contig N 50 was 64 kb, and the supercontig N 50 was 1.9 Mb. The total supercontig size was 498 Mb, remarkably close to the expected size of 475 Mb. This indicates that the PN40024 has retained few heterozygous regions.

Remaining heterozygosity was assessed by aligning all supercontigs with each other. We first selected the supercontigs more than 30 kb in size that were covered over more than 40% of their length by another supercontig with more than 95% identity. After visual inspection of the alignments, we added to this list the supercontigs more than 10 kb in size that aligned at more than 40% of their length with supercontigs identified previously. All potential cases were then inspected visually to discard potential heterozygous regions (aligning relatively homogeneously across their complete length) and retained repeated regions (with more heterogeneous alignments). This treatment identified 11 Mb of potentially allelic supercontigs. We confirmed that in most cases their coverage was about half the average of the homozygous supercontigs. [START_REF] Adam-Blondon | Construction and characterization of BAC libraries from major grapevine cultivars[END_REF] were assembled into 1,763 contigs with FPC [START_REF] Soderlund | Contigs built with fingerprints, markers, and FPC V4.7[END_REF] , v. 8. In parallel, 1,981 markers were anchored on a subset of BAC clones [START_REF] Lamoureux | Anchoring of a large set of markers onto a BAC library for the development of a draft physical map of the grapevine genome[END_REF] , among which 388 markers mapped onto the genetic map, and 77,237 BAC end sequences were obtained [START_REF] Lamoureux | Anchoring of a large set of markers onto a BAC library for the development of a draft physical map of the grapevine genome[END_REF] . Blat 32 alignments (90% identity on 80% of the length, fewer than five hits) were performed with BAC end sequences on the 3,830 supercontigs of sequences with lengths over 2 kb. The results were then filtered with homemade Perl scripts to keep only the occurrences in which two paired ends were matching at a distance of less than 300 kb and with a consistent orientation. Two supercontigs were considered linked to each other if two BAC links could be found or one BAC link and a BAC contig link. A total number of 111 ultracontigs were constructed with this procedure. Genome annotation. Several resources were used to build V. vinifera gene models automatically with GAZE [START_REF] Howe | GAZE: a generic framework for the integration of gene-prediction data by dynamic programming[END_REF] . We used predictions of repetitive regions by repeatscout [START_REF] Price | De novo identification of repeat families in large genomes[END_REF] , conserved coding regions predicted by the exofish method [START_REF] Roest Crollius | Estimate of human gene number provided by genome-wide analysis using Tetraodon nigroviridis DNA sequence[END_REF][START_REF] Jaillon | Genome-wide analyses based on comparative genomics[END_REF] , genewise 36 alignments of proteins from Uniprot 37 , Geneid [START_REF] Parra | GeneID in Drosophila[END_REF] and Snap 39 ab initio gene predictions, and alignments of several cDNA resources (Supplementary Information).

A weight was assigned to each resource to further reflect its reliability and accuracy in predicting gene models. This weight acts as a multiplier for the score of each information source, before being processed by GAZE. When applied to the entire assembled sequence, GAZE predicted 30,434 gene models.

Paralogous and orthologous gene sets. We identified orthologous genes in six pairs of genomes from four species: A. thaliana, O. sativa, P. trichocarpa and V. vinifera. Each pair of predicted gene sets was aligned with the Smith-Waterman algorithm, and alignments with a score higher than 300 (BLOSUM62; gapo 5 10, gape 5 1) were retained. Two genes, A from genome GA and B from genome GB, were considered orthologues if B was the best match for gene A in GB and A was the best match for B in GA.

For each orthologous gene set with V. vinifera, clusters of orthologous genes were generated. A single linkage clustering with a euclidean distance was used to group genes. The distances were calculated with the gene index in each chromosome rather than the genomic position. The minimal distance between two orthologous genes was adapted in accordance with the selected genomes. Finally, we retained only clusters that were composed of at least six genes for Arabidopsis and O. sativa, and eight genes for P. trichocarpa (Supplementary Table 10).

To validate the clustering quality we used a method described previously [START_REF] Simillion | The hidden duplication past of Arabidopsis thaliana[END_REF] . For each cluster we computed the probability of finding this cluster in the gene homology matrix (Supplementary Table 11). This matrix was constructed from two compared chromosomes with genes numbered according to their position on each chromosome, with no reference to physical distances.

Paralogous genes were computed by comparing all-against-all of V. vinifera proteins by using blastp, and alignments with an expected value of less than 0.1 were retained and realigned with the Smith-Waterman algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF] . Two genes A and B were considered paralogues if B was the best match for gene A and A was the best match for B. Moreover, clusters of paralogous genes were constructed in the same fashion as orthologous clusters (Supplementary Table 10).

Figure 1 |

 1 Figure 1 | Comparison between three paralogous Vitis genomic regions and their orthologues in P. trichocarpa, A. thaliana and O. sativa. Orthologous gene pairs are joined with a different colour for each of the three paralogous grapevine chromosomes 6 (green), 8 (blue) and 13 (red). a, Orthologous regions in the poplar genome are different for each of the three Vitis chromosomes, showing that the triplication predates the poplar/Vitis separation. One Vitis region recognizes two poplar segments because of a WGD in the poplar lineage after the separation. b, Orthologous regions with Arabidopsis are different for each of the three Vitis chromosomes. This shows that the Arabidopsis/Vitis ancestor had the same palaeo-hexaploid content. One Vitis region corresponds to four Arabidopsis segments, indicating the presence of two WGDs in the Arabidopsis lineage after separation from the Vitis lineage. c, Orthologous regions in rice are the same for the three paralogous chromosomes. This indicates that the triplication was not present in the common ancestor of monocotyledons and dicotyledons. The presence in rice of different homologous blocks is due to global duplications in the rice lineage after divergence from dicotyledons.

Received 5

 5 April; accepted 7 August 2007. Published online 26 August 2007. 1. Tuskan, G. A. et al. The genome of black cottonwood, Populus trichocarpa (Torr. & Gray). Science 313, 1596-1604 (2006). 2. International Rice Genome Sequencing Project. The map-based sequence of the rice genome. Nature 436, 793-800 (2005). 3. Arabidopsis Genome Initiative. Analysis of the genome sequence of the flowering plant Arabidopsis thaliana. Nature 408, 796-815 (2000). 4. De Bodt, S., Maere, S. & Van de Peer, Y. Genome duplication and the origin of angiosperms. Trends Ecol. Evol. 20, 591-597 (2005). 5. Scannell, D. R., Byrne, K. P., Gordon, J. L., Wong, S. & Wolfe, K. H. Multiple rounds of speciation associated with reciprocal gene loss in polyploid yeasts. Nature 440, 341-345 (2006).

Figure 3 |Figure 2 |

 32 Figure3| Positions of the polyploidization events in the evolution of plants with a sequenced genome. Each star indicates a WGD (tetraploidization) event on that branch. The question mark indicates that ancient events are visible in the rice genome that would require other monocotyledon genome sequences to be resolved. The formation of the palaeo-hexaploid ancestral genome occurred after divergence from monocotyledons and before the radiation of the Eurosids.

Table 1 |

 1 Global statistics on the genome of Vitis vinifera

	(a) Assembly
		Status	Number	N 50 (kb)	Longest (kb)	Size (Mb)	Percentage of the
							assembly
	Contigs	All	19,577	65.9	557	467.5	-
	Supercontigs	All	3,514	2,065	12,675	487.1	1 0 0
		Anchored on chromosomes	191	3,189	12,675	335.6	6 8 .9
		Anchored on chromosomes	143	3,827	12,675	296.9	6 0 .9
		and oriented
	(b) Annotation
		Number	Median size (bp)	Total length (Mb)	Percentage of the genome	%GC
	Gene	30,434	3,399	225.6	4 6 .3	3 6 .2
	Exons CDS	149,351	130	33.6	6 .9	4 4 .5
	Introns CDS	118,917	213	178.6	3 6 .7	3 4 .7
	Intergenic	30,453	3,544	261.5	3 4 .7	3 3 .0
	tRNA*	600	73	0.04	NS	43.0
	miRNA{	164	103.5	0 .002	NS	35.9
	(c) Orthology
		Number of orthologous proteins	Mean identity (%)
	P. trichocarpa	12,996	72.7
	A. thaliana	11,404	65.5
	O. sativa	9,731	59.8
	Common to eudicotyledons{	10,547
	Common to Magnoliophyta1	8,121

  Only one supercontig of each allelic pair was therefore conserved in the final assembly, which consists of 3,514 supercontigs (N 50 5 2 Mb) containing 19,577 contigs (N 50 5 66 kb), totalling 487 Mb. If the haploid genome size of 475 Mb is considered correct, then our final assembly contains only about 12 Mb of remaining heterozygosity, or 2.6%.A set of 30,151 bacterial artificial chromosome (BAC) fingerprints of the BAC clones of a Cabernet-Sauvignon library

Acknowledgements The sequencing of the grapevine genome was launched and carried out after a scientific cooperation agreement between the Ministry of Agriculture in France and the Ministry of Agriculture in Italy, involving l'Institut National de la Recherche Agronomique (INRA), Consiglio per la Ricerca e Sperimentazione in Agricoltura (CRA) and Friuli Venezia Giulia Region. This work was financially supported by Consortium National de Recherche en Ge ´nomique, Agence Nationale de la Recherche, INRA, and by MiPAF (VIGNA-CRA), Friuli Innovazione, Universita `di Udine, Federazione BCC, Fondazione CRUP, Fondazione Carigo, Fondazione CRT, Vivai Cooperativi Rauscedo, Eurotech, Livio Felluga, Marco Felluga, Venica e Venica, Le Vigne di Zamo `(IGA). We thank S. Cure for correcting the manuscript; F. Ca ˆmara and R. Guigo for the calibration of the GeneID gene prediction software, and the Centre Informatique National de l'Enseignement Supe ´rieur for computing resources.

Supplementary Information is linked to the online version of the paper at www.nature.com/nature.

Author Information The final assembly and annotation are deposited in the EMBL/ Genbank/DDBJ databases under accession numbers CU459218-CU462737 (for all scaffolds) and CU462738-CU462772 (for chromosome reconstitutions and unanchored scaffolds). An annotation browser and further information on the project are available from http://www.genoscope.cns.fr/vitis, http:// www.vitisgenome.it/ and http://www.appliedgenomics.org/. Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to P.W. (pwincker@genoscope.cns.fr).

The French-Italian Public Consortium for Grapevine Genome Characterization Olivier Jaillon 1 *, Jean-Marc Aury 1 *, Benjamin Noel 1 , Alberto Policriti 2,3 , Christian Clepet 4 , Alberto Casagrande 2,5 , Nathalie Choisne 1,4 , Se ´bastien Aubourg [START_REF] Urgv | Evry Genomique Ve ´ge ´tale, 2 rue Gaston Cre ´mieux[END_REF] , Nicola Vitulo 6,15 , Claire Jubin 1 , Alessandro Vezzi 6,15 , Fabrice Legeai 7 , Philippe Hugueney 8 , Corinne Dasilva 1 , David Horner 9,15 , Erica Mica 9,15 , Delphine Jublot 4 , Julie Poulain 1 , Cle ´mence Bruye `re [START_REF] Urgv | Evry Genomique Ve ´ge ´tale, 2 rue Gaston Cre ´mieux[END_REF] , Alain Billault 1 , Be ´atrice Segurens 1 , Michel Gouyvenoux 1 , Edgardo Ugarte 1 , Federica Cattonaro 2 , Ve ´ronique Anthouard 1 , Virginie Vico 1 , Cristian Del Fabbro 2,3 , Michae ¨l Alaux 7 , Gabriele Di Gaspero 2,5 ,Vincent Dumas 8 , Nicoletta Felice 2,5 , Sophie Paillard 4 , Irena Juman 2,5 , Marco Moroldo 4 , Simone Scalabrin 2,3 , Aure ´lie Canaguier [START_REF] Urgv | Evry Genomique Ve ´ge ´tale, 2 rue Gaston Cre ´mieux[END_REF] , Isabelle Le Clainche 4 , Giorgio Malacrida 6,15 , Ele ´onore Durand 7 , Graziano Pesole [START_REF] Seoighe | Genome duplication led to highly selective expansion of the Arabidopsis thaliana proteome[END_REF]11,15