Insight into incipient reproductive isolation in diverging populations of Brachionus plicatilis rotifer

The emergence of reproductive isolation is key in maintaining within- and between-species diversity and one of the initial steps of speciation. In the Iberian Peninsula, the diverging populations of the Brachionus plicatilis rotifer create an ideal system to shed light on the mechanisms that give rise to the emergence of reproductive isolation. Herein, we quantify the degree of behavioural reproductive isolation in two groups of B. plicatilis populations, namely, neighbouring populations diverging by adaptation to the local environment and populations diverging in the absence of gene flow due to geographic distance. We conduct behavioural no-choice assays to test mating reproductive isolation between these populations. The analysis shows signatures of ongoing behavioural reproductive isolation in most of the population crosses, which is more pronounced in populations with a higher level of adaptive divergence, presumably under high migration rates. Overall, this study suggests that local adaptation is associated with mating behaviour resulting in reproductive isolation.


Introduction
Reproductive isolation has strong implications for both genetic differentiation and evolutionary dynamics, affecting population structure, local adaptation, speciation and sister species coexistence (Coyne & Orr, 2004;Weber & Strauss, 2016). There is a wide range of intrinsic and extrinsic factors that could contribute to isolation. The causal loops that might occur among them make the dissection and testing of such a panoply complicated (Palumbi, 1994;Coyne & Orr, 2004).
One of the factors leading to reproductive isolation between populations may result from the difficulty of immigrants arriving, surviving and/or becoming sexually involved in the recipient population (Sexton et al., 2014;Turbek et al., 2018). If arrival rate is the prevalent factor, then a pattern of isolation by distance (IBD) is expected to emerge due to the combined, opposite effect of migration and genetic drift (Wright, 1943;Church & Taylor, 2002). On the other hand, if unmatched timings of dispersal/reproduction and survival are the prevalent factors, then a pattern Abstract The emergence of reproductive isolation is key in maintaining within-and between-species diversity and one of the initial steps of speciation. In the Iberian Peninsula, the diverging populations of the Brachionus plicatilis rotifer create an ideal system to shed light on the mechanisms that give rise to the emergence of reproductive isolation. Herein, we quantify the degree of behavioural reproductive isolation in two groups of B. plicatilis populations, namely, neighbouring populations diverging by adaptation to the local environment and populations diverging in the absence of gene flow due to geographic distance. We conduct behavioural no-choice assays to test mating reproductive isolation between these populations. The analysis shows signatures of ongoing behavioural reproductive isolation in most of the population crosses, which is more pronounced in populations with a higher level of adaptive divergence, presumably under high migration rates. Overall, this study suggests that local adaptation is 1 3 Vol:. (1234567890) of isolation by environment (IBE) is expected (Nosil et al., 2005;Wang & Bradburd, 2014). Local adaptation is likely to be involved in this second case (Nosil et al., 2002). IBD and IBE have been found in a broad range of organisms (IBD, e.g., Zeller et al., 2006;Mills et al., 2007;Campillo et al., 2011;Ventura et al., 2014;IBE, e.g., Shafer & Wolf, 2013;Martin et al., 2021) and may act together (Sexton et al., 2014;Weber et al., 2017). In both scenarios-IBD and IBE-low gene flow would promote behavioural reproductive isolation-i.e., based on mate preferences-due to genetic drift acting on genes that affect behaviour. Moreover, the genetic hitchhiking of behavioural genes linked to genes for local adaptation might play a role in the case of IBE (Hawthorne & Via, 2001;Wu, 2001, but see Parchman et al., 2013). In addition, if genetic divergence, either adaptive or nonadaptive, is associated with at least some degree of post-mating isolation, then selection for behavioural isolation would occur because it results in an optimized allocation of parental resources, i.e., avoiding the production of hybrid offspring. Thus, the degree of behavioural isolation is expected to increase with the geographic and environmental distance between populations and concomitantly with the genetic distance in neutral markers. Notably, regardless of how behavioural isolation arises, it is expected to reinforce genetic divergence in neutral markers, which illustrates the two-way causal links that might be at work.
Brachionus plicatilis Müller 1786 is a cyclical parthenogenetic rotifer that inhabits saline and brackish ponds. Its predominant type of reproduction is ameiotic parthenogenesis (Gilbert, 1963;Wallace et al., 2015;Serra et al., 2018). Its natural populations, which dwell in the water column temporarily within a year, are composed mainly of asexual (also called amictic) females. Sexual reproduction is initiated at high population densities (Stelzer & Snell, 2003) and starts with the production of sexual (mictic) daughters (Snell & Boyer, 1988). Sexual females, which produce meiotic eggs, have two mutually exclusive reproductive fates: to become male producers if not inseminated by a male while they are young or to become diapausing egg producers otherwise . Sexual and asexual females are morphologically identical, but they can be identified according to the eggs they carry.
Mating behaviour in rotifers is a three-step sequence: male-female encounter, circling and copulation. As rotifers do not secrete any sexual pheromones into the environment (Snell et al., 1995), male-female encounters occur randomly (Snell & Garman, 1986). After the encounter, the ciliated corona of the male contacts a female's body (Snell & Morris, 1993;Snell et al., 1995), and the male displays a circling behaviour, i.e., swimming around the female while maintaining contact by his corona. After several circlings, the male locates either the female's corona or the foot opening, and hypodermic insemination takes place (Gomez & Serra, 1995). The whole mating behaviour can take from tens of seconds to several minutes, and the threestep sequence can be interrupted at any time (Snell & Hawkinson, 1983;Snell & Hoff, 1987). A surface glycoprotein (mate recognition protein, MRP) present on the female's body is involved in recognition of a mate (Snell & Morris, 1993;Snell et al., 1995;Gribble et al., 2011). Although mate choice is mostly driven by males, females may exhibit a reaction by varying their resistance to male circling (Snell et al., 2007). Female susceptibility to fertilization and male capacity for fertilization decline with age (Snell & Childress, 1987). Male rotifers tend to copulate more often with young females but show no preference for sexual over asexual females (Gomez & Serra, 1996).
In the Iberian Peninsula, populations of the cyclically parthenogenetic rotifer B. plicatilis create an ideal system for studying intrinsic factors resulting in some level of reproductive isolation. The species consists of highly divergent lineages , reflecting strong founder effects with subsequent resource monopolization by the founder genotypes. Additionally, Iberian populations of B. plicatilis inhabit highly diversified habitats, from ephemeral-and highly environmentally unpredictable-to permanent ponds. They show local adaptation affecting their sexual reproduction patterns, adjusting them to variation in environmental predictability (i.e., high propensity to sex in populations with low environmental predictability Franch-Gras et al., 2017a, b). A recent study on the diversity of the gene mmr-b, which is involved in mate recognition, has shown that despite being under stabilizing selection, a correlation is found between gene diversification and differences in environmental factors when comparing B. plicatilis populations (Jezkova et al., 2022). Brachionus plicatilis belongs to a species complex, and the divergence of mmr-b among the species of the complex has also been documented (Gribble & Welch, 2012).
In this research, we tested the hypothesis that evolution of premating behavioural isolation occurs between populations of the species B. plicatilis. We hypothesized that the ecological divergence is correlated with that isolation. We used bidirectional nochoice mating tests to investigate within-species mating preferences in populations belonging to different phylogeographic groups or showing ecological divergence. The results were used to assess the consistency of different evolutionary hypotheses, and phylogeographical and ecological divergences were found to play a role in the establishment of behavioural reproductive isolation.

