An analysis of mating biases in trees.

Assortative mating is a deviation from random mating based on phenotypic similarity. As it is much better studied in animals than in plants, we investigate for trees whether kinship of realized mating pairs deviates from what is expected from the set of potential mates and use this information to infer mating biases that may result from kin recognition and/or assortative mating. Our analysis covers 20 species of trees for which microsatellite data is available for adult populations (potential mates) as well as seed arrays. We test whether mean relatedness of observed mating pairs deviates from null expectations that only take pollen dispersal distances into account (estimated from the same data set). This allows the identification of elevated as well as reduced kinship among realized mating pairs, indicative of positive and negative assortative mating, respectively. The test is also able to distinguish elevated biparental inbreeding that occurs solely as a result of related pairs growing closer to each other from further assortativeness. Assortative mating in trees appears potentially common but not ubiquitous: nine data sets show mating bias with elevated inbreeding, nine do not deviate significantly from the null expectation, and two show mating bias with reduced inbreeding. While our data sets lack direct information on phenology, our investigation of the phenological literature for each species identifies flowering phenology as a potential driver of positive assortative mating (leading to elevated inbreeding) in trees. Since active kin recognition provides an alternative hypothesis for these patterns, we encourage further investigations on the processes and traits that influence mating patterns in trees.

INTRODUCTION

Systematic deviations from random mating have major ecological and evolutionary implications for example, sexual selection (Jiang, Bolnick, & Kirkpatrick, 2013), reinforcement (Kirkpatrick, 2000; Servedio & Noor, 2003) and inbreeding (Epinat & Lenormand, 2009). A specific type of deviation from random mating is assortative mating where the probability of mating is influenced by phenotype. In large well‐mixed populations in mutation‐drift balance, assortative mating can occur without necessarily elevating the relatedness of mating pairs at all (Templeton, 2006). However, in spatially structured populations, individuals that possess identical alleles at loci that determine assortativeness may often do so because of genetic relatedness. Matings between phenotypically similar individuals may in these cases imply elevated biparental inbreeding (in which case the trait similarity is caused by alleles identical by descent). Positive assortative mating increases homozygosity, linkage disequilibrium and the variance of quantitative traits (Weis, 2005; Wright, 1921). It can also contribute to reproductive isolation between diverging populations (Bolnick & Kirkpatrick, 2012; Gavrilets, 2003; Howard, 1999). Negative assortative mating has the potential to perform the opposite, that is, to decrease homozygosity and inbreeding depression (Pusey & Wolf, 1996; Waser, 1993).

Assortative mating is well investigated across a wide array of animal taxa, where positive assortment has been shown to predominate (Jiang et al., 2013). Plants appear less well studied in this respect, although spatial assortment patterns that are produced by sessile life (impacting “who mates with whom” via biotic or abiotic pollen transfer) may combine with traits such as inbreeding avoidance and self‐incompatibility mechanisms (Bedinger, Broz, Tovar‐Mendez, & McClure, 2017; de Nettancourt, 2001; Takayama & Isogai, 2005) to produce a variety of potential negative assortment patterns. Several additional mechanisms have been suggested to promote positive assortative mating, including assortative pollination between synchronously flowering plants, repeatable pollinator behaviour favouring mating among plants of similar size or flower traits, as well as assortative fertilization after successful pollen transfer (Waser, 1993). A recent simulation study shows that synchrony of flowering time is capable of structuring populations in a manner that is identifiable at neutral loci; these show much stronger differentiation when, in addition to the commonly assumed isolation‐by‐distance, matings are nonrandom due to flowering phenology having to match between the pollen donor and the pollen recipient (Peters & Weis, 2019).

After pollen has reached the stigma, the sequence from pollen adhesion, pollen germination and pollen tube growth to fertilization features complex cellular interactions (Hiscock & Allen, 2008). Pollen that do not adhere to the stigma and fall off cannot reach the germination state. Pollen adherence as a first stage of pollen sorting is relevant at least in plants with dry stigmas (Hiscock & Allen, 2008; Luu, Heizmann, & Dumas, 1997). Timing is important in this context: plants need to have similar phenologies for fertilization to be possible at all, and the importance of timing extends to finer scales. Prior to pollen germination, delayed stigma receptivity can synchronize the fertilization chances of early and late arriving pollen (Herrero & Hormaza, 1996). Timing of stigma receptivity is thought to be a female strategy to increase the number of pollen donors and male‐male competition (Herrero & Hormaza, 1996; Madjidian, Hydbom, & Lankinen, 2012).

Polyandry is considered to be beneficial to female fitness by promoting pollen competition and alleviating pollen limitation (Pannell & Labouche, 2013). Competition among pollen is largely mediated via speed of pollen tube growth (Skogsmyr & Lankinen, 1999; Snow & Spira, 1991; Swanson, Hammond, Carlson, Gong, & Donovan, 2016), and it has been argued to mitigate effects of inbreeding, should related pollen land on a stigma (see Armbruster & Rogers, 2004 for a selfing context). Pollen competition is strongly influenced by female functional traits, such as the length of the pistil (Lankinen & Skogsmyr, 2001; Ramesha et al., 2011), the provisioning of limited nutrients for heterotrophic pollen germination and pollen tube growth (Herrero & Hormaza, 1996) as well as guidance of pollen tube growth with chemical cues (Hiscock & Allen, 2008). Studies of pollen competition typically apply different mixtures of self and outcross pollen (e.g., Lankinen, Smith, Andersson, & Madjidian, 2016). Although varying siring rates among mating pairs in pollen mixtures have been detected in an outcrossing context (Marshall, 1998; Snow & Spira, 1996), the effects of potential relatedness between mates in the absence of selfing appear understudied. Given what is known about the complex pollen‐pistil interactions, it seems likely that the lack of a pistil in gymnosperms reduces the effectivity of pollen sorting or pollen competition.

