Sex-specific spatial variation in fitness in the highly dimorphic Leucadendron rubrum

Sexual dimorphism in plants may emerge as a result of sex-specific selection on traits enhancing access to nutritive resources and/or to sexual partners. We here investigated sex-specific differences in selection of sexually dimorphic traits and in the spatial distribution of effective fecundity (our fitness proxy) in a highly dimorphic dioecious wind-pollinated shrub, Leucadendron rubrum . In particular, we tested for the effect of density on male and female effective fecundity. We used spatial and genotypic data of parent and offspring cohorts to jointly estimate individual male and female effective fecundity on the one hand and pollen and seed dispersal kernels on the other hand. This methodology was here adapted to the case of dioecious species. Explicitly modeling dispersal avoids the confounding effects of heterogeneous spatial distribution of mates and sampled seedlings on the estimation of effective fecundity. We also estimated selection gradients on plant traits while modeling sex-specific spatial autocorrelation in fecundity. Males exhibited spatial autocorrelation in effective fecundity at a smaller scale than females. A higher local density of plants was associated with lower effective fecundity in males but was not related to female effective fecundity. These results suggest sex-specific sensitivities to environmental heterogeneity in L. rubrum . Despite these sexual differences, we found directional selection for wider canopies and smaller leaves in both sexes, and no sexually antagonistic selection on strongly dimorphic traits in L. rubrum . Many empirical studies in animals similarly failed to detect sexually antagonistic selection in species expressing strong sexual dimorphism, and we discuss reasons explaining this common pattern. 20 21 22 23 24 25 26 27 28 29 30