Study populations
We selected two groups of populations of B. plicatilis inhabiting shallow ponds in eastern Spain. One group was based on known phylogeographic structure, with each population belonging to a different clade, and the other was based on known differentiation in adaptation to local conditions (Table 1). For the group of phylogeographically structured populations (PS), we chose three populations, each with a dominant haplotype belonging to divergent mtDNA clades and showing a high degree of pairwise population differentiation   (Table 2). The selected  Table 2.

Clone isolation and culture conditions
The sediment of the ponds was sampled in 2013. Diapausing eggs were extracted from the sediment in 2017 using a sugar flotation technique (i.e., , and rotifer clones were established from isolated, individually hatched eggs. As B. plicatilis belongs to a cryptic species complex, hatchlings were taxonomically identified based on restriction fragment length polymorphism (RFLP) analysis of a fragment of the mitochondrial gene COI (Campillo et al., 2005). Individual clones (30-40 per population) were maintained in the laboratory in 15 ml stock cultures and grown in a culture of the microalgae Tetraselmis suecica in artificial seawater (Instant Ocean® Sea Salt, Aquarium Systems) at 12 g l −1 salinity, f/2 enriched medium (Guillard & Ryther, 1962), and at 25°C under constant illumination (PAR: approx. 35 µEm −2 s −1 ). These conditions, hereafter referred to as the standard, were also used with minor modifications in the experiments. Every week, a small volume of each culture was transferred to fresh medium, and the rest was discarded.

Mating experiments
Pre-experimental cultures were established by placing five asexual, ovigerous females of each clone into a flask with 40 ml artificial seawater and T. suecica algae at a concentration of 2 × 10 5 cells ml −1 . The aim was to achieve exponential growth, to induce sexual reproduction and to obtain a high abundance of both males and females. Thus, the flasks were checked daily until sexual reproduction was observed (i.e., when the first males were detected, after 4-7 days).
In the mating assays, young individuals (not older than 3 h) of both sexes were used to maximize the potential for mating responses. To collect virgin newborns, eggs were detached from the female's body by vigorous shaking of the flasks. This method has been proven efficient and harmless for eggs and hatchlings (Gomez & Serra, 1996). The eggs were sorted according to their type (smaller-male-and largerfemale-eggs), collected using a micropipette, placed separately in the standard conditions with no algae addition and allowed to hatch. For each mating assay, one clone provided females (25 individuals), and a  Gomez et al. (2002; for PS). All pairwise F ST values were significant at P < 0.001. The lower 3 × 3 hemi-matrices show the geographic distance in km (see Table 1 Gomez & Serra (1995), but the observation time was increased to 10 min. All the females were placed together in 50 µl fresh medium at 25°C with no algae, and males were added sequentially. After the addition of each male, the mating behaviour was observed for 10 min and then replaced with another male. This was repeated until all six males were used. The experimental volume was kept constant, as individual males were added and removed in a small volume of medium of approximately 2 µl. Observations were conducted under a stereomicroscope (Olympus SZX10, Japan), and the number of encounters, circlings and copulations was recorded according to Gomez & Serra (1995) and Snell & Hawkinson (1983). We defined an encounter as any physical contact between the male's corona and the female's body, a circling (implying previous encounter) was recorded when the male swam around the female body at least once, and copulation occurred when the male penetrated the female's soft body parts (corona or foot opening). Circling always preceded copulations. Nine population crosses were carried out for each population group (PS and LA) using males and females in all possible pairwise combinations (i.e., three intrademic and six interdemic combinations). Each population cross was replicated five times (9 crosses × 2 population groups × 5 replicates = 90 mating assays; note that six sequential measures-one per male-were obtained for each mating assay). Replication consisted of changing the clones in the mating assays. No clone was used in more than one assay within a group (5 populations × 30 clones per population = 150 clones). The experiment lasted 4 months, and one to two, rarely three, mating assays were performed per day. A single observer (I.J.) recorded the data between March and July 2018, and a previous anonymization of the experimental cultures was adopted to prevent observation bias.

Data analysis
Variation in individual male activity within an assay, albeit included in the experimental procedure, was not considered a factor in our data analysis (i.e., the counts from six males were summed up). For each population combination, we calculated the average percentage of circlings resulting from encounters (hereafter, "circling after encounter"), the average percentage of copulations resulting from circlings ("copulation after circling") and the average percentage of copulations resulting from encounters ("copulation after encounter"), the last one being a compound of the other two.
We used a GLMM (as implemented in the lme4 package in R software, version 3.6.2) (Bates et al., 2015;R Core Team, 2019) to analyse the significance of the effects in the experimental design on copulations after encounters. Our restriction to this response was oriented to avoid dependence in our statistical analysis while favouring the biological relevance in focus. GLMM was applied separately to the PS and LA study groups. Based on the Bayesian information criterion, we assumed a binomial error distribution (assessed against a Poisson distribution with fate-copulation or not-as an additional factor). For the selected distribution, its canonical link function (logit) was used. The factors in the GLMM were (1) male-providing population (MP), (2) female-providing population (FP), (3) their interaction (MP:FP) and (4) clone pair as a random factor. Notice that the MP:FP interaction accounts for a deviation from independence (i.e., when a specific population combination deviates from the average copulations after encounters of the two populations involved). Therefore, to find a behavioural propensity against interdemic mates, a significant MP:FP interaction is necessary; however, the MP:FP effect of intrademic crosses must also be higher (i.e., relatively more copulations vs. encounters) than in interdemic crosses. Hence, to quantify this propensity between population crosses, we estimated coefficients associated with the interaction effect. These coefficients were obtained by additively decomposing the proportion of copulations (vs. encounters) observed in a mating assay. The components were (1) a (grand) mean value, (2) a component due to the male-providing population (mean deviation of that population when providing males), (3) a component due to the female-providing population (mean deviation of that population when providing females) and (4) the residual, which retains the interaction (MP:FP) and the random variation among mating assays. The decomposition follows the approach in Tortajada et al. (2009). The 45 coefficients accounting for the MP:FP interaction for each mating assay in each study group (PS or LA) were inspected to check whether the MP:FP of intrademic crosses was higher than that of interdemic crosses.
Additionally, to compare with other studies, we also computed the relative deficit-excess of interdemic copulations using a reproductive isolation index (RI, Sobel & Chen, 2014), although it presents caveats, such as that this index is affected by the indiscriminating mating activity, and mating activity might be largely sensitive to the experimental environment. RI, computed as follows: 1 − 2 interdemic copulations total copulations (range from − 1 to + 1) indicates the strength of the effective prezygotic barrier (opposing gene flow) in the experimental conditions. One-sided t test on Spearman's rank correlation coefficient using the cor.test function implemented in R software, version 4.1.2. (R Core Team, 2021) was applied to assess the correlation between our mating behaviour metrics (vectors of interaction coefficients and RI) and different variables linked to the population characteristics (genetic differentiation in neutral traits, geographic distance and differentiation in environmental predictability, the latter for LA only).

Results
From the observation of the 540 males (90 mating assays; intrademic and interdemic crosses), we recorded 4,318 encounters (LA 2,055; PS 2,263), 3,191 circlings (LA 1,610; PS 1,581) and 1,850 copulations (LA 853; PS 997). Raw data from the mating experiment can be found in Supplementary Table 1. The percentages of copulations after encounters are shown in Table 3. They averaged 48% for intrademic crosses and 42% for interdemic crosses. When comparing pairwise intrademic and interdemic crosses that share a population, a higher percentage of copulations was found in intrademic crosses in 15 out of 24 comparisons. However, only SAL (in the two groups) and HYB (from the LA group) had a higher percentage of copulations for intrademic crosses in all comparisons (LA group: SAL intrademic 69.5% vs. interdemic 19.4-45.0%; HYB intrademic 56.2% vs. interdemic 19.4-50.8%, PS group: SAL intrademic 64.4% vs. interdemic 34.7-59.0%; for disaggregated values, see Table 3). In general, circling after encounters and copulation after circling had similar qualitative patterns as copulation after encounters ( Supplementary Fig. 1). An exception was the CHI population (PS group), where the percentages of copulations after encounters were low and associated with high rates of encounter in intrademic crosses, while subsequent circlings and copulations did not show values remarkably different from those in other populations.
The results of the GLMM for copulation after encounter showed significant effects for all fixed factors [i.e., male-providing population (MP), female-providing population (FP), as well as its interaction (MP:FP) for each study group (Table 4)]. Figure 1 ranks the coefficients for an interaction effect obtained from additively decomposing the percentage of copulations after encounters. For both study groups, the average coefficient is higher for intrademic crosses than for interdemic crosses, indicating a higher propensity towards intrademic over interdemic crosses. This trend was stronger in the LA study group than in the PS group. The positive interaction associated with intrademic crosses is largely due to SAL and HYB populations. At the same time, in the PS study group, crosses between males from CHI and females from TON showed rather high positive values (Table 5). The prezygotic isolation index, RI, was in general slightly positive for all population combinations, pointing to weak behavioural reproductive isolation (Table 6). Spearman's correlation between the interaction coefficients of the two study groups and population characteristics (predictors: difference in environmental predictability, geographic distance and F ST values in neutral traits) were, except for geographic distance in the PS group, significantly negative, as expected after a preference for intrademic crosses ( Table 7). The correlation coefficients between RI and the predictors were positive (as expected for a preference for intrademic crosses), except for F ST of microsatellites in the PS group. However, these correlations were not statistically significant.

Discussion
Our results show that populations of B. plicatilis in the Iberian Peninsula evolved a low but not negligible level of behavioural mating isolation. These populations have a deep genetic divergence associated with post-glacial recolonization Campillo et al., 2011) and/or ecological divergence (Franch-Gras et al., 2017b). We documented the tendency for mating preferences using two metrics: interaction coefficients for copulations after encounters and the reproductive isolation index, RI. Both metrics are qualitatively consistent in suggesting a preference for intrademic mating. Positive interaction coefficients were mostly found for intrademic crosses, and RI had mostly positive values (five out of six comparisons), indicating a tendency for reproductive isolation. RI values are not expected to be high for populations of the same species inhabiting the same region. The RI between closely related Brachionus species ranges from 0.4 to 1.0 (Gribble & Welch, 2012), indicating incomplete behavioural isolation even between species. Overall, our results are compatible with the presence of partial behavioural isolation among our studied populations.
Despite showing the same trend, we regard our results on the interaction coefficients for successful encounters as more reliable than comparisons based on RI. Estimates of RI might be affected by differences in mating activity of the genotypes. As male-female encounters are random (Gilbert, 1963;Snell & Hawkinson, 1983), this activity depends on swimming speed, among other factors. Thus, an effect of the laboratory environment on RI is possible if, for instance, the environment favours the performance of individuals of some populations over others. Thus, following the approach in Tortajada et al. (2009), we used the interaction coefficients for encounters resulting in copulation because they are additive departures from the expectancy based on the general mating activity of the concerned genotypes. Additionally, we focused on copulations resulting from encounters because copulations have stronger implications for isolation than the other events.
Several mechanisms for behavioural mating isolation between genetically divergent populations have been recognized in the literature, namely, genetic drift acting on mating behaviour-coding genes (Kaneshiro, 1980;Carson & Templeton, 1984), indirect selection if those genes are physically linked to loci under disruptive selection and associated with local adaptation (Dobzhansky, 1937in Gavrilets, 1999Coyne & Orr, 2004) and direct selection if outbreeding fitness depression exists (Coyne & Orr, 2004). Conversely, when a preference for intrademic mating occurs, it should work as an intrinsic gene flow barrier enhancing genetic divergence. These mechanisms predict a negative correlation between genetic differentiation and behavioural mating. Our results show signatures of such a pattern in the two groups of populations  Table 7 Spearman's rank correlation (ρ) between (1) reproductive isolation metrics (the interaction coefficient after an additive decomposition of the proportion of copulation after encounters and the reproductive isolation index RI) and (2)  studied here. The pattern is supported by our Spearman correlation analysis on the interaction coefficients, suggesting that, whichever the mechanisms at work are, their effects are stronger than the stabilizing selection expected to act on the mate recognition systems (Brooks et al., 2005;Smadja & Butlin, 2009), a selection whose benefit is to not become unrecognizable for a potential partner. Our results contrast with those of Berrieman et al. (2005), who, using a single clone for each of four populations, did not find evidence for reproductive isolation between the northern and southern Iberian clades of B. plicatilis. By comparing the two experimental groups of populations, our results provide insights into the association of behavioural mating isolation with IBE and IBD. First, signatures for intrademic mating preferences are stronger in the LA group than in the PS group, despite the former populations belonging to the same phylogeographic clade  and lying in geographical proximity. This suggests a stronger association of mating behaviour with IBE than with IBD and contrasts with the strong evidence for IBD in Iberian B. plicatilis populations when neutral markers were studied . Second, we did not find a significant correlation between geographic distance (IBD) and mating isolation in the PS group, while in the LA group, mating isolation was slightly more correlated with environmental divergence than with geographic distance. With both predictors being collinear, the correlation with geographic distance in LA is likely spurious to some extent. As another piece of evidence, genetic differentiation in the mmr-b gene, which codes for the receptor responsible for mating recognition, increases with ecological distance (Jezkova et al., 2022).
There are several findings worthy of note.
(1) The populations with the lowest differentiation in neutral traits (RAS-HYB; Montero-Pau et al., 2017) showed no mating discrimination when considering the RI index or even preference for interdemic crosses in the case of interaction coefficients. The localities of these two populations are situated very close to each other (< 1 km) and are likely to frequently merge. This suggests the role of migration in blurring differentiation.
(2) Due to the historical colonization pattern, the coastal lineage and the northern lineage are the closest lineages Serra et al., 2019). This could explain why mating between TON and CHI (from the coastal and northern clades, respectively) shows lower avoidance of interdemic crosses than shown by mating between either population with SAL (the southern clade). (3) In both population groups, the SAL population drives the most distinguishable pattern of behavioural isolation, so the caveat is whether our argument is based only on one single population. We notice, however that behavioural isolation in the SAL-HYB pair is higher than that in the SAL-RAS pair, which is the expected ranking in an IBE scheme.
The association between IBE and mating isolation found herein makes the question of selection of isolation (Butlin, 1987) worthy of investigation in these rotifers. An important condition for selection for premating isolation is that the intensity of the potential gene flow must be intermediate, meeting a balance between the homogenizing effect of extensive gene flow and lack of gene flow (no need for selection for isolation; (Nosil, 2013;Yukilevich, 2012). Even if direct measures of migration rates for rotifers are scarce (Lopes et al., 2016;Moreno et al., 2019), a gradient of rates is most likely in populations such as those in the LA group.
In our research, we focused on behavioural mating isolation as a premating (i.e., low-cost) intrinsic barrier. As anticipated above, other intrinsic barriers not captured here might be involved in premating isolation in the wild. In Brachionus rotifers, the beginning of sexual reproduction is triggered by the accumulation of a crowding chemical signal in the environment Snell, 2017). While it was shown that the signal is rather conservative and sex can be cross-induced between species García-Roger et al., 2009), the threshold concentration varies even among genotypes of a single species (Franch-Gras et al., 2017a). These differences could affect the timing of sex, thereby promoting reproductive isolation between immigrants and residents via their allochronic sexual periods.
How populations differentiate, become locally adapted and eventually diverge into distinct species are persistent questions in evolutionary biology, and their answer depends upon a panoply of processes acting in opposing directions, with complex interactions due to causal feedbacks and having different strengths in relation to organisms and habitat features. The taxon of cyclically parthenogenetic rotifers harbours high species richness, and the few well-studied species in the taxon are clusters of genetically differentiated and locally adapted populations. Our results suggest that isolation may arise associated with genetic divergence and might play a role in the diversification of rotifers. We nevertheless remark that behavioural isolation is just a piece of the premating isolation puzzle.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. The research was supported by a Grant from the Valencian Conselleria de Educación, Investigación, Cultura y Deporte (AICO/2020/013), a Grant CGL2015-65422-P funded by MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe. IJ was supported by a Grant from the Valencian Consellería de Educación, Investigación, Cultura y Deporte (GRISOLIAP/2017/096).
Data availability Enquiries about data availability should be directed to the authors.

Conflict of interest
The authors declare no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.