Assortative mating and differential male mating success in an ash hybrid zone population

Abstract

Background

The structure and evolution of hybrid zones depend mainly on the relative importance of dispersal and local adaptation, and on the strength of assortative mating. Here, we study the influence of dispersal, temporal isolation, variability in phenotypic traits and parasite attacks on the male mating success of two parental species and hybrids by real-time pollen flow analysis. We focus on a hybrid zone population between the two closely related ash species Fraxinus excelsior L. (common ash) and F. angustifolia Vahl (narrow-leaved ash), which is composed of individuals of the two species and several hybrid types. This population is structured by flowering time: the F. excelsior individuals flower later than the F. angustifolia individuals, and the hybrid types flower in-between. Hybrids are scattered throughout the population, suggesting favorable conditions for their local adaptation. We estimate jointly the best-fitting dispersal kernel, the differences in male fecundity due to variation in phenotypic traits and level of parasite attack, and the strength of assortative mating due to differences in flowering phenology. In addition, we assess the effect of accounting for genotyping error on these estimations.

Results

We detected a very high pollen immigration rate and a fat-tailed dispersal kernel, counter-balanced by slight phenological assortative mating and short-distance pollen dispersal. Early intermediate flowering hybrids, which had the highest male mating success, showed optimal sex allocation and increased selfing rates. We detected asymmetry of gene flow, with early flowering trees participating more as pollen donors than late flowering trees.

Conclusion

This study provides striking evidence that long-distance gene flow alone is not sufficient to counter-act the effects of assortative mating and selfing. Phenological assortative mating and short-distance dispersal can create temporal and spatial structuring that appears to maintain this hybrid population. The asymmetry of gene flow, with higher fertility and increased selfing, can potentially confer a selective advantage to early flowering hybrids in the zone. In the event of climate change, hybridization may provide a means for F. angustifolia to further extend its range at the expense of F. excelsior.

Background

Hybrid zones, where lineages differentiable by one or more heritable traits meet and intercross, provide unique opportunities for studying the nature and dynamics of barriers to gene exchange. The evolution of these barriers can have many different outcomes, including divergence of populations leading to speciation, collapse of differentiated populations, hybrid speciation or invasion. The structure and evolution of hybrid zones depend mainly on the relative importance of dispersal, local adaptation and the fitness of hybrids [1, 2], influencing the strength of reproductive isolation. For example, with relatively high dispersal between adjacent populations, gene exchange between species is prevented only if local adaptation is sufficiently strong to eliminate hybrids. Temporal isolation is a particular ecological isolation process that can result from divergent adaptations and cause assortative mating by itself. It involves differences in reproductive timing and can be the result of behavioral or developmental schedule divergences [3, 4]. In this study we focus mainly on the role of dispersal and temporal assortative mating in shaping the mating patterns in a plant hybrid zone population, and on the relative male fitness of hybrids and parental species.

Pollen dispersal is an important factor influencing the dynamics and evolution of plant populations (e.g. [5, 6]). In particular, the frequency of long-distance dispersal events can have strong effects on the distribution of genetic diversity by connecting distant demes in metapopulations [7, 8]: in contact zones, it may have important consequences because it can break down the genetic integrity of locally adapted populations and counter-balance the strength of selection. Several processes can limit gene flow despite long-distance dispersal, and thereby increase the efficiency of response to selection. High selfing rates for example may preserve genetic integrity of well-adapted populations, but on the other hand they can also generate inbreeding depression and reduce the effect of selection [9–11]. In plants, temporal assortative mating is usually due to flowering time differences, which can also impede gene flow despite sympatry. If parental species are adapted to different habitats in parapatry, divergence in flowering time can be reinforced in the contact zone, thereby preventing gene flow and maladaptive hybridization. Examples of this can be seen in natural populations of Agrostis tenuis and Anthoxanthum odoratum at mine boundaries [12] and in A. odoratum growing under different treatments in the 150-year old Park Grass Experiment [13]. In contrast, if divergent flowering times are selected in allopatry by different environmental factors, they may overlap in sympatry if intermediate ecotones exist [3] or in the case of habitat disturbance (e.g. [14]). If hybrids do not suffer reduced fitness, the only way to maintain temporal isolation is a variation of selection regimes through the reproductive season [15]. Indeed, as flowering times are often highly heritable, assortative mating due to flowering phenology can strongly influence the response to selection, for example it can increase the rate of response to directional selection [16, 17]. Moreover, in the case of strong temporal isolation, high selfing rates can provide reproductive insurance and thus contribute to maintain this isolation.

Local scale studies involving cross-generational approaches with molecular markers are known to be useful for exploring the interactions between selection, assortative mating and dispersal in natural hybrid populations [18]. However, very few recent studies have used methods such as paternity or mating system analyses to estimate the importance of assortment and/or heterogeneity in mating success in structuring hybrid zones of plants (e.g. [19–21]) or animals (e.g. [22]). Detecting temporal assortative mating, i.e. the correlation in flowering time between pollen donors and recipients, can be accomplished at a local scale by paternity analysis [23]. Here we extend a recently developed mating model [24] to estimate the level of temporal assortative mating, along with other important parameters involved in the evolution of a hybrid zone population between two closely related forest tree species.

