High-density linkage mapping and evolution of paralogs and orthologs in Salix and Populus

SNP markers and genotyping

396 SNPs were selected and genotyped in the S₁ and S₃ mapping populations. In a first panel of 384 SNPs conversion rate was 92% (354/384) and the overall sample call rate for the approved SNPs was 97.5%. The reproducibility was 100% according to duplicate analysis of 2.3% of the genotypes. In the second 12 SNP panel, the SNP conversion rate was 83% (10/12) and the overall sample call rate for the approved SNPs was 99.0%. The reproducibility was 100% according to duplicate analysis of 39.5% of the genotypes.

Linkage mapping

After excluding SNPs that were monomorphic or contained null alleles, segregation data for 307 SNPs in 463 individuals were used to construct linkage map S₁ (see Additional file 2). 75 markers were heterozygous only in the female parent, 158 only in the male parent and 74 were heterozygous in both. In addition, 45 microsatellite markers were included, of which 25 segregated in both parents, six were maternally informative and 14 were paternally informative. Finally, 287 AFLP markers were included, with genotypic data on 89 individuals. Consensus linkage maps were constructed using all 639 markers together. In the LOD Groupings step in Joinmap, markers were grouped at a LOD threshold of 5.0 and 25 linkage groups were identified of which three were duplets and three contained only four markers. The other 19 groups represent the 19 linkage groups in poplar (Figure 1 and Figure 2). 48 AFLP markers were unlinked and excluded from further analysis. For the markers within each group, map calculations were performed in Joinmap based on recombination frequencies, which is the step when linked markers are placed on linkage groups and ordered relative to each other. For five groups (II, V, XIV, XV and XIX) no consensus LGs could be calculated due to lack of markers segregating in both parents. These groups were split into the parental groups and the markers tested separately within each parental group. A number of markers remained unlinked after the map calculations due to insufficient linkage to the other markers. These markers were tested separately within groups to examine if any of them were linked to each other. A number of small groups were defined in this way and included in the map. Ten of the groups (eight duplets and two with four markers) contained only AFLPs and since they cannot be connected to the poplar genome they were excluded from further discussion. After inspection of LG VI, marker ph27 was removed because of inconsistency in the placement of this marker compared to S₃ and poplar.

In general, the S₁ map without the AFLP markers was in very good agreement with the map including the AFLP markers. However, in LG VIII the gene order differed slightly when constructing the map with or without the AFLP markers. Since the gene order in this linkage group without the AFLP markers was in better agreement with the gene order in the poplar genome as well as with the S₃ map, we present LG VIII without AFLP markers in Figure 1 (all other LGs in Figure 1 and Figure 2 include AFLP markers). Furthermore, II-24_sa on LG II was not mapped when no AFLP markers were included. 88 markers in map S₁ showed distorted segregation ratios (P < 0.00001). XI_21_sa and R_66_sa were removed because they were distorted and had a high number of missing genotypes. The consensus S₁ linkage map with 495 linked markers aligned to the poplar physical map is shown in Figure 1 and Figure 2. The map spanned 2477 cM with an average interval between the markers of 5.0 cM. The linkage groups A to F and the female linkage groups present on LGs II, V, XIV, XV and XIX, where no consensus maps were constructed are not included in the total map length.

214 SNPs genotyped in 282 individuals were used to construct the linkage map S_3. 94 SNPs segregated in the female parent only, 79 in the male parent only, and 41 in both parents. 41 microsatellite markers were also included, 21 of which segregated in both parents, 13 were maternally informative and seven were paternally informative (Additional file 2). Consensus linkage maps were constructed using all 255 markers. For map S₃, marker groupings at a LOD threshold of 4.0 identified 27 linkage groups and seven markers remained unlinked. Similarly to S₁, some markers that did not map to the other markers within linkage groups due to insufficient linkage were tested separately. Groups Ia, Ib, V, VII, VIII, IX, XII, XIV, XV, XVII and XVIII are represented by more than one linkage group in S₃, and they were placed on the map based on positional information in S₁ (Additional file 3). The most plausible reasons why we were unable to detect linkage between the small groups and other linkage groups is likely an effect of the small mapping population and limited number of markers, and as a consequence markers flanking the gaps can be too far apart, contain too little information, or are separated by a region of high recombination rate. Distorted segregation ratios (p < 0.00001) were found for 38 markers. R_24_sa and XVI_5_pIII form a separate duplet that cannot with certainty be connected to any LG. The same applies to the group containing SB945 and VI_8_sa_pIII.

