PCR recombination artifacts

When PCR amplification is performed on groups of closely related sequences, such as the PR-10 gene family, accurate sequences are essential to distinguish between members. When we initially determined the relationships among the recovered sequences, most clades in the Bayesian consensus tree had low posterior probabilities. Several sequences occupied intermediate positions between well-defined clusters. This suggested the possibility of recombination. Recombination could have occurred in vivo through a crossing-over or gene conversion between prior existing genes, or in vitro during the PCR through strand-switching or re-annealing of incompletely amplified fragments.

Evidence supports the view that recombinant sequences were PCR artifacts. Recombination signals were abundant in the sequences obtained after the 30-cycle PCR (Figure 1a) and virtually absent when 22 cycles were employed (Figure 1b). Several GenBank sequences showed clear evidence for recombination too (Figure 1c). Putative recombinants from our experiments lacked or nearly lacked unique mutations and could be separated into two or three stretches that were identical to other sequences obtained. The GenBank recombinants shared close to 100% sequence identity to combinations of other accessions.

Several of the PR-10 sequences from B. pendula that were available from GenBank/EMBL/DDBJ also showed clear evidence for recombination (Figure 1c). However, it was not possible to do a similar check as mentioned above for the GenBank sequences, and we can only hypothesize on the presence of artifact sequences without such direct evidence. Given the regular occurrence of artifact sequences in our experiments, we maintained only those 40 GenBank sequences that were not under suspicion of recombination artifacts for further analysis.

Phylogenetic analysis: newly isolated sequences

We sequenced 404 individual clones in both directions (Table 1). Fourteen different sequences were identified in the diploid cultivar 'Schneverdinger Goldbirke', 7 of which were unique for this cultivar. Fifteen different sequences (of which 10 were unique) were identified in the diploid cultivar 'Tristis'. Approximately twice as many different sequences, namely 28, were identified in the tetraploid cultivar 'Long Trunk'; of these, 16 sequences were of B. pendula origin (8 unique), and 12 sequences were from another Betula species (all unique). In total, 32 different sequences were found.

The Open Reading Frame (ORF) of the sequences was highly conserved and the alignment was straightforward. All but one ORF contained 483 nucleotides, coding for a putative protein that is 160 amino acids long. One sequence from 'Schneverdinger Goldbirke' required the inclusion of an indel between base 388 and 389 of the consensus. This sequence was denoted as a pseudogene, since the indel introduced a stop codon at 7 codons downstream. It cannot be excluded that this pseudogene is expressed, since the stop codon was located near the 3' end. The intron position was identical in all sequences and located at codon 62, being inserted between the first and second nucleotide. Most alleles had 5' splicing sites of AG:GT, with the exception of one allele that had a GG:GT splicing site. The 3' splicing sites were AG:GC or AG:GA. This is in concordance with known motifs for plant introns. The introns were relatively AT-rich (55–65%).

We determined the relationships among the PR-10 sequences from B. pendula. For this, the 'Long Trunk' sequences that were designated to the other parental species were excluded. Excluding primer traces, 171 of the 452 aligned exon positions were variable, while 150 positions were phylogenetically informative. The consensus tree from the Bayesian analysis indicated several well-defined clusters (Figure 2). We implemented a cut-off level of 98% identity and allowed maximally two alleles per cultivar per gene to estimate the number of genes. In this way, we putatively identified ten genes in 'Schneverdinger Goldbirke', eleven in 'Tristis', and thirteen in 'Long Trunk'. Thirteen different genes were distinguished when the information was combined. Each gene was identified as such in at least two birch cultivars and was characterized by a distinct intron, in most cases of a different size between 84 and 152 bp (Figure 2). The similarity between different alleles of one gene ranged from 98.9 to 100% identity in the exons, which corresponds to 0–5 SNPs. Bet v 1.02A and Bet v 1.02B were not well distinguished in the coding sequences, but had distinguishable introns. The pseudogene PR-10.03B-p01 from 'Schneverdinger Goldbirke' was identical to the PR-10.03B allele from 'Long Trunk' except for its indel.

When all alleles from 'Long Trunk' were included, 182 variable positions were identified among the 452 aligned exon positions. 154 positions were phylogenetically informative. The topology of the consensus tree from the Bayesian analysis showed that seven specific 'Long Trunk' genes were clustered pair wise to the B. pendula genes and these are likely to represent orthologuous genes from the second Betula species (not shown). Given the high identity (up to 100%) to PR-10 sequences from B. ermanii, this species, or a close relative, is likely to represent the second parental species (unpublished data, Schenk et al.). Intron sequences of orthologuous genes mostly showed slight differences in length or base pair composition. Three genes were recovered only from the unknown parent species and five only from B. pendula, indicating that Betula species do not necessarily have the same (number of) genes.