Even if a plant is an obligate outcrosser (self‐incompatible), one can expect relatedness to play a role in fertilization. Pollen and seed dispersal often show moderate to strong spatial restriction, which results in spatial genetic structuring (SGS) of the adult population and more frequent mating among nearby relatives (Ghazoul, 2005; Vekemans & Hardy, 2004). As mentioned above (Peters & Weis, 2019), temporal variation can be equally important in structuring populations. If variation in flowering time is heritable, it can increase mating among related individuals (Daïnou et al., 2012; Ennos & Dodson, 1987; Suni & Whiteley, 2015), forming a clear causal route from assortativeness (with respect to timing traits) to relatedness among mating pairs—a route that is absent in the simplest models (e.g., Templeton, 2006) that, for conceptual clarity, assume that populations are so large and well mixed that assortativeness can operate without making relatives meet. Yet another real‐life complication is that SGS and variability in flowering time are interrelated; while flowering synchrony and spatial proximity both increase mating frequency (Ison, Wagenius, Reitz, & Ashley, 2014), variation in flowering time increases pollen dispersal and reduces SGS (Kitamoto et al., 2006).

Interestingly, studies on the effect of distance between mates on siring success and offspring vigor frequently have found there to be an optimal outcrossing distance (Schierup & Christiansen, 1996; Souto, Aizen, & Premoli, 2002; Waser & Price, 1991; Waser, Price, & Shaw, 2000). Together with the evidence presented above that mating is not necessarily random because of phenological aspects as well as competition among pollen that have landed on a stigma, these patterns raise the question of relatedness levels that plants “choose” in terms of realized fertilizations, relative to what is available to them in terms of pollen dispersal distances. As the spatial arrangement of adults and the pollen dispersal distances are accounted for, “choice” (or more generally, mating bias, Kokko, Brooks, Jennions, & Morley, 2003) has to be driven by traits of the mating pairs. These can be based on some form of “active” choice that is visible in pollen‐pistil interactions or early abortion of seeds; or they may be more classifiable as mating bias, reflecting similar phenologies of related individuals or any floral trait that pollinators make a positive association with.

Here we perform an analysis of mating biases utilizing published data on trees. We do not seek to distinguish between active choice (based on some form of kin recognition) and more general mating biases; instead, we ask whether the realized mating pattern gives overall evidence for kin‐based choices, assortative mating with respect to a trait shared by both the pollen donor and the pollen recipient, or both processes acting together.

Although the tree growth habit does not represent a monophyletic group, it is an interesting growth form with implications for reproductive strategies. Trees are known to show high genetic diversity, maintained in populations that feature strong inbreeding depression, extremely high recruit mortality and high pollen flow which reduces local adaptation while strong selection can promote local adaptation (Petit & Hampe, 2006). In trees such local adaptation by strong selection is most common at scales of tens to hundreds of kilometers (Savolainen, Pyhäjärvi, & Knürr, 2007). While pollen‐pistil interactions in trees are understudied, it appears likely that trees share the pollen precedence patterns that are commonly observed in various angiosperms. Post‐fertilization processes, important in shaping the realized mating patterns, seem better studied. In particular, in a combination of controlled and open pollination experiments, trees have shown high post‐zygotic abortion rates of premature seed (Bawa & Webb, 1984), which is thought to reflect either early acting inbreeding depression or maternal selection (Seavey & Bawa, 1986).

MATERIALS AND METHODS

Our analysis is inspired by an approach recently published by Monthé, Hardy, Doucet, Loo, and Duminil (2017), who designed a method to test whether mating events between related individuals occurred more or less often than would be expected by chance, where “chance” is estimated against the background of the spatial extent of pollen dispersal. This extent of pollen dispersal is estimated based on a sample of realized mating events irrespective of relatedness. One can then ask if there is an elevated or depressed probability for a fertilization to happen, for a given distance, if the parents are related. Effectively, this is achieved by resampling a distribution of matings assuming that relatedness has no effect (while distance between potential mates retains its observed effect) and creating a distribution of pairwise kinship coefficient (F; Loiselle, Sork, Nason, & Graham, 1995) values based on this set of assumptions. Each resampling procedure creates a sample of F values for which the mean is recorded; repeating the procedure a large number of times creates the distribution for mean F that can be expected under the null hypothesis that distance matters while relatedness does not. The observed mean F in the real data set can be compared to this expectation (see below for statistical details), to see if realized offspring are a result of more, or less, related mating pairs than expected by chance. Deviations of the observed from the expected mean F are then discussed in terms of mating biases that may result from assortative mating and/or kin recognition.