There are no major rearrangements distinguishing the two willow maps, although there are a few apparent discrepancies in marker order on some LGs. Most differences involve closely spaced markers located on small groups or at the distal ends of groups, where the exact positioning of markers is uncertain. In several cases there was a swop in the positioning of two markers with tiny genetic distances, most likely an effect of the small size of the S₃ mapping population. The S₃ consensus map, excluding the small groups S₃A and S₃B and the groups S₃Ia-2 and S₃V-2 that could not be placed in groups Ia and V respectively, consists of 221 markers aligned to the S₁ map and is shown in Additional file 3. S₃ spanned 1793 cM with an average interval between the markers of 8.1 cM.

We could link nine previously unmapped scaffolds in poplar to linkage groups in the willow genome (scaffold 64 and 134 to LG Ib + XVI, scaffold 88 to LG V, scaffold 132 to LG VIII, scaffold 899 to LG XI, scaffold 142 to XIII, scaffold 40 to LG XIV and scaffold 82 and 122 to LG XV). We were also able to merge some groups that were unlinked in Hanley et al. (2006) (IIIa and IIIb, VIb and VIc, VIIb and VIIc, IXa and IXb, XIIIa and XIIIb, XVIIa and XVIIb, XIXa and XIXb).

Comparative mapping

We found two major interchromosomal rearrangements distinguishing the karyotype of willow from poplar, involving poplar LG I and XVI and willow LG Ia and Ib. Markers located on poplar LG I are in willow placed on two different LGs; willow LG Ia aligns to the first part of poplar LG I, from 0.3 Mb to 14.8 Mb (Figure 1) and willow LG Ib aligns to markers on the second part of poplar LG I, from 16.0 Mb to 35.3 Mb. In addition, markers located on poplar XVI are in the willow map linked to LG Ib. These data support either a fusion resulting in poplar LG I or a fission resulting in willow LG Ia and Ib sometime since the divergence of the two lineages. Similarly, either a fission in the poplar lineage forming poplar LG I and XVI or a fusion in the willow lineage forming willow LG Ib must have happened since the divergence of the two species (Figure 1).

Generally, gene order within syntenic groups is very well conserved in the willow - poplar comparison. However, the data indicate a few differences between the two genomes, five of which are inversions involving several markers each. Three of these are located on willow LG Ib (Figure 1) and correspond in size to ~1.6 - 5.4 Mb in the poplar genome sequence. In poplar, these regions correspond to parts of LG XVI and I. The other two supported inversions are located on LG XI and XVIII and involve roughly half of the total lengths of the linkage groups (6.9 Mb on LG XI and 6.3 Mb on LG XVIII) (Figure 2). For LG XI, the inversion is not evident in the S₃ map, which shows the same marker order as poplar, while the inversion on LG XVIII cannot be detected in the S₃ map because the markers involved are located on a separate linkage group so the gene order cannot be determined. Apart from these supported differences there are a number of suggested gene order differences between the willow and the poplar genomes, it is however uncertain whether or not there are true differences or errors in either the poplar genome assembly or the willow linkage maps. The assembly of the poplar genome is not complete and most likely the genome sequence as it is presented today contains numerous gaps.

Sequence analyses