The two European ash species Fraxinus excelsior L. and F. angustifolia Vahl have separate distributions in France but occur in sympatry in several contact zones where they hybridize [25], although they show completely disjoint flowering phenologies in allopatry [26, 27]. We previously showed that the extent of hybridization seems to be limited by climatic variations in some regions, but intermediate conditions provide ecotones where hybrids are widespread, such as in the Loire valley in central France [25]. Here we focus on one of the populations of this hybrid zone, in which we have shown that genetic and morphological differentiation of the adult trees correlates with differences in flowering times, producing isolation by time patterns [28]. Individuals with extreme phenologies appear genetically and morphologically similar to one parental species, while intermediate flowering individuals cluster into several hybrid classes with flowering dates between those of the two species. Moreover, we showed that these intermediate flowering hybrids produced more flowers and fruits over the two years of study. If these high levels of flowering and fruiting lead to a higher fitness, we may expect that these hybrid genotypes will rapidly invade the zone, especially if assortative mating occurs and/or selfing is frequent. The two parental species are known to be outcrossing species (e.g. t _m values provide outcrossing rates between 0.94 and 1 for F. excelsior [29, 30], and between 0.95 and 1 for F. angustifolia [Fernandez-Manjarrés and Gérard, unpublished]) but to date, there have been no study examining the level of selfing in hybrids.

In this paper, we use a modified version of the mixed-mating model [24, 31, 32] to explore the relative importance of diverse forces influencing the evolution of this hybrid zone population by estimating jointly: (i) the pollen dispersal kernel and the rate of pollen immigration from outside the population, (ii) the strength of spatial and temporal assortative mating, (iii) the selfing rate within the population, (iv) the relative male fitness through mating success of the parental species and different intermediate flowering hybrids and (v) the effect of different phenotypic trait variations on male fecundity (i.e. sexual phenotype, flowering intensity, tree size and fruit production), and (vi) the effect of floral parasite infection intensity on the male fecundity. We also estimate the variation in selfing rates among phenological groups in order to assess the level of selfing in hybrids compared with the parental species. Additionally, we study the impact of the genotyping error rate assumed in the method, as it can have a strong impact on the estimation of relative mating successes [33, 34].

Results

Joint estimation of the dispersal kernel, temporal assortative mating and male fecundities

Dispersal kernel and immigration

The exponential power dispersal kernel provided here a better fit than the normal or the exponential kernel and than panmixia. Note however that the fit of the exponential power was not significantly better than the exponential kernel assuming a genotyping error rate of 0% and 0.1% (Table 1). The estimated dispersal kernel was fat-tailed: the shape parameters ($\widehat{b}$) were lower than 1 (Table 1). Assuming genotyping error rate of 0.1% or 2.5% decreased the estimates of the scale parameter of the Gaussian and exponential kernels leading to shorter dispersal distances. In contrast, high genotyping error led to fatter-tailed exponential power kernel (i.e. with a smaller shape parameter b) and also to a greater mean dispersal distance (Table 1 and Figure 1). The estimated immigration rates decreased when the rate of assumed genotyping errors increased (Table 1), but was still ~60% for the highest assumed error rate.

Figure 1

Log-plot of dispersal kernels estimated under the Gaussian (dotted lines), exponential (dashed lines) and exponential power (plain lines) models, without error rate (black) and with low error rate (dark grey) and high error rate (soft grey). All kernels were estimated under the complete model (model 1).

Table 1 Dispersal and selfing parameter estimates

Temporal assortative mating

The complete model 1, which modelled relative flowering phenology, provided a significantly better fit than the model without any effect of relative phenology when accounting for genotyping error (Table 2). The highest reproductive successes were observed when males flowered slightly earlier than the trees they fertilize ($\widehat{g}$_-1 and $\widehat{g}$_-2> g₀) or when they belonged to the same phenological group (g₀, fixed here at 1). Relative male fecundities decreased as the pollen-recipient trees flowered earlier than the pollen emitting trees. The fecundity on pollen-recipient trees that differed from 4 phenological groups (i.e. g_-4, the fecundity of F. angustifolia males on F. excelsior females) was estimated at 0. Accounting for possible genotyping errors did not substantially change the range of values of $\widehat{g}$, except for $\widehat{g}$_-1, which was estimated at a much lower value when no genotyping errors were assumed (Table 2).

Table 2 Temporal assortative mating parameter estimates

Male mating success of phenological groups

The complete model 2, where the fertility of each phenological group was estimated over all pollen-recipient trees regardless of their flowering time, provided a significantly better fit than the model without any effect of phenology (P < 0.01, P < 0.001 and P < 0.001 respectively with 0%, 0.1% and 2.5% of genotyping error). The estimated male fecundity was significantly higher for phenological group 2 than for all other groups, and the estimated fecundity of the late flowering "F. excelsior" group 5 was zero when assuming no error (Figure 2). All other parameters (dispersal, fecundity, s and m) did not change compared to model 1.

Figure 2

Estimates of relative fecundities of phenological groups, with or without accounting for genotyping error.