Phylogenetic analysis: GenBank sequences versus newly obtained sequences

We classified the genes into five subfamilies (I–V) based on identities of the coding regions and the intron. The average identity between alleles within each subfamily was 95–100%. The GenBank sequences were in part derived from RNA extractions from specific tissues. We used this tissue information to predict the transcription of the PR-10 genes (Table 2). For this, alleles of a single gene are assumed to have the same mode and location of transcription. Subfamily I and II consist of respectively four and three genes and include alleles that are homologous to the pollen-expressed Bet v 1 allergens [2, 13]. For four genes within these subfamilies, we found alleles that were 100% identical to pollen mRNA-derived sequences previously deposited into GenBank. Alleles from two other genes were 99.8% identical (1 SNP difference) to pollen mRNA-derived sequences from GenBank, which we take to predict the location of transcription for these six genes with a very high level of confidence. For only one of the genes in subfamily II, there was no mRNA-derived homologue in the GenBank database, but the high homology (97.8–98.2%) to the other genes in this subfamily suggests that this gene will be expressed in pollen as well.

Subfamily III consists of four genes, two of which have alleles that are 99.5% identical to homologues of the previously described ypr10a and ypr10b sequences, which are transcribed in roots and leaves [15]. Transcription of the other two genes in subfamily III is unknown, as is the transcription of the single gene in subfamily IV. Subfamily V consists of a single gene. One of the recovered alleles was 100% identical to the previously described ypr10c sequence that is also expressed in roots and leaves [15]. Given the fact that all sequences in subfamily III-V are less than 90% homologous to pollen-derived mRNAs and that there is no evidence of expression in pollen, we consider these as non-allergens.

Nomenclature

All recovered B. pendula alleles were added to Table 2. This table also lists all known GenBank accessions of Bet v 1 sequences that were not under suspicion of artifacts. If sequences were previously named, we included a cross-reference to the original nomenclature (Bet v 1a-n) as well as to the nomenclature of the allergen nomenclature committee (Bet v 1.0101 to 1.3101) [19]. Gao et al. [18, 20] designed a nomenclature system for the Mal d 1 to Mal d 4 allergens that differentiates between alleles from different genes. We have adopted their system in Table 2. Protein sequences that have less than 95% identity are designated according to the current allergen nomenclature. A Latin letter is added to the iso-allergen name for those genes that have more than 95% identity. Two numerals are added for each allelic variant at the protein level. Two additional numerals indicate silent mutations. For example, Bet v 1.02C02.02 refers to a silent mutation in the second protein variant of the Bet v 1.02C gene, which shares more than 95% identity with Bet v 1.02A and B. The PR-10 sequences that were not allergens (not pollen-expressed) are named in a similar fashion, except that names started with PR-10 instead of Bet v 1. If future evidence would indicate that a particular PR-10 gene is expressed in pollen after all, the name can easily be modified by replacing the PR-10 tag with Bet v 1.

PR-10 proteins and allergenicity

A high similarity between proteins increases the chance that they share epitopes, while on the other hand, a single amino acid change may influence allergenicity drastically. The high homogeneity among the PR-10 genes of B. pendula was reflected by a higher allelic variation at the nucleic acid level compared to the protein level. Hence, the 45 different genomic sequences encoded 32 different putative isoforms (Figure 3). This is consistent with the relatively large number of synonymous mutations compared to the number of non-synonymous mutations (K_a/K_s ratios) [21]. The number of synonymous (K_s) and non-synonymous (K_a) substitutions per site were calculated from pair wise comparisons of all alleles from the three B. pendula varieties. The average value was 0.080 for K_a and 0.247 for K_s, resulting in an average K_a/K_s ratio of 0.33 (n = 45). Analysis of the occurrence of non-synonymous mutations per codon indicated two interesting regions. The region between codon Asn⁴² and Ile⁵⁶ lacked non-synonymous mutations and is characterized by a phosphate-binding loop with the sequence motive GxGGxGx (Figure 3). Relatively many amino-acid differences were present beyond codon 125, especially between isoforms of the genes PR-10.04/PR-10.05 and PR-10.01 to PR-10.03 (Figure 3).