For 115 gene fragments where three sequences were present (willow, poplar paralog 1 and poplar paralog 2) alignments were constructed and used to estimate pairwise K_i between orthologous loci in willow and poplar, paralogous loci in willow and poplar and paralogous loci in poplar (Figure 3). The willow sequence, the expected part of poplar paralog1 to be amplified and the total genomic sequence of poplar paralog 2 and default parameters were used to construct the alignments. A total of 48 kb genomic DNA sequence were used for estimations of K_i. As can be seen in the distributions of orthologous K_i (Figure 3), there are sequence pairs with high K_i values (> 1). This is likely a result of erroneous alignments when gaps were present in one of the sequences, or if parts of the overlapping sequences were missing, resulting in alignments of non-homologous nucleotides. It is however clear that the majority of orthologous sequence pairs have K_i between 0.1 and 0.2. Including all alignments, mean K_i between orthologs in willow and poplar was 0.17 ± 0.01 and the median K_i was 0.12. The distributions of K_i values for paralogous loci show a peak in the number of sequence pairs around 0.3 - 0.4. Mean K_i between paralogs in willow and poplar was 0.48 ± 0.02 and between the paralogs in poplar, 0.45 ± 0.02. The corresponding medians were 0.42 and 0.37. K_i between orthologs were significantly lower than K_i between paralogs in willow and poplar (Mann Whitney U test, W = 7900, p < 0.00001) and between paralogs in poplar (Mann Whitney U test, W = 8095, p < 0.00001). K_i between the two estimates of paralogs were not significantly different (Mann Whitney U test, W = 14022, p = 0.1).

Sequence analysis of orthologus and paralogous loci in willow and poplar

Pairwise substitution rates in introns, K_i, between orthologous and paralogous loci in willow and poplar were estimated in 115 gene fragments covering a total of 47966 nucleotides. K_i between willow and poplar orthologs were approximately three times lower (median: 0.12) than between paralogs (willow - poplar paralogs: 0.42, poplar paralogs: 0.37). The two latter estimates did not differ significantly, suggesting that the two lineages have similar evolutionary rates since their divergence. The large difference in substitution rates between orthologs and paralogs was unexpected as the genome duplication and the divergence of the two genera have been suggested to have occurred within a short time span some 60 - 65 mya [9, 32]. Our data clearly dispute this conclusion and instead suggest that the genome duplication occurred much earlier than the divergence of the two lineages.

According to interpretations of available fossil data, both Salix and Populus occur in middle Eocene sediments that are around 45 million years old [7, 8]. Using this time as a proxy for the separation of the two lineages and K_i = 0.12, the rate of substitution (r) was estimated to 1.28 × 10^-9 per site and year (T = K_i/2r). Assuming that this substitution rate has been constant since the Salicoid duplication yields an estimated date of the genome duplication some 150 mya. This date is obviously much older than previous estimates, and is not consistent with other data showing that the Salicoid duplication is not shared with Arabidopsis[9]. Available data suggest that the lineages leading to Arabidopsis (Eurosids II) and Salicaceae (Eurosides I) separated 100 - 120 mya [33]. Assuming instead that the duplication occurred soon after the divergence of Eurosids I and Eurosids II at 110 mya yields an estimate of the divergence between Salix and Populus of around 40 mya. Considering the uncertainty in the fossil record of Salicaceae [8, 34] such a date is perhaps not inconceivable. Dating of duplication and divergence events is notoriously difficult due to uncertainties both in the fossil data and substitution rates. Acknowledging these uncertainties, our data still suggests that the Salicoid genome duplication occurred considerably earlier than previously suggested. In a recent study, the poplar genome duplication was estimated to have occurred 45 mya, although it was acknowledged this was an underestimation of the timing due to the much slower substitution rates in poplar compared to other species [35].

Our data indicate substitution rates in Salicaceae between 1.28 × 10^-9 and 1.68 × 10^-9 per site and year depending on choice of calibration date, which is close to estimates for conifers (0.7 - 1.3 × 10^-9[36]), but considerably lower than for annual plants with a reported range of 5 - 33 × 10^-9[37–39]. These results are in line with the generally lower rates of substitutions estimated in perennials as opposed to annuals [40].

Plant material and DNA extraction

The mapping population S₁ with 463 F₁ offspring was created by crossing a diploid hybrid male clone ('Björn'; Salix viminalis × Salix schwerinii) to a diploid S. viminalis female clone (78183). The S. schwerinii parent to Björn originates from Siberia while the S. viminalis parent as well as 78183 originates from southern Sweden. Initially, the parental clones to S₁ were selected based on variation in phenology traits, but have also shown variation in rust resistance (Melampsora spp.) and different growth traits.