Effect of phenotypic trait variations and gall attacks

Flowering intensity and diameter at breast height (DBH) had a significant effect on the relative male fecundities, regardless of the assumed error rate (Table 2): male fecundities were higher for trees with larger DBH and with larger flowering intensities. No other factor had a significant effect, except gall attacks for the assumed error rate of 0.1% (severely attacked trees had lower fecundities), while fruiting intensity had a marginally significant effect (P-value = 0.06, higher male fecundity for trees producing more fruits). The effect of the sexual type was never significant. The range of estimated relative fertilities did not change substantially when accounting for genotyping errors, except $\widehat{f}$₄ for flowering intensity, $\widehat{f}$₃ for gall attacks with a high error and all fertilities for sexual type (particularly $\widehat{f}$₁ corresponding to the relative fertility of males) (Table 3).

Table 3 Fertility parameter estimates

Male effective population density

With model 1, we estimated from equation (8) the reduction of effective male population density ($\widehat{d}$ _em /d _obs ) when only accounting for variation in phenotypic traits as 0.32 when the assumed genotyping error rate was set at 0%, 0.24 for a rate of 0.1%, and 0.27 for a rate of 2.5%.

The mean reduction of effective male population density due to phenotypic traits and temporal assortative mating was estimated at ($\widehat{d}$ _em /d _obs ) = 0.30, 0.23 and 0.27 (with an error rate of 0%, 0.1% and 2.5% respectively). Among pollen-recipients, regardless of genotyping error, the highest reductions were estimated from the pollen clouds of pollen-recipients belonging to the two latest phenological groups (4 and 5), and particularly that of group 5 (almost two-fold decrease, results not shown).

Finally, the mean reduction of effective male population density due to phenotypic traits, temporal assortative mating and non-random dispersal was estimated to ($\widehat{d}$ _em /d _obs ) = 0.09, 0.07 and 0.06 with an error rate of 0%, 0.1% and 2.5% respectively.

Selfing rates

The overall selfing rate estimated with model 1 or 2 was 10% and was slightly affected by genotyping error (Table 1). It varied among phenological groups (Figure 3) and among families, as estimated by MLTR: the mean outcrossing rate (t _m ) was significantly different from 1 within half of the families, with values ranging from 31% to 88% (i.e. selfing rates ranging from 12% to 69%). The selfing rate was close to 20% in group 2 (t _m = 82.7%, standard deviation = 7.7%). In group 3, the selfing rate was estimated at 7% (t _m = 93%, SD = 0.066): the family-level t _m values were significantly different from 1 within 3/8 families, ranging between 62.5% and 94%. The selfing rate was estimated at zero in the two latest flowering groups (4 and 5).

Figure 3

Outcrossing rates t _m estimated from family arrays sampled on mothers from each phenological group (2 to 5). Standard errors were computed from 1000 bootstrap replicates over families.

Discussion

Methodological insights

According to Araki and Blouin [33], assignment error can have a strong effect on the estimation of the relative reproductive successes of different groups, particularly by increasing type II errors (assignment of an untrue parent) when the proportion of unsampled parents is high. Avoiding these errors is generally difficult [35]. Here, we accounted for genotyping errors more realistically than previous methods [36], where a mistyped allele can be of any size compared to the original allele. Instead, we assumed that a genotyping error could be only of ± 1 repeat unit, which is more realistic for microsatellite loci.

As expected, some estimated parameter values varied with different assumed rates of error. However, the general trends remained quite similar: the estimates of dispersal parameters were relatively stable, and the ranking of relative phenology parameters, as well as fertility parameters of various phenotypic trait classes, were barely affected, except for sexual type. The main difference lay in the significance of effects, which was certainly due to the increased information provided by accounting for genotyping errors.

Similarly, the value of the estimated rate of external pollen (m) when genotyping errors were assumed may be more reliable than without errors. The very high value (~80%) obtained when these errors were ignored illustrates a particular caveat of paternity analyses using microsatellite markers: higher allele numbers generally produce higher estimated immigration rates. Indeed, using paternity exclusion methods, the apparent gene flow was estimated here between 45% and 55% with four loci, depending on the level of polymorphism of the selected loci, and increased up to 70% with six loci (Gérard et al., unpublished results). Increasing the level of polymorphism of the chosen loci improves the stringency of paternity analysis, by reducing the rate of multiple assignments and type II errors. At the same time, however, more polymorphic loci increase the risk of overestimating external pollen flow by increasing the error rate.

Homogenizing effect of long-distance pollen dispersal

We found a high level of pollen immigration, even when typing errors were accounted for (external pollen flow ~60% for the highest assumed genotyping error rate). This has also been found in other anemophilous forest tree species such as Pseudotsuga menziesii [37], or Quercus robur and Q. petraea [38]. Moreover, the best-fitting dispersal kernels were rather fat-tailed, which appears increasingly as a common feature in plants, including perennial herbs (e.g. [39, 40]), crops (e.g. [41]) and forest trees (e.g. [42, 43]). Long-distance dispersal may tend to connect distant populations, homogenizing differentiated gene pools. Here it is likely to have acted on the evolution of the hybrid zone of the two Fraxinus species, by connecting remote populations, increasing local genetic diversity [7] and possibly counteracting the effect of local adaptation if assortative mating is not strong enough. This long-distance dispersal has probably contributed to the creation and maintenance of a large-scale hybrid zone, as detected all along the Loire valley [25].