Introduction
Plant species with separate sexes are relatively uncommon (i.e. 5-6%, Renner, 2014). Separate sexes have nonetheless evolved repeatedly among flowering plants (Renner, 2014), and such transitions often given rise to the evolution of morphological differences between sexes (Geber et al., 1999;Puixeu et al., 2019). The degree of sexual dimorphism has also switched multiple times from low to high along the evolutionary history of certain dioecious plant lineages (Tonnabel et al., 2014). Both sex-specific costs of reproduction and male-male competition to access ovules have been suggested as potential forces causing the evolution of such dimorphism. These two factors could trigger sexually antagonistic selection (Delph & Ashman, 2006;Moore & Pannell, 2011), whereby selection exerts forces in opposite directions in each sex towards sex-specific optima (Cox & Calsbeek, 2009). The goal of this study is to estimate sex-specific fitness as well as the strength and form of selection acting on morphological traits in each sex, in a highly dimorphic dioecious windpollinated plant species. To do so, we combine, and adapt to the case of dioecious species, recently developed statistical methods estimating effective fecundity, a proxy for fitness, and its dependence on morphological traits, while explicitly modeling various spatial effects that could bias such estimations.
The sex-specific cost of reproduction hypothesis posits that sexes should diverge in morphology to satisfy their respective reproductive needs (Freeman et al., 1976;Delph & Bell, 2008). Such divergence can emerge when reproduction involves a stronger cost in one sex than in the other, or when the reproductive costs of each sex imply different resource 'currencies' (Freeman needed for their respective reproduction. The cost of reproduction is generally higher in females than in males, at least considering the cost per reproductive structure. However, at the scale of the whole plant, this trend is often reversed in wind-pollinated plants, which produce large amounts of pollen (Obeso, 2002;Harris & Pannell, 2008;Tonnabel et al., 2017). In some dioecious species inhabiting fire-prone environments, the cost of reproduction markedly differs between sexes because females need to maintain a canopy-stored ('serotinous') seed bank (released by fire). As water intake is necessary to maintain cones closed and prevent seed release during unfavorable period between two fires (Martín-Sanz et al., 2017), we may expect selection for enhanced water conduction to have favored a divergent plant architecture between sexes. Consistent with this prediction, the evolution of longer maintenance of cones in the canopy is indeed associated with the evolution of higher sexual dimorphism in the genus Leucadrendron (Harris & Pannell, 2010).
The evolution of sex-specific differences in vegetative traits may also originate from selection for male morphologies that are better at dispersing pollen and therefore at accessing mates (Tonnabel et al., 2019a,b). The male-male competition hypothesis postulates that male reproduction, by being mostly limited by mating opportunities, selects for males that exhibit traits enhancing their competitive abilities (Bateman, 1948;Arnold, 1994). Several studies have pinpointed the importance of male-male competition in shaping male reproductive and floral traits. These studies showed more extravagant floral displays in males than in females to attract pollinators (e.g. Bond & Maze, 1999;Elle & Meagher, 2000;Wright & Meagher, 2004;Delph & Ashman, 2006;Waelti et al., 2009;Schiestl & Johnson, 2013;Dorken & Perry, 2017), variation in male flowering phenology to track the female phenology (Delph & Herlihy, 2011;Forrest, 2014), increased pollen grain competitive performance in response to polyandry (Lankinen et al., 2017) and morphological evolution of structures that disperse pollen, which prevents the attachment of competing pollen to the pollinator (Coccuci et al., 2014). In wind-pollinated plants, sexual selection may also target vegetative traits such as plant size, branch length or the length of male flower peduncles, which can facilitate pollen dispersal (Klinkhamer et al., 1997;Eppley & Pannell, 2007;Pickup & Barrett, 2012;Harder & Prusinkiewicz, 2013;Tonnabel et al., 2019b). Wind-pollinated plants tend to evolve larger degrees of sexual dimorphism than insect-pollinated lineages, because pollinators require similarity of floral morphology to successfully transfer pollen (Tonnabel et al., 2014;Welsford et al., 2016). Wind-pollinated plants, which typically show large inter-individual variation in fertility (Schoen & Stewart, 1987;Ahee et al., 2015), may thus be particularly subject to sexually antagonistic selection.
In the presence of genetic variation for sexual dimorphism, each sex should, in principle, ultimately reach its optimal trait value, and thereby resolve sexually antagonistic selection (Lande, 1980). Yet, a shared genetic basis of traits between sexes may temporally constrain the evolution of their morphological divergence (Lande, 1980). In constant and homogeneous environments, theory predicts that, with strong positive genetic correlations between sexes, opposite directional selection gradients between sexes should emerge early during adaptation and persist for a long time before the sexual conflict is resolved (Lande, 1980;Connallon & Hall, 2016). Consequently, one would expect evidence for antagonistic selection between sexes to be relatively common. The compilation of numerous sex-specific selection gradients in animals showed, however, a large diversity of patterns, including cases of aligned selection across sexes (Cox & Calsbeek, 2009). In plants, documented patterns of sex-specific selection provided mixed support for sexually antagonistic selection: sexspecific selection gradients have been found to have opposite signs in both insect-and windpollinated species (Delph et al., 2011;Castilla et al., 2014;Tonnabel et al., 2019b) but to have the same sign in other studies (Oddou-Muratorio et al., 2018;Barrio & Teixido, 2014). More recent theory suggests that temporal and spatial variation in selection pressures may explain the lack of signal for sexually antagonistic selection, despite differences in the optimal phenotypes of males and females (Connallon, 2015;Connallon & Hall, 2016;Zafitscher & Connallon, 2017).
Estimating sex-specific selection gradients requires, first, estimating male and female individual fitness, and, second, relating trait values and fitness estimates. Using genotypes of established seedlings and their potential parents, traditional methods first achieve categorical 5 parentage assignments to then estimate individual realized reproductive successes used as fitness estimates. In the next generation methods, genotypes are combined with spatial localization of sampled individuals, through spatially explicit mating models (SEMMs), to disentangle the effect of fecundity from that of the distance between mating pairs (and the distance between mothers and seedlings) on reproductive success (e.g. Oddou-Muratorio et al., 2005). To do so, dispersal is explicitly modeled and dispersal kernels are estimated for both seeds and pollen. A Bayesian method was introduced in this framework to estimate individual male and female effective fecundities (MEMM, Klein et al., 2008 for seed sampling on mother trees; MEMMseedlings, Oddou-Muratorio et al., 2018 for seedling sampling designs). This method considers the likelihood of genotypes conditional on the position of seedlings, so it is unaffected by any process acting on the distribution of seedlings, be it the potential parents' positions, or habitat suitability and disturbances. It analyses seedling genotypes in terms of (1) overall reproductive contribution of each potential parent as determined jointly by gamete production, gamete fertilization rates, seed maturation and germination, and seedling survival until census; (2) dispersal events in terms of estimated dispersal kernels and (3) pollen or seedling competition by a mass action law. Effective fecundity refers only to relative values of the first component for each parent. It varies with, e.g., male-male competitive ability through differences in overall pollen production and their subsequent ovule fertilization abilities, but not with competitive effects dependent on the composition of competitors within the pollen cloud generated by uneven spatial distribution of mates. We here extend this methodology to dioecious species. This spatially explicit approach avoids spatial bias in effective fecundity estimation, typically generated by sampling seedlings non-uniformly with respect to the positions of their parents or by the confounding effects of heterogeneous spatial distribution of mates (Oddou-Muratorio et al., 2018). Used as a proxy for fitness, effective fecundity thus provides the expected relative reproductive success if putative mates (for male fecundity) and regeneration sites (for female fecundity) were uniformly distributed in space, and all offspring could establish and be sampled (Klein et al., 2013). It therefore attenuates the impact of stochastic effects associated with sampling methods on fitness estimates.
Relating fitness estimates to plant traits using the selection gradients methodology proposed by Lande & Arnold (1983) can further suffer from specific statistical bias in sessile organisms living in heterogeneous environments. Indeed, small-scale spatial variation in resources fundamental to plant physiology, including sex-specific reproduction, is common across a range of habitat types (Silvertown et al., 1999;Araya et al., 2011). To disentangle the fitness effect of plant characteristics (such as their ability to harvest resources, which may be sex-specific) from that of the environment (such as the spatial distribution of resources), the spatial distribution of individuals must be accounted for (Rausher, 1992). Indeed, not modeling explicitly the spatial autocorrelation of unmeasured ecological variables affecting fitness can lead to detect false-positive effects of traits on fitness, as on any other response (Guillot & Rousset, 2013). To address this problem, we fitted mixed-effect models with spatially autocorrelated random-effects, using the spaMM package (Rousset & Ferdy, 2014). To our knowledge, it is the first time that spatial effects are taken into account in the estimation of selection gradients. Moreover, spatial variation in plant density and the local sex ratio may generate spatial variation in competition for resources, which can be studied by analyzing their fixed effects on plant fitness. In conclusion, our MEMMseedlings model controls for spatial confounding effects on estimations of effective fecunditity relative to spatial sampling biases and to the spatial heterogeneity of plant distribution (potentially impacting competition among males), while our spaMM procedure controls for spatial environmental heterogeneity.
We applied our methodology to Leucadendron rubrum, a dioecious wind-pollinated serotinous shrub, endemic to the fire-prone South-African fynbos. L. rubrum displays extreme sexual dimorphism (Fig. S1), with males being typically more highly branched, having smaller leaves and taller stature than females (Harris & Pannell, 2010;Welsford et al., 2014;Welsford et al., 2016). A single recruitment pulse typically occurs after fire, killing all adult trees and releasing seeds stored in their canopy seed bank (Cowling & Lamont, 1987). This particular life-cycle allows the estimation of life-time effective fecundity by sampling seedlings only once (i.e. after the fire event). Furthermore, because recruitment is synchronized by fire, all sampled adult individuals in the population have the same age (Bond, 1984). The strong sexual dimorphism of L. rubrum has been previously hypothesized to be the consequence of sex-specific resource requirements (Harris & Pannell, 2010). Indeed, the cost of reproduction in L. rubrum is likely to differ strongly between sexes due to the cost of maintaining the canopy-stored ('serotinous') seed bank in females (Martín-Sanz et al., 2017). We hypothesized that, because of such maternal care, female fitness may be more sensitive to water availability than male fitness. Owing to these differences in resource requirement for reproduction, we therefore tested whether male and female effective fecundities (as defined above) display different spatial structure and whether the observed strong sexual dimorphism is associated with sex-specific selection gradients (Lande, 1980).

Study species and site
Leucadendron rubrum is a dioecious wind-pollinated shrub species endemic to the Western Cape of South Africa (Rebelo, 2001) where natural fires occur every 10-15 years (van Wilgen et al., 2010;Kraaij et al., 2011). L. rubrum belongs to the Proteaceae family and flowers from August to September. Seed recruitment is constrained to a short period following fires, and seeds released between fires typically fail to establish due to competition (Bond, 1984). L. rubrum typically starts flowering at 2-3 years, and seeds are retained in woody cones for several years (Harris & Pannell, 2010). Seeds therefore form a 'serotinous' seed bank, which persists until fire kills the plants, allowing cone opening and wind dispersal of fruits via their plumed perianth (Williams, 1972;Rebelo, 2001). In serotinous species, disruption of water intake to the cone (caused by broken branch or plant death) was shown to lead to seed release suggesting a water cost to keep the cones closed (Treurnitcht et al., 2016;Martín-Sanz et al., 2017). Thus, mother death or any event leading to cone opening before the advent of fire results in seed release in poor conditions for recruitment and ultimately in the loss of progeny.
The study population was located at Bainskloof pass (33°32'21.25''S 19°10'12.10''E) and was contained in a rectangle of 135x110 meters (Fig. S2). We studied all adult individuals of our focal population. Another population of L. rubrum was located at a distance of 310 meters (smallest distance found between two plants from the two populations). All adults of L. rubrum of the focal population (i.e. 86 females and 88 males) were mapped (see Supplementary Methods S1), sampled for DNA analyses and measured for several traits in 2004. In summer 2006, a fire burnt the population, killing all adults. A total of 1,265 seedlings were mapped, and their leaves were sampled in the following Fall (February 2007), four to five months after the fire. The spatial distribution of adults and seedlings is shown in Fig. S2. In one part of the study site, a ditch had been dug for construction after seeds had dispersed, so we were unable to determine the undisturbed positions of 172 seedlings located in that area which we therefore eliminated from the dataset. However, the effective fecundity of adults in this area can still be estimated without bias induced by the disturbance, thanks to the use of MEMMseedlings (see below).
When sampling seedlings, the presence of seedlings from another closely related sympatric species (Leucadendron salignum P.J. Bergins) rendered the identification of L. rubrum juveniles difficult. To ensure that only juveniles of L. rubrum were included in later analyses, we genotyped juveniles (see below for details on the genotyping protocol), as well as adults from both species.
This analysis aimed at assigning seedlings to either species and its results are described in Supplementary Method S2 and Fig. S3. We did not find evidence for the existence of hybrids between L. rubrum and L. salignum. In addition, 254 juveniles were excluded after genotyping as they belonged to L. salignum.

Measurements of adult traits in the field
For adult shrubs, we measured in 2004 three traits describing plant architecture and three traits describing leaf morphology (available at doi:10.5061/dryad.ngf1vhhst). All six traits are known to be sexually dimorphic in this species (Harris & Pannell, 2010;Welsford et al., 2014Welsford et al., , 2016 Fig. S4).
The three traits describing plant architecture were (1) plant height, (2) the first diameter defined as the greatest horizontal diameter of the plant (hereafter, canopy diameter), and (3) the second horizontal diameter defined as the diameter perpendicular to the first diameter (hereafter, second diameter). Several leaves were collected randomly along branches of each adult, dried and photographed. Pictures were analyzed to measure the three traits describing leaf morphology: (1) leaf area, (2) length and (3) width using ImageJ (Schneider et al., 2012). The number of leaves analyzed per adult ranged from 10 to 23 with an average of 20.3. In females, we counted the number of cones in the last two cohorts (cones produced in the last two seasons of cone production and maintained closed since). Because older cohorts were not counted, this measure reflects cone production rather than the maintenance of the serotinous seed bank. We did not count male cones because they were too numerous.

Microsatellite genotyping
We genotyped both adults and their progeny in our focal population (available at doi:10.5061/dryad.ngf1vhhst). For both adults and seedlings, sampled leaves were preserved in silica gel prior to DNA extraction using a modified version of the CTAB protocol (Justy et al., on an ABI3500XL sequencer. The genotypes of all adults and offspring were scored using GeneMapper at the eight microsatellite loci, which exhibited between five and 23 alleles (Table S1 for information per marker). After excluding individuals that did not amplify, our dataset contained 82 females, 85 males and 869 juveniles, corresponding to an amplification failure of roughly 8% for both adults and juveniles. For each microsatellite marker, we used CERVUS (version 3.0.7; Kalinowski et al., 2007) to estimate the non-exclusion probabilities of the first parent, the second parent and parent pairs, to test for Hardy-Weinberg equilibrium, and to compute null allele frequencies (Table S1). Non-exclusion probabilities correspond to the probabilities that the set of loci will not exclude an unrelated candidate parent (or parent pair) of an arbitrary offspring. Finally, we used NMπ (Version 1.0, Chybicki, 2017) to estimate the per-marker genotyping errors.

Joint estimation of effective fecundities and both pollen and seed dispersal kernels
We used a method that uses information about the genotype and the spatial location of adults and seedlings to jointly estimate pollen and seed dispersal kernels and the individual male and female effective fecundities -a proxy for fitness (see the Introduction). Our method builds on the Mixed Effects Mating Model, MEMMseedlings (Oddou-Muratorio et al., 2018), which models mating and dispersal events in a hermaphroditic plant population to estimate the selfing rate, immigration rates and dispersal kernels for both pollen and seeds as well as the variance in male and female effective fecundity (i.e. using random individual effects). We modified the MEMMseedlings algorithm to produce a new version that handles separate sexes, with distinct spatial distributions of male and  (1) in Supplementary methods S3), which is therefore unitless.
Given the life cycle of L. rubrum and the sampling of seedlings after the fire, our estimated effective fecundity integrates the effect of variation between individuals, either male or female, not only for pollen and ovule production, but also for pollen export, fertilization rate, seed maturation and dispersal, maintenance of seeds within the cones (degree of serotiny), adult survival until the fire, seed germination and juvenile survival until the seedling census. Seeds released after the fire and previously stored in the canopy were potentially fertilized in different years. Because adults of L. rubrum do not survive fire, the establishment of progeny after a fire thus represents their entire lifetime reproduction.
Briefly, our MEMMseedlings model combines genotypes data and spatial distribution data for both adults and offspring to estimate, in a Bayesian framework, individual male and female effective fecundities (F k and R j , respectively), the seed immigration rate (m s ), the pollen immigration rate (m p ), the rate of pollen export to the pollen cloud of non-local mothers (υ), the mean seed dispersal distance (δ s ), the mean pollen dispersal distance (δ p ), a parameter affecting the shape of the seed dispersal kernel (b s ) and that of the pollen dispersal kernel (b p ). The estimation of pollen and seed immigration rates (m s and m p ) and, the rate of pollen export to the pollen cloud of non-local mothers, (υ) depends on the actual process of immigration, but is also affected by the fact that around 8% of the parents were excluded from the analysis due to a lack of amplification.
Immigration rates therefore include the probabilities of maternity and paternity attributable to unsampled parents. Both pollen and seed immigration rates are therefore likely overestimated.
Finally, when computing Mendelian transition probabilities between seedlings and putative parents, MEMMseedlings considers genotyping errors by allowing a parent-offspring genetic discrepancy at a maximum of one locus, and at each locus, the probability to mistype any allele was fixed using the the quality of parentage assignments associated with the dataset, MEMMseedlings computes the posterior probabilities for all seedlings for which (i) both parents are known among the genotyped parents, (ii) only the mother, or (iii) the father is known, and (iv) none of the parents is known.
Estimates of dispersal kernels are based on dispersal events within the population, even if the parametric model implies dispersal in unbounded space (beyond the maximal male-female and female-seed distances found in our study population, respectively 106.4 and 139 meters for pollen and seeds dispersal). To describe dispersal within the study population, we computed, from the estimated kernels, the predicted proportion of seeds and pollen that dispersed within the population range and the predicted proportion that dispersed within an arbitrary short distance of 20m.
We estimated the model parameters using two MCMC chains of 50,000 steps and a burn-in of 10,000 steps each. We used uniform prior distributions for the parameters m s , m p , υ, δ s , b s , δ p and sampled every 20 iterations to decrease autocorrelation and averaged after concatenating the two chains. For each of these stored iterations, we also computed the variance in effective fecundity estimates among individuals. Credibility intervals at 95% were calculated for all estimated parameters, as well as the mean value across all iterations and chains (Table 1).

Sex differences in morphology, spatial distribution and analysis of number of cones
We tested for sex differences in morphology using Linear Mixed-effects Models (LMM) with spatially autocorrelated random effects. We analyzed all measured adult traits describing either plant architecture or leaf morphology. Random individual effects can be spatially autocorrelated, for instance, due to spatial variation in some ecological variables not included in the model. We fitted models with morphology as the response variable and sex as a fixed effect using the R package spaMM version 3.3.0 (Rousset & Ferdy, 2014) in R 4.0.2 (R Core Team, 2020). We fitted models with distinct residual variances for each sex. We updated the code of the spaMM package to allow for different spatial distributions of random effects between sexes (publicly available in spaMM since version 2.6.0). The classical Matérn correlation function was used to model the spatial autocorrelation of random effects, separately for each sex, as a function of the distance between individuals. The Matérn correlation model involves two parameters: the smoothness parameter (ν), and the scaling parameter for distances (ρ). We compared models fits considering (1) no spatial autocorrelation, (2) the same spatial autocorrelation for both sexes and (3) two sex-specific spatial autocorrelation structures. We used likelihood ratio tests (LRTs) to (i) first test for a sex-specific spatial structure in morphology, by comparing models with a different structure of random effects but the same fixed effects, and (ii) to test for sexual dimorphism, by comparing the selected model in step (i) to a model fitted with a similar structure of random effects, but without sex as a fixed effect. Models were fitted either by maximum likelihood, for performing LRTs between models differing in their fixed-effects structure, or by restricted maximum likelihood, for LRTs between models differing in their random-effect structures and for computing the predictions from the model fits. We also examined whether male and female plants were spatially segregated (i.e. whether the sex ratio was spatially autocorrelated) using a binomial generalized linear mixed model (GLMM) with sex as the response variable and compared model fits with or without a spatially autocorrelated random effect as described above. Finally, we tested whether the number of cones produced by females displayed a non-uniform spatial distribution by comparing models including cone number as the response variable, plant density (see below for its calculation) as a fixed factor, and either a spatial random effect or not.

Multivariate sex-specific selection analysis
Inspired by the multivariate framework of Lande and Arnold (1983), we examined in a single full model the relationship between the relative effective fecundity as the response variable and, the following explanatory variables: canopy diameter, leaf area, plant density, sex, and the interaction between each of the three former explanatory variables and sex. We fitted generalized linear mixed- effects models (Gamma GLMM with logarithm link function) to describe the variation in effective fecundity, while including spatially autocorrelated random effects. We checked that a model predicting effective fecundity using a Gamma distribution with a log link performed better than models assuming a different candidate distribution (compare gamma log link: cAIC=18.7 with gaussian with a link identity: cAIC=139; or gaussian with a log link: cAIC=126 -cAIC is a metric similar to the traditional AIC, except that it measures prediction performance conditionally on the realization of the random effects, Vaida & Blanchard, 2005). We followed the same procedure as described for the analysis of sexual dimorphism (see previous section) to fit the models using the R package spaMM and to identify the best fitting random spatial structure (i.e. comparing a sexspecific spatial structure, a common spatial structure and a lack of it). We were thus able to compare the patterns of spatial variation of residual effective fecundity between males and females, once the effects of local density and variation in morphology were included as fixed effects in our global model. Spatial autocorrelation in effective fecundity may then reflect spatial variation in unmeasured environmental variables (such as water availability) or unmeasured plant traits.
Different spatial structures between sexes therefore inform on the differential sensitivities of their effective fecundities to variation in such unmeasured variables. We took into account the variation We chose to investigate the selection gradients for two morphological predictors only (canopy diameter and leaf area), because several of the six morphological traits were highly correlated (Table S3). Notably, all three traits describing plant architecture were strongly correlated with each other (i.e. plant height, canopy diameter and second diameter, We compared the performance of local plant density measured at different scales in explaining variation in male and female effective fecundity (see Supplementary Method S5). We did not transform plant density into z-scores since this predictor is not an individual trait and thus its associated slope should not be considered as a selection gradient. To compare models explaining variation in effective fecundity with different scales used to compute plant density, we used cAIC.
The best model included the density of plants in a quadrat of 12x12m around the focal individual.
We therefore retained this metric to compute the values of plant density used in our all subsequent analyses.
To test for the significance of selection gradients, we first compared the fit of the full model whether there was any effect of the focal predictor. Second, we compared the full model to models in which only the interaction between sex and one of the three variables of interest had been removed. This allowed us to test whether the effect of the predictor was different between sexes. We also built separate GLMMs for each sex, predicting effective fecundity from the three predictors of interest. We used these sex-specific models to test the effect of each predictor in each sex, if and only if, the interaction between sex and a given predictor was significant. We compared the fit of our full model, which considers linear selection gradients only, to the fit of a similar model including quadratic and correlational selection terms for canopy diameter and leaf area, and to model fits including one quadratic or correlational term at a time.
We also estimated selection gradients for the same morphological traits (canopy diameter and leaf area), using as a proxy for fitness, not the estimated effective fecundity but the actual number of cones empirically counted on female plants. We predicted the number of female cones from traits using a Poisson GLMM with spatially autocorrelated random effects and plant density as a covariate in addition to the two focal traits. Finally, we predicted the female effective fecundity from the number of empirically counted cones, also with spatially autocorrelated random effects and plant density as a covariate.

L. rubrum showed strong sexual dimorphism in morphology
Females had significantly smaller canopy diameters (i.e., first and second diameter measures; Figs. 1a, S4c, Table S4), but they were not clearly shorter (Fig. S4a, Table S4). Females had also leaves with a larger area than males (Fig. 1b, Table S4). Correlated measures of leaf morphology (Table   S3) showed similar sexual dimorphism, as females displayed longer and wider leaves than males ( Fig. S4e,f, Table S4).
Spatial structure in the canopy diameter did not clearly differ between sexes, as the fit of a linear mixed model including sex-specific spatially autocorrelated random effects did not produce a likelihood significantly higher than the fit of the model with the same distribution of spatial random effects in both sexes (Table S4). However, a spatially structured random effect considerably improved the fit over the models without it (Table S4; Fig. S5). In contrast, leaf area showed neither a sex-specific spatial structure nor a general spatial structure, (Table S4). Similarly, we found no significant spatial structure for sex ratio (Table S4; Fig. S6).

Dispersal occurred on a smaller spatial scale for seed than for pollen
For both pollen and seed dispersal kernels, our analysis revealed fat-tailed dispersal kernels (i.e. b s and b p lower than one; Fig. 2 and Table 1). Seed and pollen immigration rates were of the same order of magnitude (11% and 15% for seed and pollen respectively). Seed dispersal occurred on a smaller spatial scale than pollen dispersal: the mean estimated dispersal distance of seeds and pollen were respectively of 10.6 and 11,041 meters. These estimated dispersal kernels predict that ca.
distance-dependent pollen dispersal nevertheless explained our data better than modeling a uniform distribution of pollen dispersal distances (Fig. S7, see Supplementary Methods S3 for a description of their comparison). Although we estimated a high probability for the pollen to travel long distance, we estimated that about three quarters of the seedlings had a genotyped father in the population (1m p =0.85 for seedlings with a known mother or υ=0.81 for seedlings with an unknown mother), which is similar to the proportion found for mothers (1-m s =0.89; Table 1). Estimations of dispersal kernels with the NMπ algorithm yielded similar parameter estimation, yet with a notable shorter estimate of mean pollen dispersal distance (Table 1 vs. S2). The estimates of dispersal kernels with NMπ were robust to the inclusion of anisotropy in dispersal events (results not shown).

Autocorrelation in effective fecundity occurred on a smaller spatial scale for males than females
Effective fecundity estimations were carried out using a set of 8 microsatellite markers containing between 6 and 24 alleles, showing non-exclusion probabilities of parent pairs ranging between 0.12 to 0.58 and genotyping error rates ranging from 0.9% to 5.5% (Table S1). The combination of the genotype data and plant spatial distribution data provided information to assign two parents to ~88% of seedlings ( Figure SM3 in Supplementary Methods S3).
We detected a clear sex-specific spatial structure for effective fecundity (Table S5). We found different spatial variation in effective fecundity for each sex, with coarse-grained and finegrained spatial effects for females and males, respectively (Fig. 3a vs. 3b). These results were robust to the removal of plants with standardized effective fecundity greater than four standard errors, which only slightly affected the previous conclusion (Table S5; Fig. S8). Several plants with large effective fecundity were found in the disturbed area (Fig. 3). The number of empirically counted cones on female plants also displayed a significant spatial autocorrelation (χ 2 =2238, df=3, p<0.0001; Fig. S9). Cone number was significantly correlated with relative effective fecundity but the effect size was small (using standardized cone number, β=0.286, χ 2 =6.02, df=1, p=0.0141) and the two spatial distributions were not fully aligned (Fig. 3a vs. S9).

Selection for higher leaf area and wider canopies similar in both sexes
Our spatially corrected selection gradient approach revealed that the leaf area was negatively associated with effective fecundity (Table 2 and S5; Fig. 4), with similar slopes in both sexes as shown by a non-significant interaction between leaf area and sex (Table S5). Yet, the leaf area was positively correlated to the number of counted cones on female plants (β=0.171; LRT: χ 2 =42.7, df=1, p<0.001). Larger canopy diameter was significantly associated with higher effective fecundity (

The effect of plant density on effective fecundity is sex-specific
The effect of plant density on effective fecundity (Fig. S6) differed between sexes, as revealed by a significant interaction between sex and plant density (Table S5; Fig. 4c). Plant density was negatively associated with male effective fecundity, but female effective fecundity showed no association with plant density (Table S5). The effect of the interaction between sex and plant density was however only marginally significant when the plants with the highest effective fecundity were removed in the robustness analysis (Table S5; Fig. S10).

Novel methods for dealing with spatial bias affecting selection estimates in plants
Technical and methodological improvements in parentage assignations now allow for estimation of plant fitness in natural populations from genetic data, and provide the link between fitness and plant traits through selection gradients analyses (e.g. Meagher & Thompson, 1986;Burczyk & Prat, 1997;Burczyk et al., 2006;van Kleunen & Burczyk, 2008). We have developed a methodology that estimates effective fecundity in dioecious plants while accounting for biases associated with their spatial distribution. Beyond the interest of documenting spatial patterns of seed and pollen dispersal, the addition of this spatially explicit component to classical parentage methods (Jones, 2010) improves the estimation of effective fecundity in the presence of confounding effects, such as spatial bias in sampling descendants, spatial variation in the intensity of mate competition triggered by a non-uniform distribution of mates, or border effects (Oddou-Muratorio et al., 2018). We investigated effects of traits on effective fecundity by classical selection gradient methods (Lande & Arnold, 1983), in which we explicitly modeled the effect of spatially correlated unmeasured environmental factors on effective fecundity. This newly developed framework will be particularly suited to estimating selection in natural populations, given that spatial biases are typically difficult to avoid regarding both sampling and uncontrolled factors.

No contemporary sexually antagonistic selection despite strong sexual dimorphism
The signal we found for selection of larger canopy diameters in both sexes may indicate a 'budget effect' of plant size, where larger plants acquire more resources that can be reallocated to gamete production (Delph & Ashman, 2006). In females, the number of counted cones was positively related to canopy diameter. However, cone number only poorly explained effective fecundity. The spatial distribution of cone number and female effective fecundity were moreover not fully matching, suggesting that processes occurring after cone production, such as cone maintenance or mother plant survival affect the female effective fecundity. In males, selection for wider canopies could be linked to flower production if both sexes are subject to similar ontogenetic constraints. We found evidence for similar selection for smaller leaf area in both sexes. Smaller leaves were previously shown to be correlated in L. rubrum with thinner, and more numerous, branches and less efficient water transport from roots to branch apex (Harris, 2007). Selection for smaller leaves in males may reflect selection for greater number of inflorescences held on more flexible branches, a trait long hypothesized to enhance pollen dispersal (Klinkhamer et al., 1997). Smaller leaves may also represent a decreased mechanical hindrance to pollen dispersal. Selection for smaller leaves in females, however, contradicts our expectation of selection for enhanced water transport to the cones.
The number of cones counted on females was furthermore positively related to leaf area. The fact that leaf area relates to the number of cones, and to effective fecundity in opposite ways, suggests that any positive effects of leaf area on fecundity through female cone production may have been masked by trade-offs with other key life history components for serotinous plants (e.g. adult survival until the fire). In the absence of positive genetic correlation between sexes, sexual conflicts may be resolved once each sex reaches its respective optimum. We however did not find evidence for stabilizing selection in the study population, which one could expect if each sex was at its optimum with sufficient genetic variation in the studied traits.
Only 17% of studies estimating selection gradients in animals identified sexually antagonistic selection (Cox & Calsbeek, 2009). The paucity of evidence for sexually antagonistic selection, with which our study concurs, is inconsistent with the idea that genetic correlations between sexes should maintain sexually antagonistic selection over long periods of time (Lande, 1980). To explain this inconsistency, both theory and experiments have suggested that temporal or spatial ecological changes may result in variable pattern of selection acting on males and females, with both sexes displaying trait values remaining far away from their ever changing ecological optimum (Kokko & Rankin, 2006;Delph et al., 2011;Sheridan & Bickford, 2011;Long et al., 2012;Berger et al., 2014;Connallon, 2015;Connallon & Hall, 2016). developments showed that positive inter-sexual covariance for resource acquisition traits could also impede the identification of sexually antagonistic selection (Zafitschek & Connallon, 2017).
An alternative explanation to the lack of detected antagonistic selection is that antagonistic selection is present but masked by a positive correlation between a locally varying unmeasured ecological factor and both effective fecundity and a focal morphological trait (Price et al., 1988). A large body of theoretical work also suggests that the combination of a strong sexual dimorphism with a lack of sexually antagonistic selection found in L. rubrum may result from adaptation to changing ecological conditions, which cause patterns of selection between sexes to align. A longer fire return interval, or low resource availability, are predicted to select for increased resource allocation to plant survival and a lower allocation to cone maintenance (Tonnabel et al., 2012). That  Grant & Grant, 2002;Reimchen & Nosil, 2002;Reed et al., 2013;Acker et al., 2015), including both temporal and spatial variation in the direction of sexually antagonist selection pressure on sexually dimorphic traits (Fargevieille, 2016).

Sex-specific spatial distribution of effective fecundities
The observation that reproductive costs differ between males and females (e.g. Antos & Allen, 1990;McDowell et al., 2000;Harris & Pannell, 2008;van Drunen & Dorken, 2012) has been pivotal to discussions on the evolution of sexual dimorphism in plants (Freeman et al., 1976 Bell, 2008). Resources key to male and female reproduction commonly display small-scale variation in the wild (Silvertown et al., 1999;Araya et al., 2011). Therefore, it is a simple corollary of the sex-specific cost of reproduction hypothesis that male and female fitness should often exhibit different spatial patterns in natural populations, as found in our study population. These sex-specific spatial patterns of fitness variation call for future studies relating small-scale variation in key resource types (Silvertown et al., 1999;Araya et al., 2011) and plant fitness in both males and females.

Only male effective fecundity was affected by density
Male effective fecundity was negatively associated with plant density while no association was found for females. The negative effect of plant density on male effective fecundity might be triggered by increased competition over nutritive resources affecting pollen production, by negative effects of a closed canopy on pollen dispersal, or by competition effects at the seedling stage affecting their offspring. The lack of an effect of plant density on female reproduction suggests that, either they are less affected than males by competition with other plants, or that the negative effects of competition are counter-balanced by positive effects of reproducing in a high density patch. In both cases, it suggests different reproductive needs and ecology in males and females. Plant density was shown to negatively affect both male and female effective fecundity in a wind-pollinated tree (Oddou-Muratorio et al., 2018), and was also shown to impede pollen dispersal in a wind-pollinated herb (Tonnabel et al., 2019b). We note however that this sex-specific effect of density on effective density was only marginally significant when removing individuals with large fecundity.

Pollen and seed dispersal kernels typical of plant dispersal behaviour
Our spatially explicit method allowed the estimation of dispersal kernels, which revealed a fat-tailed seed dispersal kernel in the anemochorous L. rubrum. Most seeds dispersed close to the mother plant, but some fraction dispersed much further. Similarly, a meta-analysis including species from various plant families, continents, vegetation types and growth forms found a predominance of fattailed seed dispersal kernels (Bullock et al., 2017). scarce, but they typically also indicate fat-tailed kernels in both insect- (Austerlitz et al., 2004;Oddou-Muratorio et al., 2005; but see Matter et al., 2013) and wind-pollinated species (Austerlitz et al., 2004;Goto et al., 2006;Gaüzere et al., 2013;Geber et al., 2014 but see Ahee et al., 2015). Our estimated pollen dispersal kernel showed a markedly fat-tailed distribution, whereby a large proportion of pollen was able to disperse over large distances; similar large distance pollen dispersal was reported in both wind-pollinated and animal-pollinated species (e.g. Devaux et al., 2005;O'Connell et al., 2007). Given the large estimates of pollen dispersal distances and the short distance to the nearest population, the low estimates of pollen immigration are unexpected, especially given the same order of magnitude as seed immigration rate. This discrepancy is nevertheless consistent with other studies showing that the amount of long distance dispersal inferred by spatially explicit parentage modeling is not always congruent with the amount predicted by dispersal kernels inferred from local dispersal events (Chybicki & Oleksa, 2018;Hardy et al., 2019). Such inconsistencies may emerge when extrapolation of dispersal kernels does not properly account for an increased probability of encountering obstacles between populations.
In conclusion, we found sex-specific variation in fitness in a natural population of a highly    ). Filled lines correspond to the posterior mean dispersal kernels obtained by averaging parameters of the concatenation of two Markov chains of 50,000 steps (burn-in phase of 10,000 steps). Grey lines illustrate the uncertainty around the averaged dispersal kernel and correspond to the kernels estimated on each iteration of the MCMC. Both dispersal kernels are represented within the minimal and maximal distances existing in our population between females and seedlings for seeds (A.) and between males and females for pollen (B.), i.e., respectively, the maximal female-seedling and male-female distances in our population. The extrapolation of dispersal kernels beyond these limits are not represented in the plots.    Figure S2: Map of the study site showing distribution of females (red circles), males (blue triangles), juveniles (green crosses), non-genotyped females (red full circles) and non-genotyped males (blue full circles). One non-genotyped male that was located nearby the population is not represented in the map. Figure S3: Results of a genetic analysis of the sampled individuals using the software STRUCTURE (Pritchard et al., 2000). Four types of samples were tested: juvenile plants of the third undetermined morphological group and, adult plants of L. salignum and of L. rubrum and juvenile plants of L. rubrum. Four genetic groups were selected by this analysis and are represented here by four different colors (i.e. blue, orange, pink and yellow). Individuals are displayed on two different panels (A. and B.) only for the sake of readability but all individuals belong to the same dataset described in Supplementary Methods S2.   Figure S5: Prediction of canopy diameter in males (A.) and females (B.) as predicted by a linear mixed-effect models including sex as a fixed effect as well as one spatially autocorrelated random effect for each sex. Figure S6: Density for all individuals (A.), for males (B.) and for females (C.) as numbers of individuals in quadrats of 12x12 meters. Figure S7: Comparisons of conditional likelihood between the model presented in the main text accounting for spatial structure of pollen dispersal (in green) and a model neglecting such a spatial structure of pollen dispersal (in brown) as a function of Bayesian iterations which were run for 50,000 steps. These simulations were run using uniform prior distributions with the following intervals [0  Figure S8: Replicate of Fig. 3 after the removal of individuals with effective fecundity greater than four standard errors (i.e. five females and one male). See Fig. 3 for legend details. Red crosses correspond to the locations of plants removed in the sensitivity analysis. Figure S9: Spatial prediction of the number of empirically counted cones as predicted by a generalized linear mixed-effect model including all fixed effects (canopy diameter, leaf area, plant density) as well as one spatially autocorrelated random effect. Figure S10: Replicate of Fig. 4 after the removal of individuals with effective fecundity greater than four standard errors (i.e. five females and one male). See Fig. 4