Previous research has identified isoforms with varying IgE reactivity within mixtures of pollen. In fact, this is also true for pollen from a single tree. Ferreira et al. [14] made a distinction between isoforms with high, intermediate and low IgE-binding activity. The high IgE-binding isoform Bet v 1a (X15877) clustered with alleles of gene Bet v 1.01A and was 100% identical at the protein level to the Bet v 1.01A01 isoform. The intermediate IgE binding isoforms Bet v 1c and f (X77265, X77268) clustered with alleles of the genes Bet v 1.01C and Bet v 1.02C and differed by one amino acid from the Bet v 1.01C01 and the Bet v 1.02C02 isoforms. The low IgE binding isoform Bet v 1d (X77266) clustered with alleles of gene Bet v 1.01B, and was 100% identical at the protein level to the Bet v 1.01B01 isoform. Nine amino acids have been identified that affect the allergenicity of the Bet v 1 proteins on the B-cell level [12, 14]. These are marked in Figure 3. The tetraploid 'Long Trunk' contained several isoforms with unique amino acid substitutions due to its putative hybrid origin. This greatly enlarges the variation in putative IgE-reactivity among isoforms of this cultivar.

Two major T-cell binding epitopes have been identified, which are positioned between amino acids 112–123 and 142–156 [22]. Activation of T lymphocytes with the Bet v 1_142–156 epitope induced cytokine production (IL-4, IL-5). The major T-cell epitope Bet v 1_142–156 has an amino acid variation at position Thr¹⁴³. At this position the genes Bet v 1.02B and Bet v 1.02C code for Ala¹⁴³ instead. We identified an additional variation in the Bet v 1.01D gene from the tetraploid 'Long Trunk', namely the presence of Arg¹⁵⁰ instead of Ser¹⁵⁰. The T-cell epitope Bet v 1_112–123 is more variable and contains amino acid variations in the genes Bet v 1.01B, Bet v 1.01C, Bet v1.02C and Bet v 1.02D (Figure 3).

PCR artifacts

When PCR amplification is performed on groups of closely related sequences, such as the PR-10 gene family, accurate sequences are essential to distinguish between members. We used Pfu polymerase, which has proofreading functionality and reduces the number of base substitution error rates. However, previous research has shown that Pfu polymerase may generate more and more complex recombination artifacts than Taq polymerase [23, 24] through incomplete primer extension and re-annealing to a different template [25], or strand switching between different templates [26]. Reducing the number of cycles [24] was an efficient solution to lower the amount of artifacts.

We identified PCR recombination artifacts in several of our sequences using the computer program Phylpro [27]. The comparison of independent PCRs enabled sequence validation and exclusion of both recombination and base substitution artifacts. The high occurrence of recombination artifacts in our experiments (27–46% after 30 PCR cycles) is not uncommon. For example, Wang and Wang [28] report 32% recombination artifacts after 30 cycles of PCR amplification for 16S rRNA genes. The occurrence of recombinant sequences within a mixture of approximately 15 different sequences is expected to be high, because almost all recombinations are detected. Half of the recombinations would, for example, remain undetected if only two templates are present. The presence of 13 genes also increases the amount of template compared to a single copy gene.

Using the same analysis to detect recombination in GenBank/EMBL/DDJB accessions, we observed a recombination signal in 22 out of 62 (35%) accessions. These are not necessarily all PCR artifacts, as some recombinations may have occurred in vivo during evolution of the genes. However, we found only one true recombinant in our dataset and no evidence of past recombination events in the comparison between sequences from two different species (within the hybrid). This indicates that the occurrence of in vivo recombination is probably rare.

Base substitution error rates for mixtures of non-proofreading Taq and proofreading Pfu are approximately 5.6 × 10^-6 under optimal conditions [29]. This error rate is ~2–4 fold higher when only Taq is used [29, 30]. Without confirmation in independent PCRs these errors can not be excluded and it is very likely that Bet v 1 sequences with base substitution errors have been deposited into GenBank/EMBL/DDJB. Therefore, not all published isoforms will be clinically relevant. In addition, the clinical relevance of the isoforms will be influenced by their expression levels. If multiple allergen isoforms exist, there is a risk for selecting a recombinant isoform with low IgE-reactivity as a diagnostic tool, or even selecting isoforms that resulted from PCR artefacts. We therefore strongly suggest the use of primers that are highly specific for one gene [18] or, preferentially, the application of multiple independent PCRs to facilitate sequence validation in future sequencing work on allergens.

PR-10 and Bet v 1 genes