Isolating role of assortative mating and selfing

Spatial assortative mating

Contrary to long-distance dispersal, short-distance dispersal may limit gene exchange at a local scale. The estimate of the shape parameter b of the exponential power kernel was not much smaller than 1 (the value for an exponential kernel) and the estimated mean dispersal distance was quite low (100 <$\widehat{\delta }$ < 150 m, depending of the typing error rate). These findings are consistent with a previous study we conducted in the same population showing that co-flowering trees were patchily distributed in space [28]. Thus, even if long-distance dispersal takes place at a non-negligible rate, a substantial number of reproductive events inside the stand may occur within the patches. Short-distance may produce pronounced spatial genetic structure, as previously detected in this population [28].

Temporal assortative mating

The homogenizing effect of long-distance dispersal may also be counteracted by temporal assortative mating. Even with incomplete assortative mating, a large part of inside-stand reproductive events occurred within the same phenological group, which is strengthened by the high level of selfing. These patterns may contribute to maintain isolation by time patterns [28] and potentially increase the rate of response to selection [16, 17].

One of our main results is that we detected no crosses between the parental species, as expected considering their phenologies. Thus, intermediate flowering hybrids could represent bridges to gene flow in the contact zone, as described for example for Asclepias species [44]. Incidentally, the absence of mating events observed here between the two distinct species raises the question of how the first hybrids were produced. This could be due to exceptional events of hybridization between the two species, which occur at too low a rate to be detected in our study. For example, exceptionally hard winters, in which F. angustifolia individuals may have flowered much later than usual, as they require particular heat and chill conditions to flower [26], could have favoured the first hybridization events.

Selfing rates

High selfing rates may also contribute to counteract the homogenizing effect of long-distance dispersal, and thus help maintain temporal isolation and increase the success of well adapted genotypes if they do not suffer strong inbreeding depression. Indeed, it has been shown that a mutant exhibiting high selfing rates can invade locally stable outcrossing populations despite strong inbreeding depression [45]. The relatively high rates of self-fertilization of hybrids may provide an advantage for them in colonizing the hybrid zone, as the parental species show very limited selfing: the estimated rates vary between 0 and 6% for F. excelsior [29, 30], and between 0 and 5% for F. angustifolia (Fernandez-Manjarrés and Gérard, unpublished). Given the low selfing rates of parental species, the occurrence of relatively high rates within intermediate flowering groups is surprising. The hypothesis of pollen limitation can be excluded for three reasons: (i) the number of individuals and the mean flowering intensities were higher in the groups where high selfing rates were detected [28], (ii) the reduction in effective male population density due to fertility parameters and phenology was smaller in pollen clouds of mothers from group 2 and 3 ($\widehat{d}$ _em /d _obs ) = 0.31 and 0.33 respectively), and (iii) mothers showing the highest selfing rates were not those that had the strongest reduction in their number of effective sires (results not shown). A possible explanation for this increased selfing may be a breakdown of epistasis in intermediate flowering hybrids, caused by linkage disequilibrium between alleles at loci involved in self-rejection mechanisms that co-evolved independently within the two species.

Male mating success of hybrids and parental species

The fertility of hybrids (i.e. the intermediate phenological groups), and in particular individuals from group 2, is higher than the fertility of either parental species (Groups 1 and 5). Differences in some phenotypic traits may participate in this increased fertility.

We found in our previous study [28] that intermediate flowering hybrids (groups 2–4) produced more flowers and had a large DBH (Table 4). Here we show that larger tree diameters and higher flowering intensities increased male fertility (Table 3), as expected for wind-pollinated species [46], and also detected in other tree species (e.g. [24]). This may contribute to increase the fertility of group 2.

Table 4 Summary of phenotypic trait values in the phenological groups

Surprisingly, we did not detect any effect of the sexual type on fertility, probably due to the similar values of siring success of different classes of hermaphrodites (i.e. individuals with a majority of staminate flowers did not sire more seeds than perfect hermaphrodites). Thus male mating success did not depend on the relative proportions of staminate vs. perfect flowers. Nevertheless, a higher male fertility of males relative to hermaphrodites was detected in controlled crosses [47] and natural populations [30] of F. excelsior. This was retrieved here by the very high relative male fecundity of pure males estimated with a typing error rate but it should be nuanced by the very low frequency of males in the population (2%).

Optimal sex allocation may also contribute to increase the male fertility of hybrids. Indeed, trees producing many seeds also had the highest relative male fecundity. Our results are consistent with classical predictions of sex allocation theory, i.e. a constant optimal sex allocation for simultaneous hermaphroditism, where individuals simultaneously produce male and female gametes [48, 49]. Nevertheless, many plant species seem to exhibit a gradual shift in sex allocation, and thereby in functional gender, with increasing size [46, 50, 51]. A positive correlation between male and female reproductive success has been detected in very few cases only [52–54]. Here we present a case where individuals with a high male mating success also show a high female success. Further long-term work is required to confirm that this observation remains true through time and for other hybrid populations.