The mapping population S₃ was made by crossing the diploid S. viminalis male clone (81084) with the diploid S. viminalis female clone (78195) to produce 282 F₁ offspring. 78195 originates from Southwest of Sweden and 81084 originates from Southeast of Sweden. The parental clones to S₃ were selected based on variation in gall midge resistance.

Genomic DNA was extracted from frozen young leaves with the FastDNA Kit (MP Biomedicals) according to the protocol provided with the kit and DNA concentrations were determined with a Nanodrop spectrophotometer (Nanodrop Technologies).

SNP marker development and genotyping

SNP primers (Additional file 2) were designed to amplify fragments of genes at distances of approximately every 1 Mb on each linkage group (n = 19) of the poplar genome (The poplar genome website: http://genome.jgi-psf.org/Poptr1_1/Poptr1_1.home.html). In order to minimize amplification of multiple loci, genomic regions with little or no homology to other regions were selected. In effect, this means that single copy loci were favoured over those with multiple copies. Primers were positioned in exonic regions and designed to amplify across introns of approximately 500-1000 bp. Primers were also positioned to amplify parts of candidate genes for drought tolerance [41], rust [42–44] and insect resistance [45] (Additional file 1 and 2). 698 primer pairs were designed for the purpose of this study and PCR reactions were run with DNA from the parents of the mapping populations, of which 456 (65%) gave clean PCR products suitable for sequencing. 79 primer pairs were taken from Hanley et al. (2006) and tested in the parents as well, of which 70 resulted in clean PCR products. The PCR products were directly sequenced at Macrogen Inc. (Macrogen, Seoul, South Korea).

The primers were tested in each of the four parents in 15 μl volume PCR reactions containing 10 ng genomic DNA, 1 × PCR buffer II (Applied Biosystems), 2.5 mM MgCl₂ (Applied Biosystems), 0.2 mM dNTP mix (Fermentas), 0.5 μM of each primer (Invitrogen), and 0.5 U AmpliTaq Gold DNA polymerase (Applied Biosystems). The PCR reactions were run on a MyCycler thermal cycler (Bio-Rad Laboratories) with a PCR profile consisting of 10 min denaturation at 95°C followed by 35 cycles of 30 s denaturation at 95°C, 30 s annealing at 55°C and 1 min extension at 72°C with a final 10 min 72°C step. Amplification success was determined by agarose gel electrophoresis. For the first 96 primer pairs, a poplar sample (Populus trichocarpa) was also included to test if failure to amplify a product in willow was due to differences in nucleotide sequence between willow and poplar and not a result of unsuccessful primer design.

The willow PCR products were cleaned with 1 μl of a mixture of Exonuclease I (New England BioLabs) and Shrimp Alkaline Phosphatase (SAP) (Fermentas) for every 5 μl of PCR product before they were sent for sequencing of both strands with the forward and reverse PCR primers. Sequences were aligned using the DNA Baser V. 2 software. Contigs were then examined for the presence of informative SNPs for mapping. The SNPs were genotyped in the S₁ and S₃ mapping populations using the Golden Gate Assay [46] from Illumina (San Diego) at the SNP Technology Platform, Uppsala University http://www.genotyping.se. We applied a nomenclature for the SNP markers with roman letters referring to the poplar linkage group the primers were designed to amplify, followed with an arbitrary number and the letters 'sa'. In case markers were informative only in one of the mapping populations or different SNPs were used in S₁ and S₃ populations, the suffix I or III was added to the name.

Microsatellite and AFLP markers

PCR primers for microsatellite loci were taken from literature [47–51] and tested in the four parents (for details of the markers see Additional file 2). The PCR products were genotyped to check for levels of polymorphisms. The microsatellite loci selected for mapping were multiplexed and separated on ABI3730XL instruments at the Uppsala Genome Center. The majority of the microsatellite markers were PCR amplified by the Uppsala Genome Center, although some were amplified in-house. The in-house PCR reactions were run according to the protocol as described above for the SNP markers. The PCR reactions run at the Uppsala Genome Center were amplified in 10 μl volumes containing 10 ng genomic DNA, 1 × PCR buffer, 5% DMSO, 0.2 mM dNTP mix, 0.02 μM forward primer, 0.2 μM reverse primer, 0.2 μM M13-tail, 0.005 U Hot Star DNA polymerase and 1.25 μg BSA. The touchdown PCR profile consisted of 5 min denaturation at 95°C, six cycles of 30 s denaturation at 94°C, 30 s annealing at 58 - 53°C (-1°C per cycle) and 30 s extension at 72°C, seven cycles of 30 s denaturation at 94°C, 30 s annealing at 53 - 50°C (-0,5°C per cycle) and 30 s extension at 72°C, 30 cycles of denaturation at 94°C, 30 s annealing at 50°C and 30 s extension at 72°C and finally extension at 72°C for 7 min. Microsatellites were scored using the Peak Scanner software V 1.0 (Applied Biosystems). The genotypic data from an AFLP analysis was included from an earlier study on a subset of 89 individuals from the S₁ cross [21].

Map construction

Map construction was performed with the JoinMap 3.0 software [52]. Linkage maps were constructed using a LOD (logarithm of odds threshold) of 5.0 for map S₁ and a LOD of 4.0 for map S₃ to determine marker groupings. The LOD scores were selected based on the number of groups formed for each LOD score. The Kosambi mapping function was used for map construction with the following JoinMap settings: Rec = 0.4, LOD = 1.0, Jump = 5. The resulting linkage maps were drawn using the MapChart 2.1 software [53]. We did not use the ability of the software to force problematic markers on to the map and used the map resulting after the second round of map building. Mendelian segregations of the different markers were tested with chi-square analyses. We constructed S₁ maps both with and without the AFLP markers but present the map including AFLPs.

The locations of the SNP markers developed in this study (markers with prefix 'sa') on the physical map of the poplar genome were placed according to the positions of the genes where the primers were designed. Since we did not have access to the target positions in the poplar genome of the markers taken from Hanley et al. (2006), the positions of these markers on the poplar physical map were determined by BLAST searching the poplar genome with the willow gene fragments.

Database searches and Sequence analyses

BLAST searches against the Populus trichocarpa genome database on the NCBI web site were performed with the willow gene fragments that had been sequenced for the purpose of SNP identification in this study. BLAST searches were also performed with the poplar genomic sequences of the genes where the primers were positioned (paralog 1), hence the expected ortholog to the willow sequences. Usually we found two strong hits, of which we assumed the second best BLAST hit was the paralogous gene copy resulting from the Populus whole genome duplication (paralog 2). Since when designing the primers, genes or genomic regions with multiple copies were avoided as much as possible, a second BLAST hit was many times not present. When a second BLAST hit was present, the positions of the gene fragment of the willow sequence, the poplar paralog 1 and 2 in the poplar genome were compared, to determine if the expected orthologous gene fragment had been amplified in willow or, if the paralog had been amplified.

When three sequences were present for a gene-fragment; i.e. the willow sequence, the poplar paralog 1 and paralog 2, alignments were constructed using Clustal W [54] in the Alignment Explorer tool in MEGA4 using default parameters [55]. With the aim of comparing the number of substitutions in introns (K_i) between orthologs and paralogs, pairwise K_i was estimated in orthologous loci between willow and poplar, paralogous loci between willow and poplar and paralogous loci in poplar using MEGA4 with the maximum composite likelihood method. When the willow gene-fragment showed higher sequence homology to poplar paralog 2 and the location in the willow map correspond to the position of poplar paralog 2, willow and paralog 2 were regarded as orthologs. Frequency distributions of K_i values for each pairwise comparison are presented in Figure 3. K_i values were statistically compared between the three groups by non-parametric Mann-Whitney U tests in Mini-Tab 15. Medians and means (± standard error of the means) are presented.