The PR-10 gene family of Betula pendula was shown to encompass at least thirteen genes. This is a conservative estimate since we used strict inclusion criteria. The distinction between genes is supported by the presence of a distinct intron. Each gene was identified as such in at least two birch cultivars. We attributed previously described GenBank sequences to these genes. Alleles from ten identified genes had previously been described [2, 10–13], while we identified three new genes. The genes are grouped into five subfamilies, based on sequence homologies in the ORF and intron. Differences in transcription coincide with the division between subfamilies. An organization of PR-10 genes into subfamilies was also reported for Malus domestica [18] and for Pinus monticola [5].

A striking feature of the PR-10 isoforms in B. pendula is their homogeneity, which may extend to other Fagales species, such as alder (Alnus glutinosa) and hazel (Corylus avellana). The intron has conserved 3' and 5' splicing sites and is always located at codon 62, as is e.g. reported for C. avellana and M. domestica [11]. High homogeneity may result from strong purifying selection or from concerted evolution. The presence of low K_a/K_s ratios among the isoforms suggests the occurrence of purifying selection. Evidence for concerted evolution is present in the overall gene tree of the PR-10 family [17]. Concerted evolution causes genes to evolve as a single unit with members exchanging genetic information through gene conversion and unequal crossing-over. Tandemly arranged genes may have high conversion rates [31], while this is a prerequisite for the occurrence of unequal crossing-over. Most PR-10 genes in apple map to two loci and are arranged in a duplicated cluster [18]. This organization may be a common feature for PR-10 genes. However, as pointed out by Nei and Rooney [32], the molecular mechanism of gene conversion is not well understood, and the model of birth-and-death evolution of genes may also explain the evolution of the PR-10 gene family. The presence of pseudogenes, although at a low frequency, is therefore of particular interest. An analysis of a species at an intermediate evolutionary distance, such as C. avellana, would be useful to clarify which mechanisms determine the evolution of PR-10 genes and to investigate a possible recent radiation of PR-10 genes.

The birch genome contains at least seven pollen-expressed genes that encode a mixture of Bet v 1 isoforms with varying IgE-reactivity. Swoboda et al. [13] found that pollen mixtures from multiple trees contain multiple Bet v 1 isoforms. We identified 14 different Bet v 1 isoforms in the three cultivars, nine of which are entirely new. The IgE-reactivity has been tested for several isoforms using recombinant proteins [14]. The allergenicity of the new isoforms can be examined in the future by expressing the isoforms as recombinant proteins and use these in a SPT or T-cell activation tests. Ferreira et al. [14] divided the Bet v 1 isoforms into three groups according to their IgE-reactivity and confirmed the division between high, moderate, or low IgE-reactivity in a Skin Prick Test (SPT). One high and one low IgE-reactive isoform from their analysis were 100% identical to isoforms that we have obtained from a single tree, while two intermediate IgE-reactive isoforms differed only by one amino acid from the alleles of two other identified genes. This strongly suggests that isoforms of different IgE-reactivity are in fact alleles encoded by different genes. Thus, each examined cultivar has the genetic background to express a mixture of isoforms with a high, moderate, and low IgE-reactivity. We plan to confirm this at the protein level in the near future.

Nomenclature

The nomenclature of Bet v 1 raises several issues. The first isoforms were termed Bet v 1a to Bet v 1n by Swoboda et al.[13], but these have subsequently been renamed and incorporated into the official database of the allergen nomenclature committee [19]. This database currently lists 37 allergen isoforms that have been termed Bet v 1.0101 to 1.3101. However, as can be seen from the list of known isoforms in Table 2, not all published isoforms have been added to this list, even though several of these isoforms were obtained from pollen mRNA. On the other hand, isoforms which have been recovered only from mRNA from roots and leaves are included as allergens (Bet v 1.1101 to Bet v 1.1301). Also, several of the described isoforms are highly suspicious as we observed clear recombination signals. As a result, the list is a random series of alleles that belong to different genes and has no biological basis. Similar problems were described for the Mal d 1 and Mal d 3 allergens by Gao et al. [18, 20]. To allow differentiation between alleles from different genes for the Bet v 1 alleles we have adopted their system.

Allergenicity of birch trees