The higher relative fertility of group 2 may be partly influenced by its lower mean level of gall attacks, compared with those of the other groups [28]. As gall mites (Eriophyes fraxiniflora) infect male ash flowers [55], it is expected that high attack levels would influence the relative male fecundity. Indeed, low levels of gall attacks seemed to have little effect but high rates strongly decreased relative fecundity.

Asymmetry of pollen flow and evolution of the hybrid zone

Flowering phenology also generates asymmetrical pollen flow: trees are quite successful in fertilizing the individuals with a flowering period that overlaps their own but starts later. In contrast, they show quite limited success on the individuals that flower earlier than they do, even if the flowering times overlap. This may result from the protogyny of the Fraxinus species: early flowering trees may participate in reproduction mainly as pollen donors. Mite attacks may also contribute to the observed asymmetry of gene flow by pollen since late flowering trees are more infected than early flowering trees [28]. They may also favour hybridization, by reducing the male fertility of late flowering trees. This asymmetry may provide a demographic advantage to early flowering hybrids because of the greater distances of pollen dispersal compared to seed dispersal in forest trees. On the other hand, although F. angustifolia individuals (Group1) had poor fertility during the year of study (which is likely due to higher susceptibility to late winter frosts, also affecting their fruit production), they may benefit from favourable years and contribute significantly to young seedling generations (Gérard et al. in preparation). Moreover, F. angustifolia shows reduced dormancy compared to F. excelsior [56], possibly conferring another demographic advantage. Thus, the invasive potential of F. angustifolia through hybridization in this region may be highly modified in the case of global warming.

Further work, however, is still needed to assess the relative fitness of seeds produced by different types of individuals, either selfed or outcrossed. Moreover, as pollen dispersal seems to occur over large distances, a larger sampling effort is perhaps needed to get more accurate estimates of the dynamics of the metapopulation in this hybrid zone. As the evolution of partially cross-fertile plant communities is greatly influenced by the strength of assortative mating and demographic characteristics [57], theoretical work is also needed to better understand the interaction of short- and long-distance dispersal and assortative mating, as well as environmental fluctuations (e.g. climate) in plant hybrid zones.

Conclusion

Temporal and spatial assortative mating limit gene flow in this hybrid zone population, even if long-distance dispersal should tend to counter-act their effects. Gene flow between parental species does not occur and intermediate flowering hybrids apparently represent bridges to gene flow between them. Early flowering hybrids, which have the highest male mating success, show optimal sex allocation which, with increasing selfing rates, can potentially confer to them a selective advantage in the hybrid zone. Moreover, temporal assortative mating could contribute to increasing the rate of response to selection by limiting gene flow between different classes of individuals.

The asymmetry of gene flow coming from early flowering pollen donors into late flowering recipients is probably a key factor involved in the dynamics and evolution of this hybrid population. If climate warming allows F. angustifolia not to suffer from winter frosts, the presence of hybrids could contribute to extending its range through this asymmetry. This study has strong implications for understanding the dynamics of forest hybrid zones and for the management of forest diversity in a climate change context.

Methods

Focus species and sampling

Common ash (F. excelsior L.) and narrow-leaved ash (F. angustifolia Vahl) are closely related species [58] with contrasted distributions across Europe: F. angustifolia has a Southern Mediterranean distribution, whereas F. excelsior occurs at more Northerly latitudes. They are in sympatry in several regions in France, such as the Loire and Saône valleys [25, 59]. Both species are post-pioneer forest trees with a colonizing behavior and a discontinuous spatial distribution. They require abundant water, especially F. angustifolia which is often found in low elevation riparian forests, particularly in central France [59, 60]. Both are protogynous, although anther dehiscence occurs while the stigma is still receptive, and pollen and fruits are wind-dispersed. Fraxinus excelsior has a complex trioecious breeding system: sexual types vary across a continuum from pure male individuals to pure females, with all kinds of hermaphrodites in between [47, 55, 61]. Much less is known regarding the mating system of F. angustifolia, which rather seems to be androdioecious [28]. The species differ in their phenology: initiation of flowering occurs between mid-March and early April for F. excelsior, and between mid-December and late January for F. angustifolia, with larger among-year variation [26, 27].