We compiled microsatellite genotype data of adult trees with physical coordinates and of seed arrays from published studies. We used 20 published data sets for which paternity analysis of seed arrays was available; our data sets were found with the following procedure. We searched once for “pollen dispersal” and once for “paternity analysis” on Dryad Digital Repository (http://datadryad.org), which returned 96 and 147 results, respectively (retrieved on 8 December 2017). These results were screened for genotype data which included at least six nuclear microsatellite loci from seed arrays of putative or known mothers, with the additional criterion that the coordinates for adult trees had to be known (note that in five out of the 20 data sets the seed were collected from the ground and were tested for genotype matching with the putative mother). Data sets with <50 adult trees, <5 seed arrays, or <200 seed genotypes were excluded. In a further stage of analysis, we excluded data sets for which paternity analysis (see below) resulted in <20 outcrossing events.

We also included a study on Jacaranda copaia despite it narrowly missing one of our criteria (being genotyped only at five loci), because we considered a total of 98 alleles in the adult population to contain sufficient information content. We found two additional suitable data sets in the Supporting Information from articles published in PLoS One, and upon request, we received additional six data sets directly from authors (see Table 2 and Acknowledgements).

Table 2

Overview of investigated tree species, data used and assessment of assortative mating in trees. F: Kinship coefficient (Loiselle et al., 1995). The column “Location within Table 1” indicates potentially positive assortative mating (Table 1, Column A), no indication of assortative mating (Table 1, Column B) and potentially negative assortative mating (Table 1, Column C) as well as biparental inbreeding (Table 1, Row a) and no biparental inbreeding (Table 1, Row b)

Species	Family	Breeding system	No. of trees	No. of seed	No. of assigned outcrossing events	Selfing rate (%)	Observed mean Fij	Predicted mean Fij	Percentile observed mean Fij of predicted mean Fij	Percentile zero Fij of predicted mean Fij	Location within Table 1		Reference	Data source	Data‐Subset used
											Column	Row
Glionnetia sericea	Rubiaceae	Monoecious	210	622	122	8	0.1836	0.0252	Above 100%	1%	A	a	Finger et al. (2014)	Dryad
Quercus robur	Fagaceae	Monoecious	135	684	603	7	0.1286	0.0615	Above 100%	Below 0%	A	a	Moracho et al. (2016)	Dryad
Dysoxylum malabaricum	Meliaceae	Monoecious	235	566	508	1	0.0743	0.0257	Above 100%	Below 0%	A	a	Ismail et al. (2012)	Dryad
Vateria indica	Dipterocarpaceae	Monoecious	85	259	132	13	0.0633	0.0208	Above 100%	2%	A	a	Ismail et al. (2014)	Authors
Cariniana legalis	Lecythidaceae	Monoecious	65	600	481	8	0.1346	0.0044	Above 100%	18%	A	b	Tambarussi, Boshier, Vencovsky, Freitas, and Sebbenn (2015)	Authors
Entandrophragma cylindricum	Meliaceae	Monoecious	239	484	303	10	0.0586	0.0062	Above 100%	17%	A	b	Monthé et al. (2017)	Dryad
Erythrophleum suaveolens	Fabaceae	Monoecious	88	239	82	9	0.0317	0.0000	99%	50%	A	b	Duminil et al. (2015)	Dryad	DRC
Quercus petraea	Fagaceae	Monoecious	162	545	352	2	0.0170	0.0046	100%	14%	A	b	Gerber et al. (2014)	PlosOne	France
Oenocarpus bataua	Arecaceae	Monoecious	185	318	256	1	0.0120	0.0004	98%	48%	A	b	Ottewell et al. (2012)	Dryad
Fraxinus excelsior	Oleaceae	Heterodichogamous	219	500	304	0	0.0162	0.0193	30%	0%	B	a	Semizer‐Cuming, Kjær, and Finkeldey, (2017)	PlosOne
Acer pictum	Aceraceae	Heterodichogamous	97	1,041	593	8	0.0088	0.0045	84%	2%	B	a	Shang et al. (2012)	Dryad
Castanopsis sieboldii	Fagaceae	Monoecious	145	486	355	2	0.0251	0.0186	93%	Below 0%	B	a	Nakanishi et al. (2012)	Dryad
Sorbus domestica	Rosaceae	Monoecious	164	1,101	325	13	0.0360	0.0260	92%	Below 0%	B	a	Kamm et al. (2009)	Authors
Baillonella toxisperma	Sapotaceae	Monoecious	87	230	50	8	0.0357	0.0126	88%	27%	B	b	Duminil et al. (2016)	Authors
Phoenix canariensis	Arecaceae	Dioecious	103	616	173	0	0.0028	0.0018	58%	35%	B	b	Saro et al. (2014)	Dryad
Shorea xanthophylla	Dipterocarpaceae	Monoecious	170	456	86	2	0.0011	0.0278	9%	8%	B	b	Kettle et al. (2011)	Authors
Jacaranda copaia	Bignoniaceae	Monoecious	207	285	161	1	0.0009	0.0025	41%	36%	B	b	Vinson et al. (2015)	Dryad
Abies pinsapo	Pinaceae	Monoecious	162	598	262	29	0.0006	0.0002	9%	47%	B	b	Sánchez‐Robles et al. (2014)	Dryad
Prunus lannesiana	Rosaceae	Monoecious	78	221	121	0	−0.0201	0.0052	2%	33%	C	b	Shuri et al. (2012)	Dryad	B
Bertholletia excelsa	Lecythidaceae	Monoecious	134	338	63	0	−0.0326	0.0146	0%	18%	C	b	Baldoni et al. (2017)	Authors

If a study contained data on multiple populations that met our data requirements, we included only the population for which we achieved the highest assignment rate in the paternity analysis. This resulted in a total of 19 data sets of angiosperm and of one gymnosperm tree species (12 from Dryad, two from PLoS One and six directly from authors).

Common among the selected studies is the general aim to characterise pollen dispersal and the aim to sample at least in a core area exhaustively all adult individuals. In 14 of the studies the motivation is to investigate mating patterns such as levels of inbreeding in small isolated, fragmented or selectively logged populations. The remaining six studies were conducted in long‐term undisturbed natural populations. See Table S1 for a species‐specific classification of the disturbance.

For all data sets we calculated pairwise kinship coefficients (F, Loiselle et al., 1995) and pairwise geographic distances for all adult trees, using software SPAGeDi 1–5a (Hardy & Vekemans, 2002). SPAGeDi software has been designed to characterize the association between genetic and spatial distances based on genotype data of mapped individuals. The most widely used metric of individual genetic distance implemented in SPAGeDi is Loiselle's kinship coefficients (F) which quantifies the correlation in the frequencies of homologous alleles, p and p, in pairs of individuals i and j. This kinship coefficient is estimated aswhere the first term denotes the expected value of F, with k = n(n–1)/2 being the number of possible pairwise distances between n individuals. The second term corrects for bias due to limited sample size and ensures that F has an expected value of zero for a population in Hardy–Weinberg equilibrium. This kinship coefficient is relative to an arbitrary sample which implies that negative kinship coefficients naturally occur between some individuals that are less related than random individuals. In addition, we computed with SPAGeDi the mean regression slope of the kinship coefficient over the natural logarithm of the distance and the mean kinship of the adult pairs within the first distance class . From that, we calculated the Sp‐statistic which is a comparable measure of SGS intensity defined as (Vekemans & Hardy, 2004).

Thereafter we applied a likelihood approach for paternity assignments with cervus 3.0.7 (Kalinowski, Taper, & Marshall, 2007; Marshall, Slate, Kruuk, & Pemberton, 1998). The definition of the 95% confidence level of the DELTA score, that is, the difference between logarithm of odds (LOD) score of the most likely parent and the second most likely parent, was based on a simulation of paternity using the following settings: 10,000 offspring were simulated, based on adult allele frequencies, with all adult trees serving as candidate fathers (allowing for selfing). The proportion of sampled fathers was set to 0.95, the proportion of loci typed was based on the observed value of the adult genotypes, the proportion of mistyped loci was set to 0.01, the minimum number of loci genotyped was set to six (except for Quercus robur genotyped at six loci for most but not all individuals, and Jacaranda copaia genotyped at five loci, we set this value to five). Unassigned seed are assumed to represent pollen dispersal from unsampled pollen sources outside the study areas. For one species (Shorea xanthophylla) we received only the cervus output and not the raw offspring genotype table.

From the SPAGeDi output we then extracted the pairwise distances and kinship coefficients (F) between the known mothers and the assigned fathers (excluding selfing events). Based on the pairwise distances of these observed mating events, we created a histogram of effective pollen dispersal distances (separately for each species under its own study‐specific spatial arrangement of adults) with a target value of 12 breaks. This was done in r 3.4.1 (R Core Team, 2017); note that the r algorithm can result in different numbers of break points. Setting the number of breaks higher than 12 did not lead to a qualitative change of the results (data not shown) while it increased the instances of distance classes to not contain any adults, especially if the species in question had a fragmented distribution.

The density distribution of the histogram was then used as the study‐specific probability distribution for pollen dispersal distances. Using this distribution, we sampled adults as potential mating pairs at similar distances as the observed mating pairs for each study‐specific spatial genetic structure and arrangement of adult trees in a given area. This sampling procedure makes pairs form at equivalent distances as in real life, but irrespective of any information on relatedness. Using a discrete probability distribution of pollen dispersal instead of a continuous distribution ensures the existence of alternative potential mating pairs within a specific distance class of an observed mating event (a precondition for the resampling to work). We then calculated the mean F for these distance‐predicted mating pairs. For obtaining the possible range of predicted mean F, we repeated the sampling 10,000 times. The observed mean F was considered to deviate significantly from the null hypothesis when it was above the 97.5 percentile or below the 2.5 percentile of the predicted mean F values.

It was also possible that the observed F value is within the 2.5–97.5 percentile range with a largely positive distribution of predicted values for F. If this occurred, with the “largely positive” quantified with the criterion that the value of F = 0 lies below the fifth percentile of the predicted mean F values (i.e., at least 95% of the predicted F values were positive), then we considered this to constitute evidence for elevated biparental inbreeding caused simply by the SGS of the population together with the pollen dispersal patterns.

We also considered that matings may not be independent within one fruit, and repeated the sampling considering only unique mating events.

Compared to the original approach described by Monthé et al. (2017), our procedure samples potential mating pairs across the whole population, without a focus on specific focal mother trees for which one is deriving potential mates. Accordingly, we also deviate from their approach in which they sampled potential fathers at a distance from the mother 0.9 and 1.1 times the observed pollen dispersal distance; we considered all observed distances to be possible in principle for any focal individual.

Although the observed pairwise F values show a continuous distribution around a peak often close to zero, they are not necessarily normally distributed (see Supplementary Information for an overview of the pairwise F values of all observed mating pairs, Figure S1, and for the unique mating pairs, Figure S2, for each species). However, for all species, the mean of this distribution is close to the modal value, that is, the most frequent kinship, of realized mating pairs.

RESULTS

The paternity assignment rate ranged from 18.6% to 90.3% while the selfing rate ranged from zero to 24.6% (assuming that the unassigned seed are from outcrossing events). The lowest exclusion probability of parent pairs was 0.998 for the species B. toxisperma. Excluding the selfing events, this resulted in 50–603 mating pairs of which between 16 and 410 were unique mating events. Because it is common to include multiple matings for direct estimates of pollen dispersal distances (Austerlitz et al., 2004), our presentation focuses on the results on all mating pairs, but we provide the results for both approaches in the Supporting Information (Table S1). The results based on unique mating pairs are qualitatively similar; note that the majority of species are one seeded or only one seed per fruit was genotyped, which increases independence of the mating events.

The detected mating patterns can be classified according to the observed mean F relative to the distribution of the predicted mean F values into three cases (Table 1, Columns A, B and C): (a) Observed mean F is high relative to the predicted values; (b) no deviations of the observed mean F from the predicted values; and (c) observed mean F is low relative to the predicted values.

Table 1

Mating pattern classification with a graphical illustration for sorting the results

		A	B	C
		Observed mean Fij above the 97.5 percentile of predicted mean Fij values	Observed mean Fij between the 2.5 and 97.5 percentile of predicted mean Fij values	Observed mean Fij below the 2.5 percentile of predicted mean Fij values
a	Fifth percentile of predicted mean Fij values is above zero	Increased inbreeding potentially indicating positive assortative mating Inbreeding caused by restricted pollen dispersal	No increased inbreeding beyond pollen dispersal based predictions Inbreeding caused by restricted pollen dispersal	Reduced inbreeding potentially indicating negative assortative mating Inbreeding caused by restricted pollen dispersal
b	Fifth percentile of predicted mean Fij values is below zero	Increased inbreeding potentially indicating positive assortative mating No inbreeding caused by pollen dispersal pattern	No increased inbreeding beyond pollen dispersal based predictions No inbreeding caused by pollen dispersal patterns	Reduced inbreeding potentially indicating negative assortative mating No inbreeding caused by pollen dispersal patterns

Columns indicate whether the observed mean F is above (A), within (B) or below (C) the pollen dispersal predicted value of F. Rows indicate whether SGS and pollen dispersal distance leads to elevated bi‐parental inbreeding (a) or not (b). None of the investigated species met the conditions for being classified under Column C, Row a. The exemplary histograms depict two possible distributions of predicted mean F values (Rows a and b). The dashed grey lines indicate the upper 97.5 percentile or the lower 2.5 percentile and the solid grey lines indicate the fifth percentile. The triangle and diamond icons indicate different exemplary positions of an observed mean F value relative to the predicted mean values. The different icon formats represent the mating pattern classifications and correspond with their use in Figure 1.

Figure 1

Boxplots of predicted values of mean F with median, upper and lower quartiles and observed mean F indicated as diamond or triangle icons. The whiskers represent the full range of the predicted mean F values. The vertical dashed line at zero indicates the expected mean F under random mating. Reverse triangles indicate significantly elevated observed mean F relative to the predicted mean F values indicating positive assortative mating (Table 1, Column A). Diamonds indicate no significant deviation of the observed mean F relative to the predicted mean F values (Table 1, Column B). Upward triangles indicate significantly lower observed mean F relative to the predicted values indicating negative assortative mating (Table 1, Column C). Black icons indicate biparental inbreeding determined by pollen dispersal distances with significantly elevated predicted mean F relative to zero F (Table 1, Row a). Grey icons indicate no significant deviation of the predicted mean F values from zero and consequently no biparental inbreeding induced by pollen dispersal pattern (Table 1, Row b)

Additionally, we classify whether the predicted mean F is largely positive or not (Table 1, Rows a and b): (a) Predicted mean F is largely positive (high relative to zero F) and (b) predicted mean F does not deviate from zero.

Species classified under Case A exhibit more frequent mating among related individuals than predicted by the observed pollen dispersal pattern. Although this can result from either kin preferences or actual assortative mating, we call this “positive assortative mating” for simplicity (see Section 4 for our views of the biological plausibility of assortativeness vs. kin recognition). Under Case B the observed F follows what is expected under the observed pollen dispersal, with no additional information about assortativeness and/or kin recognition. Species meeting the condition of Case C have less frequent mating among relatives than predicted by the observed pollen dispersal pattern, implying negative assortative mating and/or active inbreeding avoidance (again, we use “negative assortative mating” below for simplicity). The cases a and b indicate whether the observed pollen dispersal distances, together with the spatial arrangement of adult trees, predict biparental inbreeding to occur in the first place.

Nine species have an observed mean F among mating pairs above the 97.5 percentile of the distance‐predicted mean F values, where elevated inbreeding suggests some form of positive assortative mating (Table 1 Column A). Of these nine species, six have an observed mean F value which is above the entire range of the predicted values. Furthermore, in four out of these nine species, the F of zero is below the fifth percentile of the predicted values (Table 1, Column A, Row a). These four cases exhibit in addition to a pollen dispersal independent effect also a component of biparental inbreeding driven by pollen dispersal distance and adult SGS. Four species only show a deviation of the predicted mean F values from zero but without any substantial deviation from the observed mean F (Table 1, Column B, Row a). These four species experience biparental inbreeding determined by pollen dispersal distance and adult SGS but with no genome‐wide effect of potential positive assortative mating. Five species show no deviations of the observed mean F from the predicted values and the predicted values do not deviate from zero F (Table 1, Column B, Row b). In two species the observed mean F lies below the 2.5 percentile of the predicted values which do not deviate from zero F (Table 1 Column C, Row b). These two species indicate that outbreeding that associates with negative assortative mating can also occur. No species was detected to have negative assortative mating together with elevated biparental inbreeding (Table 1, Column C, Row a).

The Sp‐statistic revealed that all except three species show significant SGS (Table S1). These are the two species which show reduced kinship driven by negative assortative mating (Table 1 Column C, Row b) and one species (Abies pinsapo) which shows no inbreeding caused by pollen dispersal patterns but a tendency for negative assortative mating with 91% of the predicted mean kinship coefficients higher than the observed mean kinship. It is noteworthy that the observed and expected F values correlate (R = .482, p = .0316, two‐sided test). This correlation is reassuring, as it indicates that spatial arrangement of the adult trees and the pollen dispersal pattern (factors used to derive the expected F) appear to have effects in real life (the observed F).

Table 2 provides a species‐specific overview of the results of the paternity analysis, of the assessment of assortative mating and of the biparental inbreeding. Figure 1 provides a graphical illustration of the observed mean F in relation to the predicted mean F conditional of the pollen dispersal distance.

DISCUSSION

Our investigation demonstrates that, in trees, mating among related individuals is commonly more frequent than what would be expected by the observed pollen dispersal distances. There are at least two, nonmutually exclusive, reasons why this might be the case. First, it could be that trees actively favour related pollen to fertilize their ovules, either prior to fertilization (in the pollen‐pistil interactions) or afterwards (via abortion of seeds). We cannot exclude this option for any of the cases we have labelled as “positive assortative mating” (see Fisher, 1941; Kokko & Ots, 2006 for why it may even be adaptive to favour related mates despite some inbreeding depression). Simultaneously, it is worth reflecting whether apparent positive assortative mating is indicative of actual positive assortative mating, a form of nonrandom mating where some phenotypic trait covaries with genome‐wide genetic similarity and simultaneously impacts the likelihood that a mating pair forms. While our method does not directly comment on which phenotypic trait(s) might be responsible, it allows us to distinguish whether biparental inbreeding simply occurs at the rate predicted by the observed pollen dispersal distances (related trees are more likely to mate due to spatial genetic structure), or if an additional mechanism of positive assortative mating is at work.

The most widely investigated heritable trait capable of driving assortative mating in plants is variation in flowering phenology (Ennos & Dodson, 1987; Fox, 2003; Weis, Nardone, & Fox, 2014). While our data set does not include information on phenology, it appears useful to list what is known about phenological variation for each of the species we consider. For seven out of the nine species which show pairwise F elevated beyond distance‐based predictions (Table 1, Column A) we found direct or indirect evidence for variation in flowering time within populations with the potential to contribute to the observed assortative mating, assuming that timing has a heritable component (leading to an expectation of similar phenology in related individuals). In particular Quercus robur (Moracho, Moreno, Hampe, & Jordano, 2016), Dysoxylum malabaricum (S. A. Ismail, unpublished data), Oenocarpus bataua (Ottewell, Grey, Castillo, & Karubian, 2012) and Quercus petraea (Bacilieri, Ducousso, & Kremer, 1995) have pronounced temporal flowering variation within populations. For three species the evidence for within population flowering variation is more indirect: Populations of Erythrophleum suaveolens have a flowering season of around 4 months which clearly exceeds individual flowering period of around 2 months (Duminil et al., 2015). In this species elevated inbreeding due to positive assortative mating has recently been confirmed in a population from Cameroon (Hardy et al., 2019). Glionettia sericea and Cariniana legalis have a very long flowering season with a predominant flowering over 7 months (Finger, Kaiser‐Bunbury, Kettle, Valentin, & Ghazoul, 2014), or asynchronous supra annual flowering with a flowering season of around 4 months respectively (Engel & Martins, 2005); while the evidence is indirect, this is in line with the assumption that there may be heritable variation in flowering phenology. Currently available information on flowering phenology of Entandrophragma cylindricum is ambiguous: monthly observations over 3 years suggest that the species flowers largely synchronously (Fétéké et al., 2016), but flowering among individuals and years has also been stated as irregular (Lourmas, Kjellberg, Dessard, Joly, & Chevallier, 2007). The possibly most synchronously flowering species for which we detected positive assortative mating is Vateria indica. In this biannually flowering species all individuals within a population start flowering within 1 week for around 2 months, while over large scales even the flowering years can alternate (S. A. Ismail, personal observation 2009 to 2012).

Nine species show no deviation of mating pair relatedness from distance‐predicted expectations (Table 1, Column B). Four of these species show asynchronous flowering, while the remaining five flower rather synchronously. Variable flowering phenology has been reported for Acer pictum (Shang, Luo, & Bai, 2012), Sorbus domestica (Kamm et al., 2009) and Phoenix canariensis (Saro, Robledo‐Arnuncio, González‐Pérez, & Sosa, 2014). For Baillonella toxisperma we infer variation in flowering time because the flowering season of 2 months clearly exceeds the individual flowering time of around 1 month (Duminil et al., 2016). Synchronous flowering has been reported for Castanopsis sieboldii (Yumoto, 1987), Shorea xanthophylla (Sakai et al., 2005), Jacaranda copaia (Vinson, Kanashiro, Harris, & Boshier, 2015) and Fraxinus excelsior (Gérard, Fernandez‐Manjarres, & Frascaria‐Lacoste, 2006). In these four species, the available information on “flowering synchrony” is a generic qualitative statement and the temporal resolution of this classification remains vague. A special case is Abies pinsapo where fertilization is synchronized by delaying pollen germination up to more than a month after pollination (Arista & Talavera, 1994) which we interpret as functional flowering synchrony. Although much of the information above remains qualitative, it fits the general expectation provided by the model of Peters and Weis (2019): temporal structuring of populations can leave signatures even at neutral loci. To see if a more specific prediction is true – differences between loci in F depending on whether they are linked to the locus causing the phenotype that underlies assortativeness (Peters & Weis, 2019) – would be interesting to test, however, none of the cases above yield direct information on the loci that impact flowering phenology (which may also be a polygenic trait).

Both species, Prunus lannesiana and Bertholletia excelsa, which show negative assortative mating (Table 1, Column D) have been reported to exhibit variation in flowering phenology: In P. lannesiana synchrony in flowering among mating pairs was found to be a significant factor for predicting mating pair fecundity (Shuri et al., 2012). In Bertholletia excelsa the main flowering occurs over a period of 4 months while the individual flowering period is only three to eight weeks (Maués, 2002).

Interestingly P. lannesiana and B. excelsa, which show elevated outbreeding, potentially indicative of negative assortative mating, are self‐incompatible (Maués, 2002; Shuri et al., 2012). Given that both species show variation in flowering time within populations, one might have predicted the opposite (positive assortativeness). Negative assortative mating might indicate a combination of weak heritability in phenological traits, making it easier for any a priori positive assortativeness expectation to be overridden by self‐incompatibility mechanisms that also reduce mating among close relatives (de Nettancourt, 2001), or by other inbreeding avoidance mechanisms such as strong post‐zygotic selection against inbred offspring (Ghazoul & Satake, 2009). Negative assortative mating may also be adaptive under stabilizing selection because it reduces the production of less fit phenotypic extremes (Kondrashov & Shpak, 1998). Whether populations with elevated outbreeding due to negative assortative mating indicate effective inbreeding avoidance and/or stabilizing selection remains to be studied in more detail.

Finally, there are two special cases, the heterodichogamous species A. pictum and F. excelsior, where some individuals produce first either male or female gametes and then switch their sexual function in the course of the flowering season. The resulting more frequent mating among reciprocal sexual morphs is commonly presented as a case of disassortative mating mediated by phenotypic dissimilarity (Bai, Zeng, & Zhang, 2007; Gleiser, Chybicki, González‐Martínez, & Aizen, 2018; Gleiser, Verdú, Segarra‐Moragues, González‐Martínez, & Pannell, 2008; Shang et al., 2012). When considering genetic similarity of mating pairs at neutral loci, we find no indication of either negative or positive assortative mating for these two species, while both show spatially determined biparental inbreeding. This indicates that in these two species the loci determining the dichotomous sexual morphs do not reflect the genetic dissimilarity at the neutral loci investigated in this study.

As a whole, our results support the view that asynchronous flowering within populations, or any other similarly functioning trait, can yield positively assortative realized mating patterns, which in turn can increase levels of inbreeding detectable at neutral loci. This is in line with theoretical expectations showing that variation in reproductive phenology causes isolation by time with qualitatively similar effects as isolation by distance (Hendry & Day, 2005; Peters & Weis, 2019). Such temporal structure at neutral loci has been detected in natural populations (Daïnou et al., 2012; Suni & Whiteley, 2015). If spatial genetic structure combines with spatially structured phenotypic similarity, and the phenotypes impact “who mates with whom” (mating biases), trait similarity of mating pairs will yield matings between related individuals. As suggested for the herbaceous Brassica rapa, spatial clustering of flowering phenology can amplify the spatial aggregation of genetic similarity (Ison & Weis, 2017), and together with restricted pollen dispersal, elevates biparental inbreeding.

However, out of the nine species with no elevated inbreeding beyond expected values (Table 1, Column B), six show variation in flowering time, which indicates that variable flowering phenology does not inevitably result in more frequent mating of more related individuals. In addition, even if the population appears to flower asynchronously, there may be no phenologically induced assortative mating because the flowering schedule is determined by an array of components such as flowering date, peak flowering and duration of flowering which do not straightforwardly combine to reflect pairwise mating chances (Weis et al., 2014).

In future studies of plant mating patterns, it would be highly informative to have more detailed records on flowering phenology (including assessment of heritability of this trait). This would allow evaluating the extent to which flowering time and pollen dispersal distances predict the observed pairwise mating success. If predictions match the observed patterns, then active forms of kin recognition or other “choices” acting prior or after fertilization are not needed to explain the mating system. If phenology does not succeed in explaining observed matings, one can turn to further processes or traits influencing the mating probability of specific pairs. Examples of such potential factors are size, flower morphology or selective pollen sorting of the pistil (Waser, 1993), but empirical evidence remains scarce. In this respect intraspecific variation of floral traits seem to be the best investigated potential alternative driver of positive assortative mating, because floral display size and nectar reward have been shown to influence pollinator preferences (e.g., Dudash, Hassler, Stevens, & Fenster, 2011; Fenster, Cheely, Dudash, & Reynolds, 2006) and can increase efficiency of pollen transfer (Armbruster, Antonsen, & Pélabon, 2005; Campbell, Waser, & Price, 1996). However, increased floral attractiveness does not necessarily result in enhanced outcrossing, insofar as it instead promotes spatial fidelity of pollinators, enhancing pollination between flowers on the same plant which is problematic for self‐incompatible species (Brys & Jacquemyn, 2010; Harder & Barrett, 1995; Klinkhamer & de Jong, 1993).

We now turn to the limitations of our approach.

A constraint of our study and one that it shares with the study by Monthé et al. (2017) is that the estimation of pollen dispersal has to be based on realized mating events. In other words, to enter the data set, pollen must not only travel a specific distance and land on the stigma but also pass all the hurdles, from appropriate phenology (Fox, 2003; Weis, 2005) through pollen‐pistil interactions (Herrero & Hormaza, 1996) to post‐zygotic selection (Seavey & Bawa, 1986). When observed pollen dispersal distances are assessed based on mature seeds, they will reflect distances where pollen succeeded in completing all the steps. Our resampling procedure of the mating pairs might then not fully represent all potential mating pairs, which introduces a bias if pollen that have travelled specific (particularly short or particularly long) distances are at a relative disadvantage. As the very question we are addressing relates to processes that could cause such disadvantages, it is therefore important to evaluate the potential direction of any error caused by this bias.

In plant populations, biparental inbreeding is more pronounced among short distance pollination events because of common SGS (Vekemans & Hardy, 2004). At the same time, post‐zygotic selection is generally thought to select against inbred offspring (Ghazoul & Satake, 2009; Seavey & Bawa, 1986). If spatially restricted pollen dispersal combines with inbreeding avoidance, this would reduce the frequency of short distance matings and increase realized pollen dispersal distances relative to primary pollen dispersal. If we assume that post‐zygotic selection acts against related pollen, then our ability to detect positive assortativeness would be lowered by this bias. At the same time, it has to be acknowledged that inbreeding avoidance cannot be presumed; some of our results appear to suggest the opposite, as our analysis cannot exclude that trees favour relatives as pollen donors. Selection for avoidance might not occur if inbreeding depression is not strong enough to override the significant inclusive fitness benefit of inbreeding (Kokko & Ots, 2006). This benefit can also be phrased as increased transmission of genes identical by descent (Fisher, 1941) and is well recognized in studies of evolutionary transitions from self‐incompatibility to self‐compatibility (e.g., Charlesworth, 2006; Cheptou & Mathias, 2001; Holsinger, 2000). All this means that it is not easy, in the absence of information on loci impacting kin recognition and/or phenology (or other traits involved in assortativeness), to estimate which part of the mating bias resulted from inbreeding avoidance or preference, and which part reflects assortativeness based on “other' traits”. Should kin recognition be based on specific alleles at a particular locus (as opposed to organisms being able to assess genome‐wide relatedness), the genetic consequences of these two processes become essentially indistinguishable—to the extent that “kin recognition” could now equally well be interpreted as a special case of assortative mating.

Another limitation of our study is that the estimated pollen dispersal distances do not fully represent the entire pollen dispersal kernel. This bias occurs because the seed not assigned to any candidate father are likely to reflect pollen inflow from outside the study area, where the candidate fathers have not been sampled, and these unsampled fathers predominantly reflect distant locations. For all but three species, these fathers can be expected to be less related than the sampled fathers as reassured by the Sp‐statistic. This should, however, not introduce a serious bias to our analysis due to two reasons: First, being able to sample these fathers would reduce the expected mean F under random mating relative to the mean F of the observed matings. Second, systematic reversals of our findings (from positive to negative assortative mating or vice versa) would require that patterns that occur from fathers beyond the boundaries of the study areas deviate strongly from how sampled (unrelated or related) fathers' pollen behaves on stigmas. A father from outside a study area remains unassigned regardless of his relatedness to trees that grow within an area. There may be some bias if only those fathers succeed from far away that they are, for their distance, unusually capable to be accepted by the stigma, but the net effect of such a bias is difficult to ascertain. As study areas and number of observed mating pairs appeared large enough to contain potential mating pairs of varying relatedness, we consider it unlikely that they were too small to yield qualitatively robust results. Nevertheless, it has to be kept in mind that we cannot capture potential siring advantages of genetically very similar or dissimilar pollen sources outside of the study area. Consequently, the findings for the seven species with assignment rates below 50% (Table S1) are only valid at the investigated scales.

In conclusion, pollen dispersal patterns and spatial genetic structuring of adult plants are influenced by the spatial arrangement of adults (Vekemans & Hardy, 2004). Our results suggest that finding elevated inbreeding in populations of conservation concern should be interpreted against the possibility that inbreeding reflects positive assortative mating with respect to phenology (or other traits). Although this does not prove the absence of negative fitness consequences of inbreeding at the population level, it at least shows that extant populations have been able to persist with their current mating system. If assortativeness and/or kin preferences routinely give rise to inbreeding, the consequent level of inbreeding depression might have been reduced by purging of the mutational load (Keller & Waller, 2002). Complacency may be ill advised; however, new, elevated levels of inbreeding can result if habitat fragmentation alters an established balance. Logging or habitat fragmentation can alter realized pollen dispersal distances and impact levels of inbreeding, which can then be an important concern for managing forest genetic resources (Ismail et al., 2012; Robledo‐Arnuncio, Alia, & Gil, 2004; Vinson et al., 2015). Further, fragmented small populations' problems with pollen limitation and consequent reproduction can be aggravated by variation in flowering time (Ison & Wagenius, 2014).

Our study quantifies a mating pattern where genetic similarity can be above, similar or below of what would be predicted based on pollen dispersal distances. While the timing of flowering is a strong contender for the most likely candidate trait for driving the detected instances of positive assortativeness, future studies could usefully investigate the effect of this trait relative to other potential candidate traits for mating patterns. A detailed look on phenology can help studies go beyond any dichotomous view of selfing and outcrossing, and we thus encourage further studies looking at patterns of positive assortative mating as a potential factor influencing the levels of biparental inbreeding in trees.

AUTHOR CONTRIBUTIONS

S.A.I. conceived the study and performed the analysis. S.A.I., and H.K. wrote the paper.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.