The exact isoform composition of the three cultivars differed due to allelic variations. This may result in differences in allergenicity between cultivars. However, an exact copy of the most allergenic allele, Bet v 1a, was present in all three cultivars. Quantity measurements on expression indicate that Bet v 1a is the dominant isoform in pollen [13]. Given that the diversity of Bet v 1 isoforms within a single tree is larger than the diversity between the examined B. pendula cultivars, a characterization of Bet v 1-type isoforms should be done in other Betula species as well. Investigations towards sensitization and immunotherapy should anticipate that patients are exposed to a mixture of Bet v 1 isoforms of different IgE-reactivity, even if pollen originates from a single birch tree. Differences in allergenicity between birch trees may also result from variation in allergen content. Variation in allergen content has e.g. been shown for apple [33] and olive pollen [34].

Many Betula species and Betula hybrids have higher ploidy levels (tetraploid, hexaploid, and even octaploid) than B. pendula and are likely to contain increased numbers of allergen isoforms, as we found in the tetraploid cultivar 'Long Trunk'. For example, the tetraploid B. pubescens is dispersed throughout Europe, while other exotic birch species are increasingly introduced as cultivars, contributing to a larger allergen pool. However, interspecific Betula hybrids, which have a higher ploidy level, may also pose a potential source of hypoallergenic trees. Especially hybrids between less related species often display a reduced fertility, which may result in a reduced or aborted pollen production.

If breeding for hypoallergenic trees is implemented, approaches should take into account that the Bet v 1/PR-10 genes may be clustered and differences in allergenicity between clusters of genes may be used to guide breeding efforts. Clustering may be determined by mapping studies [18] and by screening and partial sequencing of a genomic library. In the genome sequence of Populus [35] we indeed can observe an organization of PR-10 genes into clusters. Other approaches to generate hypoallergenic trees may include the search for sterile or low pollen producing trees, or the application of RNA interference technology, which proved to be successful to silence the Mal d 1 allergens in apple [36] without phenotypic abnormalities. One issue that needs the be resolved is that although birch PR-10/Bet v 1 proteins have been suggested to act as plant steroid carriers [37], the exact conditions under which transcription is induced are still unknown. Given the abundance of Bet v 1 in birch pollen, silencing may affect pollen viability. However, for breeding of hypoallergenic cultivars that are propagated vegetatively this would be considered a welcome side-effect.

Plant material

The natural distribution range of B. pendula Roth (syn. B. verrucosa) covers almost the whole of Europe. Several B. pendula cultivars have been bred, including interspecific hybrids (also referred to as interspecies hybrids) between B. pendula and other Betula species. We collected young leaves from three B. pendula cultivars in the collection of PPO Boskoop (WUR, the Netherlands), namely 'Long Trunk', 'Schneverdinger Goldbirke', and 'Tristis'. Fresh leaf samples were sent to Plant Cytometry Services (Schijndel, The Netherlands) and screened by flow cytometry to estimate the ploidy level. Diploid (B. pendula) and tetraploid (B. pubescens) controls were included. The cultivars 'Schneverdinger Goldbirke' and 'Tristis' were diploid, while the cultivar 'Long Trunk' was tetraploid. The latter is likely to be an interspecific hybrid between B. pendula and a second, unknown, Betula species. The alleles recovered from 'Long Trunk' were either assigned to B. pendula based on sequence and intron similarity or were considered specific for this cultivar. The specific 'Long Trunk' alleles were analyzed separately. DNA was extracted using the DNeasy Plant Mini kit (Qiagen) according to the manufacturer's instructions.

PCR, cloning, and sequencing

PR-10 alleles were amplified from birch DNA with primers designed after two cDNA sequences (X15877, X77601). The primers were complementary to the regions around the start and stop codons (shown in bold); BpI-For: 5'-AATCTCTCAGGCCATCATG GGTG-3', BpI-Rev: 5'-TA GTTGTAGGCATCGGAGTGTGC-3', BpII-For: 5'-ATCTCAGGTGATCATCATG GGTG-3', and BpII-Rev: 5'-TA GTTGTAGGCATTTGGGTGTGC-3'.

PCR amplification with both primer pairs was performed with mixtures consisting of 2 μl dNTP (1 mM), 2 μl 10× Reaction buffer, 0.8 μl MgCl₂ (25 mM), 1.2 μl forward primer (10 pmol/μl), 1.2 μl reverse primer (10 pmol/μl), 0.11 μl of a 1:9 mixture of Pfu polymerase (Stratagene) and Taq polymerase (Goldstar)(5U/μl), and 20–80 ng template DNA. H₂O was added to obtain a total volume of 20 μl. PCR mixtures were subjected to the following conditions: initial heating step at 95°C for 15 minutes, denaturation at 94°C for 30 s, annealing at 50°C for 45 s, and extension at 72°C for 60s. A final extension step of 10 min at 72°C was added after 22–30 cycles. Given the observation of recombination among the recovered sequences, we subsequently varied the number of PCR cycles at intervals of 2 cycles. The minimum number of cycles was established by visual inspection of the amplification products on agarose gel at 22 for the BpI primer pair and at 24 for the BpII primer pair. Originally, 30 PCR cycles had been used. We repeated the experiment at 22–24 cycles to ensure that amplification was in its linear phase.