The study site is an unmanaged, naturally regenerated population covering almost 7 ha, located at Saint-Dyé-sur-Loire (latitude 47° 39' 10" N, longitude 1° 28' 40" E, elevation 78 m.a.s.). It is along the Loire river in central France in a zone of sympatry between the two species (see [28] for all details on the stand along with a map). The population was composed of 269 flowering adult trees, and we observed a unimodal distribution of flowering dates, which ranged between the extreme phenologies of the two species. All trees were assigned to one of five phenological groups according to their flowering date (i.e. bud flush, just before pollen offset, 1 group = 2 weeks). These five phenological classes were validated by a Canonical Discriminant Analysis on morphological traits and microsatellite allelic frequencies, and by linkage disequilibrium estimation within and among groups ([28] and Table 4). We also measured the diameter at breast height (DBH), and grouped them into five discrete classes (< 40 cm to > 160 cm, with 40 cm intervals). We finally measured for each individual two fitness components, the flowering and fruiting intensity (each as a visually-assessed 5-level class variable), as well as the intensity of floral gall attacks (also as a visually-assessed 5-level class variable: class 0 includes trees without any gall whereas in class 4, more than 80% of flowers were infected) (see Table 4). For each individual, we also recorded its sexual type, according to one of the four following categories: MM, MH, HM, HH, which correspond either to pure males (MM) or to three types of hermaphrodites, with either a majority of male flowers (MH) or a majority of hermaphroditic flowers (HM) or only hermaphroditic flowers (HH). We followed the classification of Morand-Prieur [30], who also observed three additional categories of trees (three types of females) that we did not observe here.

We harvested 432 seeds on 27 fruiting trees (16 seeds per tree) in the autumn of 2003. The distances between the sampled mother trees ranged from 14 to 1775 m, with a mean of 538 m. We tried to sample individuals of each phenological group, but the earliest flowering trees (first group) did not produce sufficient numbers of fruits. Seeds were rehydrated and sterilized as described in Raquin et al. [62], and embryonic tissues were dehydrated in a 1:1 ethanol-acetone solution. Total DNA was extracted from dried tissue using a DNeasy^® 96 Plant Kit (Qiagen).

All flowering individuals in the population and all seeds were genotyped at eight microsatellite loci: M-230, FEMSATL 4, 8, 10, 11, 12, 16 and 19 [63, 64], under conditions described by Morand et al. [65] and Gérard et al. [28]. The theoretic exclusion probability over the eight loci was 0.9995 for the 269 flowering adult trees.

Mating model analysis

Paternity analysis allows one to gauge the heterogeneity in male fecundity by estimating the effect of ecological and/or phenotypic variables on this fecundity: for example the effect of the floral phenotype [66], of the inflorescence morphology [67] or of the floral phenology [37]. Here we adapted a mating model developed by Oddou-Muratorio et al. [24], which stems from the neighbourhood model [31, 32]. This method allows one to estimate jointly the dispersal curve, level of pollen immigration, selfing rate and the heterogeneity in male fertility. It avoids type I and type II errors occurring in categorical paternity assignment (for a review see [68]).

Modeling pollen clouds

The model considers that each offspring o sampled from a given mother j _o can result either from: (i) self-fertilization (with probability s), (ii) cross-pollination by a male located outside the study population (with probability m), or (iii) cross-pollination by a male sampled within the study area (with probability 1- m- s).

The probability that an offspring o of mother j _o results from self-fertilization and has genotype g _o depends on inheritance probabilities only:

P _self (o, j _o ) = sT(g _o |$g_{j_o}$, $g_{j_o}$)     (1)

where T(g _o |$g_{j_o}$, $g_{j_o}$) is the Mendelian segregation probability [69] of the offspring genotype g _o given the mother genotype $g_{j_o}$.

The probability that an offspring o of mother j _o results from cross-pollination by a male located outside the study area, and is of genotype g _o , depends on the allelic frequencies in the "outside" populations as follows:

P _outside (o, j _o , AF) = mT(g _o |$g_{j_o}$, AF)     (2)

where AF corresponds to the allelic frequencies in the immigrant pollen cloud entering the study population. These external allele frequencies were inferred here from the retrieved paternal gametes of all offspring without compatible male parent within the study population, using the software MLTR 3.0 [70].

The probability that an offspring o of mother j _o results from cross-pollination by a male l with genotype g _l sampled within the set Ψ of males occurring in the study population, and has the genotype g _o depends on the contribution of male l to the pollen cloud as follows:

$P_{inside}(o,j_o,\Psi )=(1-m-s){\displaystyle \sum _{l\in \Psi }{\pi }_{j_o,l}({\theta }_d,\text{F},\text{G})T(g_o|g_{j_o},g_l)}\left(3\right)$

where π _jk (θ _d , F, G) is the expected proportion of pollen emitted by male k in the local pollen cloud that fertilizes female j, θ _d is a vector of dispersal parameters, F is a vector of fecundity parameters for all phenotypic classes of males (including parasite attack intensity). In addition to Oddou-Muratorio et al. [24], we included here a vector G containing the probabilities of mating between flowering classes (see below). The proportion π _jk depends (i) on the spatial location of trees (dispersal), via the probability p _jk that a pollen grain emitted by male k travels to mother j, computed from the dispersal kernel given below, (ii) on the phenotypic variables, via the fecundity function Φ _k of male k and (iii) on flowering phenology (assortative mating) via the probability γ _kj that a pollen grain is emitted by male k when the ovules of mother j are receptive, following the equation:

${\pi }_{jk}=\frac{{\Phi }_k{\gamma }_{kj}{p}_{jk}}{{\displaystyle \sum _{l\in {\Psi }_M}{\Phi }_l{\gamma }_{lj}{p}_{jl}}}\left(4\right)$

The p _jk 's, Φ _k 's and γ _kj 's are computed respectively from the dispersal kernel and the F and G vectors (see below).

Model for the dispersal kernel

We modeled pollen dispersal using a dispersal kernel p(.;x, y) that describes the probability for a pollen grain emitted in (0,0) to participate in the pollen cloud fertilizing a tree in (x, y) [7, 42]. We used the family of exponential power functions:

$p({\theta }_d;x,y)=\frac{b}{2\pi a^2\Gamma (2/b)}\exp \left[-{\left(\frac{\sqrt{x^2+y^2}}{a}\right)}^b\right]\left(5\right)$

where Γ is the gamma function [71]. θ _d includes the parameters a and b for an exponential power dispersal kernel, but only a for a Gaussian or exponential dispersal kernel (for which b is set at 2 or 1 respectively, see [42]). The shape parameter b affects the tail of the dispersal function and the scale parameter a is homogeneous to a distance [72]. The mean distance δ traveled by a pollen grain under the kernel p(a, b) is given by:

$\delta =a\frac{\Gamma (3/b)}{\Gamma (2/b)}\left(6\right)$

For b < 1, the dispersal kernel is fat-tailed [72], i.e. it decreases more slowly at long distance than an exponential kernel, implying the possibility of long-distance dispersal. For b > 1 (e.g. b = 2 for the Gaussian model), the dispersal is thin-tailed, which implies much less long-distance dispersal events. For 1 <b < 2, the dispersal kernel is leptokurtic (i.e. a distribution with a high peak and fewer long-distance events than the Gaussian one) and it is platykurtic for b > 2 (i.e. a flat-topped distribution) (see [42]).

For each dispersal kernel, i.e. for each set of dispersal parameters θ _d = (a, b), the probability p _jk used in equation (4) that a pollen grain emitted by male k travels to mother j is given by p _jk = p(θ _d ; x _k -x _j , y _k -y _j ). We thus assume that pollen emitted by each sampled male tree disperses according to the same dispersal kernel p(θ _d ; x, y).

Model for male fecundities

Male fecundity was modeled following Oddou-Muratorio et al. [24], considering discrete classes for both phenotypic traits and gall attacks. Denote $f_c^i$ is the relative fecundity of a male belonging to the c^th class of a given trait i. The fecundity Φ _k of a male k is thus given by $f_{c(k)}^i$ where c(k) is the phenotypic class to which male k belongs. We did not consider any interaction effect between phenotypic traits and gall attacks: for example, the fecundity of a male belonging to class c(k) of flowering intensity (trait 1) and suffering a level e(k) of gall attacks (trait 5) is thus simply: ${\Phi }_k=f_{c(k)}^1\ast f_{e(k)}^5$. All relative fecundities $f_c^i$'s of the different phenotypic classes of males are stored in the vector F.

Models for flowering phenology

We modeled the flowering phenology in two different ways. In model 1, we modeled the relative flowering phenology considering the phenological differences between pollen-emitting and pollen-recipient trees to assess the strength of assortative mating. In that case, we define g _d as the relative male success of a tree at siring offspring on another tree when the phenological difference between the two trees is d. The relative male fecundity γ _kj of tree k on tree j is thus g_d(kj), where d(kj) is the difference in phenology between pollen-emitting tree k and pollen-recipient tree j, d(kj) being zero if k and j flower simultaneously, positive if the pollen-emitting tree k flowers before the pollen-recipient tree j and negative if k flowers after j. The vector G contains all g _d 's.

In model 2, we applied to the phenology the approach described above for other phenotypic traits to assess the overall relative male mating success of the different phenological groups on all mothers of the population. In that case, the vector G, containing the relative fecundities of the different phenological groups, is defined in the same manner as F.

Joint estimation of the dispersal kernel, male fecundities and phenological assortative mating

We used a maximum likelihood approach to estimate jointly the vector of dispersal parameters (θ _d ), the vector of the relative male fecundities of the different phenotypic classes and infection intensities (F), the vector of the relative fecundities of different classes of flowering time or flowering time differences (G), the selfing rate (s) and pollen immigration rate (m). The log-likelihood function of all observed progenies collected from females, given the above model is given by:

$\log L({\theta }_d,\text{F},\text{G},s,m)={\displaystyle \sum _{o=1}^o\log }\left[P_{self}(o,j_o)+P_{inside}(o,j_o,\Psi )+P_{outside}(o,j_o,AF)\right]\left(7\right)$

Equation (7) assumes that all fertilization events are independent of each other.

Statistical analyses

All analyses described below were performed using Mathematica 5.0 (Wolfram Research Inc.). Notebooks are available upon request from PRG.

Fits

All fits were achieved by maximizing the log-likelihood of equation (7) following a quasi-Newton algorithm. For the exponential power dispersal kernel parameters, we estimated the mean dispersal distance δ and the shape parameter b rather than a and b. To ensure numerical convergence, parameters m and s were transformed through a logit function [m = e^m'/(1+e^m') and s = e^s'/(1+e^s')] and male fecundity parameters were transformed through $f_{c(k)}={10}^{{f^{\prime }}_{c(k)}}$.

Confidence intervals

We obtained the 95% likelihood-profile confidence intervals for δ and b, by plotting contour plots of the likelihood function [73]. For the vectors of parameters F and G (male fecundities and flowering phenology), and for the parameters m and s, we derived 95% confidence intervals by computing the asymptotic Gaussian variance-covariance matrix, which is the inverse of the Fisher's information matrix (i.e. the opposite of the expectation of second-order partial derivatives of the log-likelihood function with respect to all couples of parameters) [73]. As these parameters were estimated through a transformation function, we first computed a symmetric Gaussian asymptotic interval for the transformed parameters (m', s' and f') by using a delta method [73]. We then obtained asymmetric confidence intervals for the parameters of interest (m, s, and f) by the reverse transformation of the bounds of the intervals.

Tests

We tested for the effect of dispersal kernel, spatial and temporal non-random mating, and of the variation in phenotypic traits and gall attacks on male mating success, by building several corresponding nested models from equation (4) and using likelihood-ratio tests (LRT) in a type III approach [24].

Effect of genotyping error

To account for possible genotyping error, we re-computed the Mendelian transition probabilities allowing an error of ± 1 microsatellite repeat unit for genotypes of all males at all loci (see Appendix). All models described above were then used with the modified log-likelihood accounting for a fixed mistyping rate, either a low rate of 0.01%, or a relatively high rate of 2.5%.