To obtain the A-tailing that facilitates the ligation procedure, five additional cycles were run on 1–4 μl of PCR product with Taq polymerase (Goldstar). PCR conditions were similar as described above. PCR products were purified with the MinElute PCR Purification Kit (Qiagen). Purified samples were ligated into the pGEM-T easy Vector (Promega) and established in Escherichia coli XL1 Blue competent cells (Stratagene) according to the manufacturer's instructions. White colonies were picked from agar plates and grown overnight at 37°C in freeze medium. PCR-based screening was performed with vector-specific M13 primers. PCR products were purified with Sephadex G-50 (Millipore). The DYEnamic™ ET Terminator Cycle Sequencing Kit (Amersham) was used for the sequence reaction. Sequence products were analyzed on a 96-capillary system (ABI 3730xl).

Genomic B. pendula sequences have been submitted to GenBank as DQ296566–DQ296598 and DQ325525–DQ325535, and the specific 'Long Trunk' sequences as DQ296599–DQ296610.

GenBank sequences

Sixty-six PR-10 sequences were obtained from GenBank/EMBL/DDJB by searching with MegaBLAST for entries from B. pendula with more than 60% sequence identity to Bet v 1a (X15877). The search results included nineteen genomic DNA (gDNA) sequences: Z72429–Z72438[11]; AJ001551–AJ001557 (Cvitanich and Larsen, direct submission); AJ289770 and AJ289771 (Pellinen et al, direct submission). The accessions Z72435-8 were singletons with a relatively low homology to our and other GenBank Betula sequences and were excluded from the analysis. Forty-seven cDNA sequences were found: X15877[2]; X77200, X77265–X77274, X81972, X82028[13]; X77599–X77601[15]; Z80098–Z80106 (Larsen, direct submission); AF124837–AF124839[12]; AJ002106–AJ002110, and AJ006903–AJ006915[10].

Phylogenetic analysis

Nucleotide sequences were aligned using CLUSTALW [38]. Bayesian phylogenetic analysis was performed with MrBayes 3.1.1 [39]. The maximum likelihood model employed 6 substitution types, with base frequencies set to the empirically observed values. Rate variation across sites was modeled using a gamma distribution. The Markov chain Monte Carlo search was run twice with 4 chains for 1,000,000 generations. Topology and model parameters were sampled every 100th generation and used to estimate model parameters and to determine the posterior probabilities of clades. The first 100,000 generations were discarded as "burn in". The outgroup was composed of PR-10 sequences from Prunus armeniaca (AF020784), P. avium (U66076), Pyrus communis (AF057030), and Malus domestica (X83672, Z72425, Z72427). To date, these are the closest non-Fagales relatives of the PR-10 family from B. pendula. We confirmed the phylogenetic analysis by constructing a neighbour-joining tree with Kimura two-parameter distances. Bootstrapping was carried out with 1,000 replicates. The results were similar for both analyses; therefore, the results from the neighbor-joining analysis are not shown.

Based on variation among Mal d 1 alleles of a single locus compared to alleles of different loci [18], we predefined a cut-off level of 98% identity in order to identify clusters that encompass alleles of the same gene. A limitation of the applied method is that genes that are homozygous and differ less than 98% from each other remain undetected.

Recombination was detected and visualized with Phylpro 1.0 [27]. This program computes the correlation coefficient of pair wise distances between the target sequence and all other sequences on both sides of a sequence position. We used a sliding window of 40 base pairs.

Isoforms and allergenicity

Nucleotide sequences were aligned codon-by-codon. We analyzed general selection patterns at the molecular level using DnaSp 4.00 [21]. The number of synonymous (K_s) and non-synonymous substitutions (K_a) per site were calculated from pair wise comparisons with incorporation of the Jukes-Cantor correction. Nucleotide data were translated into amino acids. Putative isoforms were analyzed for their potential allergenicity by screening those amino acid positions that have been identified as influencing the IgE binding [12, 14] and T-cell activation [22].