Effective male population density

The reduction in effective male population density can have strong consequences: it may for example increase drift in natural populations and influence the rate of accumulation of beneficial and deleterious mutations [24, 74]. In fact, unequal male fecundities and asynchronous flowering may drastically reduce the effective male population density [74]. Using male fecundities estimated as described above, we thus estimated the reduction in effective male population density caused by heterogeneity in male fecundity under the complete model following Oddou-Muratorio et al. [24]. This density is defined by d _em = N _em /A, where N _em is the effective male population size (i.e. the inverse of the probability that two pollen grains come from the same male) and A is the area covered by the study population. The reduction in effective male population density can be expressed as:

$\frac{d_{em}}{d_{obs}}=\frac{1}{N}\frac{{\left({\displaystyle \sum _{k\in \Psi }\left({\displaystyle \prod _{i=1}^n{f}_{c(k)}^i}\right)}\right)}^2}{{\displaystyle \sum _{k\in \Psi }{\left({\displaystyle \prod _{i=1}^n{f}_{c(k)}^i}\right)}^2}}\left(8\right)$

where d _obs is the observed male population density, $f_{c(k)}^i$ is the fecundity of male k belonging to the c^th class of the phenotypic trait i, and n is the number of traits under consideration. We also estimated the reduction in effective male population density due to dispersal features and phenology, including in the above sums all p _jk 's and all g_d(kj)'s.

Estimation of selfing rates within phenological groups

We performed a mixed-mating model analysis using the software MLTR 3.0 [70], to assess the relation between selfing rate and flowering time, by estimating the mean outcrossing rates t _m in each phenological group (from the second to the fifth, because we had no seeds for the first group) from the multilocus genotypes of the seeds harvested on the mother trees, using an EM algorithm. Standard errors were computed by performing 1000 bootstrap replicates using families as resampling units. We also estimated the family-level t _m values (i.e. the mean outcrossing rates among seeds harvested from each mother-tree) [70].

Appendix

Modelling the genotyping error rate

If α was the error rate and if the paternal allele of the offspring was unambiguously known, we assigned as father with a probability α, a male homozygous for an allele that differ by 1 repeat unit from the offspring, and with a probability 1–2α, a male homozygous for the same allele.

This means that in equation 8, the standard transition probabilities T(g _o |$g_{j_o}$, g _l ) defined by Meagher (1986) were replaced by $T(g_o|g_{j_o},g_l)={\displaystyle \prod _{loci}{T}_{loc}(g_{o,loc}|g_{j_o,loc},g_{l,loc})}$, where for each locus loc, and with the notations from [70]:

$T_{loc}\left(\{o_1,o_2\}\left|\right\{m_1,m_2\},\{f_1,f_2\}\right)=\frac{D_{o_1}^{m_1{m}_2}{\tilde{D}}_{o_2}^{f_1{f}_2}+D_{o_2}^{m_1{m}_2}{\tilde{D}}_{o_1}^{f_1{f}_2}}{\left(1+{\delta }_{o_1{o}_2}\right)},$

with δ _ij = 1 if i = j and 0 if i ≠ j, $D_o^{m_1{m}_2}=\left(\frac{{\delta }_{o{m}_1}+{\delta }_{o{m}_2}}{2}\right)$ is the Mendelian contribution from the mother, without any mistyping, and the contribution from the father accounting for possible mistyping:

${\tilde{D}}_o^{f_1{f}_2}=\left(\frac{(\alpha {\delta }_{o{f}_1+1}+\alpha {\delta }_{o{f}_1-1}+(1-2\alpha ){\delta }_{o{f}_1})+(\alpha {\delta }_{o{f}_2+1}+\alpha {\delta }_{o{f}_2-1}+(1-2\alpha ){\delta }_{o{f}_2})}{2}\right)$