Mitonuclear discordance and patterns of reproductive isolation in a complex of simultaneously hermaphroditic species, the Allolobophora chlorotica case study.

Historical events of population fragmentation, expansion and admixture over geological time may result in complex patterns of reproductive isolation and may explain why, for some taxa, the study of mitochondrial (mt) and nuclear (nu) genetic data results in discordant evolutionary patterns. Complex patterns of taxonomic diversity were recently revealed in earthworms for which distribution is largely the result of paleogeographical events. Here, we investigated reproductive isolation patterns in a complex of cryptic species of earthworms in which discordant patterns between mt and nu genetic lineages were previously revealed, the Allolobophora chlorotica aggregate. Using four nu microsatellite markers and a fragment of the cytochrome c oxidase subunit I mt gene, we carried out a parentage analysis to investigate the mating patterns (i) between individuals belonging to two divergent mt lineages that cannot be distinguished with nu markers and (ii) between individuals belonging to lineages that are differentiated both at the mt and nu levels. Amongst the 157 field-collected individuals, 66 adults were used in cross-breeding experiments to form 22 trios based on their assignment to a mt lineage, and 453 obtained juveniles were genotyped. We showed that adults that mated with both their potential mates in the trio produced significantly more juveniles. In crosses between lineages that diverged exclusively at the mt level, a sex-specific pattern of reproduction characteristic to each lineage was observed, suggesting a possible conflict of interest concerning the use of male/female function between mating partners. In crosses between lineages that diverged both at the mt and nu levels, a high production of cocoons was counterbalanced by a low hatching rate, suggesting a post-zygotic reproductive isolation. Different degrees of reproductive isolation, from differential sex allocation to post-zygotic isolation, were thus revealed. Lineages appear to be at different stages in the speciation process, which likely explain the observed opposite patterns of mitonuclear congruence.

INTRODUCTION

Over the past decades, an increasing number of studies investigating phylogenetic and phylogeographic hypotheses have reported discordant patterns between nuclear (nu) and mitochondrial (mt) DNA markers (Toews & Brelsford, 2012). Causes of mitonuclear discordance are numerous and may result from incomplete lineage sorting, sex‐biased dispersal, asymmetrical introgression, natural selection or Wolbachia‐mediated genetic sweeps (Toews & Brelsford, 2012). But most often, the distinct mode of transmission between mtDNA and nuDNA, i.e. maternal versus biparental, is sufficient to explain that they respond differently to demographic fluctuations (Després, 2019). Indeed, the mtDNA is a haploid cytoplasmic non‐recombinant genome maternally transmitted and, therefore, has a four‐fold lower effective population size than the nu genome. Thus, when a population is fragmented into small isolates, such as during glacial periods, genetic drift will act more strongly on mtDNA, which will diverge more rapidly than the nu genome, which will retain a greater variability. Such periods of allopatry between isolates may lead to the formation of distinct genetic lineages that, upon subsequent geographic and/or population expansion, will eventually come into secondary contact. The genetic outcome of this secondary contact will depend on the mechanisms of reproductive isolation that have occurred during the allopatric phase. In the absence of reproductive isolation, nu genomes will fully recombine, whereas divergent mtDNA haplotypes will be retained. In the presence of reproductive isolation mechanisms, either speciation processes during the allopatric phase are already at an advanced stage and prevent outcrossing between lineages, or they are incomplete, and then heterogeneous gene flow can be observed between lineages due to selections against hybrids (Després, 2019).

Earthworms are simultaneous hermaphrodites that may reproduce sexually or by parthenogenesis. Their distributions have been largely shaped by paleogeographical and paleoecological constraints as a result of their hypothesized ancient age and low dispersal capabilities (Bouché, 1972; Novo et al., 2011). Over the past decades, studies that have used molecular techniques to explore earthworm systematics, taxonomy and phylogeography found an unprecedented cryptic diversity (e.g. Dupont et al., 2011; King et al., 2008; Novo et al., 2010; Taheri et al., 2018) and, in particular, in a few sexually reproducing species, cases of mitonuclear discordance were described (e.g. Dupont et al., 2016; Torres‐Leguizamon et al., 2012). Using such closely related earthworm species to investigate mechanisms of reproductive isolation (e.g. Jones et al., 2016; Marchán et al., 2018) could, therefore, provide new insights to explain these discrepancies. In that context, the Allolobophora chlorotica complex of earthworm species is a good model. It is a sexually reproducing species that was previously described as an aggregate of several closely related lineages (Dupont et al., 2016) but shows discordant patterns between mt and nu DNA markers. Allolobophora chlorotica exists as two colour morphs, green and pink (King et al., 2008; Satchell, 1967). The green morph represents a single species composed of two divergent mt lineages that cannot be distinguished with nu markers (i.e. L2 and L3). The taxonomic status of the pink morph is less clear, even though it is believed to be composed of at least three species corresponding to three mt lineages (i.e. L1, L4 and L5) also distinct at the nu level (Dupont et al., 2011, 2016). Field and laboratory observations showed that the green morph tends to be more common in wet soils and the pink morph in dry soils (Lowe & Butt, 2007; Satchell, 1967). Thus, although it is often stated that these lineages have a sympatric distribution (i.e. they co‐occur in a region), some of them at least do not live in syntopy (i.e. they do not use the same habitat, Rivas, 1964). These distinct ecological preferences led Lowe and Butt (2007) to suggest that soil moisture acts as a prezygotic barrier to intermorphic mating in syntopic populations of these simultaneous hermaphrodite earthworms that do not self‐fertilize (Dupont et al., 2011). Breeding experiments additionally reported mechanisms of postzygotic isolation between colour morphs. The viability of cocoons resulting from the crossing between the two colour morphs is severely restricted and the male offspring from the backcrossing of hybrids with pure bred morphs are sterile (Lowe & Butt, 2008). The genotyping of individuals collected in two natural populations also supported the hypothesis of a process of reproductive isolation between the L1 lineage of the pink morph (referred to as ‘pink’ in the rest of the paragraph) and the L2/L3 lineages of the green morph (referred to as ‘green’ in the rest of the paragraph; Dupont et al., 2016). Dupont et al. (2016) showed a high level of congruence between the assignment based on the mt COI and the nu microsatellites, suggesting that hybridization is uncommon between pink and green morphs. However, they also reported cases of introgression, with a few individuals having a pink mtDNA haplotype but assigned to a nu cluster grouping all the green individuals (Dupont et al., 2016). These individuals could result from multiple generations of unidirectional hybridization between pink females and green males. Each backcross with a green individual would dilute the proportion of pink nu allele by half until the population had overwhelmingly accumulated nu green alleles but retained the maternally inherited pink mtDNA haplotype. The authors also recorded one individual having a green mtDNA haplotype but assigned to the nu cluster grouping the majority of the pink individuals, suggesting that unidirectional hybridization could also happen from mating between green females and pink males (Dupont et al., 2016). These events of unidirectional hybridization reinforce the hypothesis of sterility of the male function in pink–green crosses.

Here, we combined cross‐breeding experiments and a parentage analysis based on the sequencing of the cytochrome c oxidase subunit I mt gene and nu microsatellite markers, to explore the reproductive isolation mechanism within this intriguing complex of closely related earthworm lineages, the A. chlorotica aggregate. First, we investigated whether earthworms of this species complex may have offspring from multiple mates during the same mating period. Second, we examined mating patterns (i) between divergent mt lineages which could not be differentiated by nu markers in previous studies (L2 and L3; Dupont et al., 2011, 2016) and (ii) between lineages that diverged both at the mt and nu levels (L1 and L2/L3), in order to detect evidence of reproductive isolation processes in progress. Third, we examined whether one function (male or female) is preferentially used in mating between these lineages.

MATERIAL AND METHODS

The work described below corresponds to (i) the analysis of cocoon production and hatching rate and (ii) the parentage analysis of the hatched offspring, of A. chlorotica individuals collected in the field. The setup of the cross‐breeding experiments described below is based solely on the mtDNA assignment of the individuals collected in the field. To assign the parentage of their offspring, we used mt markers and we additionally confronted these results with the assignment based on the use of nu markers (i.e. microsatellites). This genotyping, using microsatellite markers, also made it possible to characterize the genetic structure of the field population at a fine scale.

Field sampling and genetic analysis of collected specimens

The sampling was carried out at a farm site (Walton Hall Farm, Preston, UK) in a field grazed by cattle. The soil is an alluvial, sandy clay with a pH of 8.3. At this site, A. chlorotica was found in high‐density clusters (i.e. patches). Regardless of their colour morph, we collected a total of 157 A. chlorotica individuals: 46 juveniles in 2012 and 71 juveniles and 40 adults (i.e. clitellate) in 2015. From each individual, we removed the last segments of the caudal section and preserved it in 96% ethanol before DNA extraction using the NucleoSpin® 96 Tissue kit (Macherey‐Nagel). It is worth noting that the ablation of the last segments of the caudal section of earthworms does not affect their survival due to their regeneration capacity (e.g. Xiao et al., 2011).

Assignment of the field‐collected individuals to Allolobophora chlorotica mitochondrial lineages

To assign the collected individuals to A. chlorotica mitochondrial lineages, we amplified and sequenced the fragment of the cytochrome c oxidase subunit I mt gene (COI), proposed as a standard DNA barcode for animals (Hebert et al., 2003). Individuals collected in the field population in 2012 were processed at the Canadian Centre for DNA barcoding (CCDB), in the context of the global earthworm barcoding campaign (EarthwormBOL, Rougerie et al., 2009) and as part of the International Barcode of Life Initiative (iBOL). PCR amplifications and DNA sequencing were performed according to the standard protocols used in CCDB (Hajibabaei et al., 2005). For the individuals collected in the field population in 2015, the COI gene was amplified using the primer pair described in Folmer et al. (1994). DNA sequencing was performed by Eurofins Genomics company and we manually aligned the sequences using the BioEdit program (Hall, 1999). We inferred the mt lineage of the individuals using the identification engine of BOLD (Barcode of Life Data Systems – https://www.boldsystems.org/).

Genetic structure of the field population

To explore the mt genetic variation in the field population, we estimated haplotype frequencies using the software DNASP 6.0 (Rozas et al., 2017). We then examined the relationships amongst haplotypes using a haplotypic network constructed using the software NETWORK 10 (www.fluxus‐engineering.com/sharenet.htm). The network was obtained by applying the reduced median algorithm (Bandelt et al., 1995) and was postprocessed using maximum parsimony calculations to reduce the number of superfluous network links.

To explore the nu genetic variation in the field population, we genotyped individuals at four highly polymorphic microsatellite loci, Ac127, Ac170, Ac418 and Ac476, as described in Dupont et al. (2011). We amplified the loci by polymerase chain reaction (PCR) in one multiplex set and in 12 μl reactions using 10 ng of DNA and the Qiagen ® Multiplex Kit, according to the manufacturer's protocol. The migration of the PCR products was carried out on an ABI 3130 xl Genetic Analyzer using the LIZ500 size standard (Applied Biosystems); alleles were scored using GeneMapper 5 software (Applied Biosystems). We then calculated all basic genetic parameters including allele frequencies, number of alleles (N all), polymorphic information content (PIC) for each locus and the observed (H o) and expected (H e) heterozygosity using CERVUS v.3.0.7 (Kalinowski et al., 2007; Marshall et al., 1998). We tested for the null independence between loci from statistical genotypic disequilibrium analysis using Genepop V4.4 (Rousset, 2008). The significance of a deviation from the Hardy–Weinberg equilibrium including a Bonferroni correction and null allele frequencies were estimated using CERVUS v.3.0.7.

The admixture model of the STRUCTURE software (Pritchard et al., 2000) was used to identify potential hybrids in the field population by modelling cluster assignments for K = 1–5 clusters. We made 10 independent runs for each K to confirm consistency across runs. In all simulations, we performed a burn‐in period of 10 000 iterations and 1 000 000 Markov chain Monte Carlo iterations. To determine the most likely value of K, we used the ΔK method of Evanno et al. (2005) implemented in STRUCTURE HARVESTER (Earl & Vonholdt, 2012). We combined the results from the 10 replicate runs into one output using CLUMPP software (Jakobsson & Rosenberg, 2007).

Laboratory cross‐breeding experiment

We kept the 71 juveniles collected in 2015 in individual pots and placed them in an incubator until they reached adulthood, in order to have virgin adults. We performed cross‐breeding experiments in trios to test for multiple paternity and reproductive success between individuals from divergent lineages. On the basis of the available individuals from the different lineages, we used a subset of 66 adults to form 22 trios that we labelled from letter A to V (Table 1). We formed six types of trios with up to two adults per lineage. In detail (and see Table 1), these trios were composed of: (A–E) two adults of L2 and one adult of L3, (F–J) two adults of L3 and one adult of L2, (K–L) three adults of L2, (M–N) three adults of L3, (O–R) two adults of L2 and one adult of L1, (S–V) two adults of L3 and one adult of L1. Trios composed of adults from two lineages (A–E, F–J, O–R, S–V) were replicated 4–5 times and trios composed of adults of the same lineage (K–L, M–N) were replicated twice (see Table 1). We collected the cocoons produced by each trio monthly for four months and kept them in an incubator until hatching. Upon hatching, the juveniles were fixed in ethanol before we performed DNA extraction, as described above.

TABLE 1

For each trio, the assemblage of adults according to their mitochondrial lineage, the number of cocoons and juveniles produced and their hatching rate, plus the number and proportion of juveniles genotyped

Cross‐breeding type		Offspring			Genotyping		Summary per cross‐breeding type
Trio	Assemblage	# Cocoons	# Juveniles	% Hatching	# Genotyped juveniles	% Genotyped juveniles	# Cocoons (mean ± SD)	% Hatching (mean ± SD)
A	L2/L2/L3	58	51	87.9	22	43.1	52.2 ± 9.3	75.7 ± 10.1
B		63	49	77.8	21	42.9
C		39	32	82.1	22	68.8
D		53	35	66.0	21	60.0
E		48	31	64.6	22	71.0
F	L2/L3/L3	9	6	66.7	6	100.0	33.6 ± 17.4	72.8 ± 14.3
G		45	24	53.3	21	87.5
H		26	23	88.5	20	87.0
I		34	24	70.6	18	75.0
J		54	46	85.2	21	45.7
K	L2/L2/L2	28	21	75.0	21	100	40.0 ± 17.0	78.8 ± 5.4
L		52	43	82.7	21	48.8
M	L3/L3/L3	25	15	60.0	14	93.3	27.0 ± 2.8	71.4 ± 16.1
N		29	24	82.8	20	83.3
O	L1/L2/L2	30	20	66.7	20	100.0	50.2 ± 15.8	60.0 ± 6.6
P		68	44	64.7	24	54.5
Q		55	30	54.5	24	80.0
R		48	26	54.2	25	96.2
S	L1/L3/L3	70	38	54.3	22	57.9	58.2 ± 10.8	58.2 ± 10.2
T		64	47	73.4	23	48.9
U		46	24	52.2	21	87.5
V		53	28	52.8	23	82.1
		997	681	68.9 ± 12.5	452	73.3 ± 20.1

Parentage assignment of the offspring from the cross‐breeding experiment using COI mitochondrial marker

We determined the mt lineage of each hatched juvenile that is the lineage of the parent that used its female function during mating, using three methods, each specific to each cross‐breeding type. (1) For the trios of purebred of L2 or L3 (trios K to N), the mt lineage of the hatched juveniles is the one that characterize the trio, since the adults are of the same mt lineage. (2) For the hatched juveniles from trios composed of L1/L2/L2 and L1/L3/L3 adults (trios O to V corresponding to a mix between the two colour morphs), we combined High‐Resolution Melting (HRM) analysis to DNA Barcoding (Bar‐HRM), as described by Baudrin et al. (2020). This method was shown to discriminate L1 and L2/L3 lineages but not L2 and L3 (Baudrin et al., 2020). In brief, the HRM analysis of the putative parents was carried out in triplicates, at the same time as for the juveniles, using MeltDoctor HRM Master Mix (Applied Biosystems) according to the manufacturer protocol in 10 μL reaction volume and using EwD/EwE primers developed by Bienert et al. (2012). (3) For the hatched juveniles produced from the remaining trios, composed of a mix of L2 and L3 adults of the green morph, we carried out a PCR screening using COI species‐specific primers as described by King et al. (2010). The COI‐AchL2A‐F5 + COI‐AchL2A‐R3, COI‐AchL2B‐F3 + COI‐AchL2B‐R3, COI‐AchL3‐F2 + COI‐AchL3‐R2 and COI‐Ach L1–F4 and COI‐Ach L1–R2 were first used to check the method in triplicates on all putative parents. This method allowed us to assign 100% of the putative parents to their correct COI lineage (previously determined by sequencing). It also revealed that all L2 parents could be amplified using L2 B‐specific primers, making the use of L2A‐specific primers unnecessary. We amplified the DNA of the hatched juveniles in monoplex using two couples of primers, according to the lineage of their putative parents (L2A and L3, L2A and L1 or L3 and L1). Amplifications were performed in 15 µl, containing 5× Flexi Reaction Buffer, 0.125 mM of each dNTP, 1.5 mM MgCl2, 0.5 U of GoTaq® Flexi DNA Polymerase (Promega), 0.2 µM of each primer and 1 µl of extracted DNA. PCR cycling conditions were 94°C for 3 min, followed by 35 cycles at 94°C for 30 s, 57°C (L1, L2A and L3 primers) or 58°C (for L2B primers) for 45 s, 72°C for 60 s and a final extension at 72°C for 10 min. We visualized the results on 2% agarose gels. The method used to determine the COI lineage of each hatched juvenile is available in the Supporting information. Finally, we controlled the assignments by sequencing (Supporting information) the COI fragment of a subset of 84 hatched juveniles, using the method described above.

Parentage assignment of the offspring from the cross‐breeding experiment using microsatellite markers

We genotyped 453 of the 681 hatched juveniles using four microsatellite markers, as described above. For the two trios F and M that produced few juveniles, we genotyped all the juveniles from the trio F (N juveniles produced = 6) and 14 juveniles of the trio M (N juveniles produced = 15). For trios that produced at least 20 juveniles, we genotyped on average 21.7 ± 1.6 juveniles per trio. Detailed information on the number of juveniles produced and genotyped in each trio is presented in Table 1. Because genotyping errors are an important source of problems for parentage analysis, we duplicated all the PCR results of the adults that were putative parents in the cross‐breeding experiments.

The CERVUS v.3.0.7 software was used (i) to calculate the combined non‐exclusion probabilities over loci, for first parent (NE‐1P), second parent (NE‐2P), parent pair (NE‐PP), unrelated individuals (NE‐I) or siblings (NE‐SI), (ii) to determine the confidence of parentage assignments using a likelihood‐based approach to assign parental origin combined with simulation parentage analysis and (iii) to perform parentage assignment. Likelihood score ratios (LOD) estimate the likelihood that the candidate parent is the true parent divided by the likelihood that the candidate parent is not the true parent. Before proceeding to the parentage assignment, simulations were run in CERVUS to determine the distribution of the critical values of LOD scores for 80% and 95% confidence levels. The following simulation parameters for 100 000 offspring were chosen: ‘candidate parents’ 3, ‘prop. Sampled’ 1, ‘prop. loci typed’ 0.85 and ‘prop. loci mistyped’ 0.01. Confidence levels obtained from simulations were used for true paternity screening of the offspring.

Reproductive strategies

To investigate differences in mating strategies between divergent mt lineages which could not be differentiated by nu markers, that is mating between the L2 and L3 lineages of the green morph, we studied reproduction first amongst the types of trios and then at the level of each parent. These analyses focused on data from trios A to N. We used linear models to investigate differences in cocoon production, hatching rates and the number of hatched juveniles amongst the following trios: L2/L2/L3 (trios A to E), L2/L2/L2 (K and L), L2/L3/L3 (F to J) and L3/L3/L3 (M and N). At the level of each potential parent (of either L2 or L3 lineage), regardless of the trio in which they were involved, we studied the number of juveniles produced by each earthworm and tested for differences according to the lineage of the adults and the number of mates, using a linear model. We also tested whether these lineages consistently preferred to use a function (male or female) in such cross‐breeding. This was modelled as a binomial response with, for each parent, a two‐vector variable with the number of offspring produced using the female function and the number of offspring produced using the male function. Note that these analyses were limited to the subset of juveniles we genotyped (see above).

We followed the same procedure to investigate mating strategies between lineages that diverged both at the mt and nu levels, that is the reproductive isolation between L2 or L3 (of the green morph) and L1 (of the pink morph) lineages. First, we used linear models to test for differences in cocoon production, hatching rates and the number of hatched juveniles, between trios including one L1 individual (trios L1/L2/L2, O to R, and trios L1/L3/L3, S to V) to trios with no L1 (trios A to N). Then, focusing on trios including a L1 parent (trios O to V), we studied the parentage assignment of the juveniles. We specifically tested for differences in the number of juveniles produced from breeding between adults of L2 or L3, to the number of juveniles produced from cross‐breeding between an adult of L2 or L3 and an adult of L1, using a linear model. Last, based on the subset of juveniles produced from the crossing between adults of L2 or L3 and L1, we tested whether these lineages consistently preferred to use a function (male or female) in such cross‐breeding using a Pearson's Chi‐square test.

All analyses were performed in R version 4.0.2 (R Core Team, 2020).

RESULTS

Genetic structure of the field population

We examined and compared the population genetic structure of the A. chlorotica agg individuals sampled in the field using their individual assignment based on the mitochondrial COI and the nuclear microsatellites.

Mitochondrial genetic variation

We obtained a total of 154 COI sequences from the 157 juveniles and adults sampled in the same habitat. Amongst them, we identified 20 haplotypes (Genbank accessions MZ930411‐30) amongst which 9, 5, 4 and 2 haplotypes belonged, respectively, to the L1, L2, L3 and L4 mitochondrial lineages (Figure 1). No haplotype was identified as belonging to the L5 lineage. In the field population, the majority of the individuals belonged to the L2 (50.65%) and the L3 (35.71%) lineages, whereas 10.39% belonged to the L1 lineage and only 3.25% to the L4 lineage (Figure 2).

FIGURE 1

Haplotype network showing the frequency of each COI haplotype belonging to the L1 (H1–H9), L2 (H10–H14), L3 (H15–H18) and L4 (H19–H20) mitochondrial lineages of the Allolobophora chlorotica complex in the Preston population and their relationships. The small grey circles indicate inferred steps not found in the data set. Connecting lines show mutational pathways between haplotypes. More than 50 mutational steps are indicated

FIGURE 2

Nuclear clusters and mitochondrial lineages for the field population. Each vertical bar represents one individual. The first row refers to the mitochondrial lineage, the second to the estimated nuclear composition based on the multilocus microsatellite genotype (STRUCTURE software)

Nuclear genetic variation

The statistics for the four microsatellite loci used in the study are given in Table 2. Overall polymorphism was high with the number of alleles (N all) varying from 10 to 25 alleles per locus, H e from 0.649 to 0.886 and PIC from 0.605 to 0.873. Null allele frequencies ranged from 0.0306 (Ac476) to 0.1178 (Ac418). The STRUCTURE analysis supported the presence of two nu clusters (A and B) within the data set composed of the 157 individuals sampled in the field (∆K = 501 for K = 2; Figure 2). Overall, the nu cluster A corresponds to the L1 mt lineage and the nu cluster B to the L2 and L3 mt lineages, with a few exceptions. Specifically, the case of the L4 mt lineage was ambiguous since amongst the five L4 individuals, two were assigned to the nu cluster A and three to the nu cluster B. This result might be explained by the lower number of microsatellite loci used here, by comparison with Dupont et al. (2016; 4 vs. 5 microsatellite loci). However, this has no consequence on the present study since no L4 individual was used in the cross‐breeding experiment. The association was also not categorical for five individuals from the L1, L2 and L3 lineages, whose mt haplotype did not correspond to their nu cluster. Such cases have been previously described as resulting from introgression (Dupont et al., 2016). Specifically, two individuals of the L1 mt lineage were assigned to the nu cluster B and three individuals belonging to the L2 mt lineage were assigned to the nu cluster A. Note that these latter three introgressed individuals were used in three different trios (D, E and O, Table 3). Last, three individuals of L2 and one individual of L3 showed admixture above 10%, suggesting recent hybridization. The L3 hybrid individual which we assigned at 60% to the nu cluster B and at 40% to the nu cluster A was used in trio B (Table 3).

TABLE 2

Genetic diversity parameters of the four microsatellite loci used in the study and combined non‐exclusion probabilities over loci (Combined NE) for first parent (NE‐1P), second parent (NE‐2P), parent pair (NE‐PP), unrelated individual (NE‐I) or siblings (NE‐SI)

Locus	N all	H o	H e	PIC	NE‐1P	NE‐2P	NE‐PP	NE‐I	NE‐SI	HW	F (Null)
Ac127	20	0.739	0.855	0.838	0.447	0.286	0.116	0.036	0.333	*	0.0752
Ac170	25	0.747	0.886	0.873	0.370	0.228	0.076	0.023	0.314	*	0.0817
Ac418	12	0.519	0.649	0.605	0.750	0.577	0.382	0.167	0.468	***	0.1178
Ac476	10	0.683	0.739	0.697	0.659	0.483	0.292	0.109	0.409	NS	0.0306
Mean – Combined NE	16.75	0.672	0.782	0.753	0.08178	0.01146	0.00099	0.00001	0.20046	ND	ND

*Significant at the 5% level; ***Significant at the 0.1% level.

Abbreviations: F (Null), estimated frequency of null alleles; H e, expected heterozygosity; H o, observed heterozygosity; HW, exact test of departure from Hardy–Weinberg equilibrium; N all, Number of alleles; ND, Not done; NS, not significant; PIC, polymorphic information content.

TABLE 3

Nuclear clusters and parentage assignments

Trio	Nuclear cluster of the adults	Cross without offspring			Individual without offspring
		No.	Lineages	Cluster	No.	Lineage	Cluster
A	3B	0	–	–	0	–	–
B	2B –1 hybrid 40%A	0	–	–	0	–	–
C	3B	1	L3–L2	B–B	0	–	–
D	2B – 1 introgressed 100% A	1	L2–L2	B–A a	0	–	–
E	2B –1 introgressed 100% A	1	L2–L2	B–A a	0	–	–
F	3B	2	L3–L2, L3–L3	B–B	1	L3	B
G	3B	0	–	–	0	–	–
H	3B	1	L3–L2	B–B	0	–	–
I	3B	0	–	–	0	–	–
J	3B	0	–	–	0	–	–
K	3B	0	–	–	0	–	–
L	3B	1	L2–L2	B–B	0	–	–
M	3B	2	L3–L3	B–B	1	L3	B
N	3B	0	–	–	0	–	–
O	1B – 1 introgressed 100% A – 1A	1	L2–L2	B–A a	0	–	–
P	2B – 1A	1	L1–L2	A–B	0	–	–
Q	2B – 1A	2	L1–L2	A–B	1	L1	A
R	2B – 1A	2	L1–L2	A–B	1	L1	A
S	2B – 1A	2	L1–L3	A–B	1	L1	A
T	2B – 1A	0	–	–	–	–	–
U	2B – 1A	2	L1–L3	A–B	1	L1	A
V	2B – 1A	2	L1–L3	A–B	1	L1	A

For each trio: the assemblage of adults according to their nuclear cluster such as inferred from the STRUCTURE analysis, the number (No.) of pairs of adults which produced no juveniles and the mitochondrial lineage and the nuclear cluster of each potential parent, and, the number, lineage and cluster of the adults that produced no juveniles.

Indicate the introgressed adults.

Parentage assignment

The simulations in CERVUS resulted in high assignment rates, of 95% and 98% of parental pair assignment, at strict (95%) and relaxed (80%) levels, respectively. Only 2% of simulated offspring remained unassigned. Values of combined non‐exclusion probabilities for the first parent, second parent and parent pairs were low, at 0.08178, 0.01815 and 0.00001, respectively. Results per trio are summarized in Table 3 and detailed in the Supporting information. Most of the adults from the trios had offspring, except two L3 individuals in the trios F (L2/L3/L3) and M (L3/L3/L3) and five L1 individuals in the trios Q and R (L1/L2/L2) and S, U and V (L1/L3/L3). In several trios (A, B, G, I, J, K, N, T), all possible crosses produced offspring, highlighting that the individuals frequently mated with the two available partners. It is worth noting that individuals who only had viable offspring with one partner could have reproduced and had cocoons (but not viable) with both. Moreover, in the three trios comprising an introgressed individual, one cross in each trio involving this individual did not produce any viable offspring (Table 3).

Overall, trios produced 45.32 ± 3.41 cocoons (mean ± SE), of which 30.95 ± 2.54 juveniles hatched (mean ± SE). The number of cocoons that hatched per trio was proportional to the total number of cocoons it produced (estimate = 0.64 ± 0.08, t = 7.60, p < 0.001, R 2 = 0.73, Table 1). Based on the parentage assignment of the subset of juveniles genotyped, we noted that seven adults did not produce any juveniles (Table 3).

To investigate differences in mating strategies between divergent mt lineages which could not be differentiated by nu markers, we focused on trios A to N, which involved adults of the L2 and L3 lineages only. We looked at the reproduction amongst trios and at the individual level of the parent. Amongst trios, we found that trios with a majority of L2 adults (L2/L2/L3) produced significantly more juveniles than the trios with a majority of L3 adults (L2/L3/L3) (F 1,12 = 5.42, p = 0.038, Figure 3). The L2/L3/L3 trios produced a significantly lower number of cocoons (estimated difference ± SE = −17.00 ± 7.11, t = −2.39, p = 0.034, Figure 3), but they showed a similar rate of hatching than the L2/L2/L3 trios (estimated difference ± SE = −0.04 ± 0.06, t = −0.69, p = 0.51, Figure 3). At the individual level of the parent, that is regardless of the trio in which they were involved, we found no difference in the number of juveniles produced between adults of the L2 and L3 lineages (F 1,38 = 0.35, p = 0.56). However, the adults who mated with both their potential mates produced significantly more juveniles (F 2,38 = 10.51, p < 0.001). Adults who mated with a unique partner produced 10.43 ± 1.37 juveniles (mean ± SE), whereas adults who mated with both their potential mates produced 15.23 ± 1.70 juveniles (mean ± SE). Note that the number of partners did not vary between lineages (χ2 = 3.30, df = 2, p = 0.19). In such cases of cross‐breeding between L2 and L3 lineages, the L2 parent preferentially used their male function, whereas the L3 parent used their female function (χ2 = 6.84, df = 1, p = 0.009).

FIGURE 3

Boxplot of (a) the total number of cocoons produced, (b) the hatching rate and (c) the number of juveniles, according to each type of trio. In hatched boxes, the values are for the mixed trios of lineages. In light and dark green are the trios of L2 and L3 (of green morph) and in pink are the trios including an adult of L1 (of the pink morph)

To investigate mating strategies between lineages that diverged both at the mt and nu levels, we focused on the eight trios including one L1 (trios O to V). When comparing the reproductive rates of trios including one L1 (trios O to V) to trios with no L1 (trios A to N), we observed that trios that involved one L1 adult produced significantly more cocoons (estimate ± SE = 14.04 ± 6.55, t = 2.14, p = 0.04), but their hatching rate was significantly lower (estimate ± SE = −0.15 ± 0.045, t = −3.44, p < 0.01), leading to similar numbers of juveniles than in trios without L1 (estimate ± SE = 1.84 ± 5.40, t = 0.34, p = 0.74, Figure 3). The parentage assignment of juveniles in trios including L1 further showed that the cocoons that hatched were mainly from mating between adults of L2 or L3 (at 81.3%, Table 3). In five of the trios (Q, R, S, U, V), the L1 did not produce any viable juveniles. The number of juveniles produced by adults of L2 and L3 was significantly higher when they reproduced with an adult of the same lineage than with an L1 adult (estimate ± SE = 16.38 ± 2.22, t = 7.39, p < 0.001). It is difficult to draw any solid conclusion about sexual function preferentially used in crossing between L1 and L2 or L3 because only three of the eight L1 adults produced juveniles from such crossing. Still, based on the available data, 28 of the 34 juveniles were produced using the female function, suggesting that in such crossing L1 adults would preferentially reproduce using their female function (χ2 = 9.90, df = 2, p = 0.007).

Note that all our results hold, whether we considered that the introgressed individuals belong to the lineage defined by their mt haplotype or by their nu genotype.

DISCUSSION

The genetic study that we carried out here, on individuals belonging to different mt lineages but sampled in the same habitat (i.e. in syntopy) and on their offspring, sheds new light to understand reproductive isolation patterns within the A. chlorotica agg. On one hand, we confirmed that individuals from the two divergent mt lineages L2 and L3 cannot be distinguished using nu markers, and on the other hand, we revealed that the sex function used in the L2–L3 crosses was specific to each lineage. Moreover, we provide additional evidence that mechanisms of post‐zygotic reproductive isolation between different but closely related species (L1 and L2/ L3) are at play. Note that our interpretations of the reproductive isolation patterns are in part restricted to the subset of juveniles that have hatched and for which we know the parentage. Still, potential interpretation biases are limited since the juveniles genotyped from each trio were selected at random and that the proportion of genotyped juveniles per trio was relatively high (73.3 ± 20.4 juveniles ranging from 42.9% to 100%). We are, therefore, confident that what was observed is representative of what could have been observed from the genotyping of all individuals. Another limitation of this study, due to the effort required to set up these cross‐experiments, is the low replication of the type of trios, with N = 4 or 5 for between lineages trios and N = 2 for within lineages trio. Despite these limitations, we were able to produce the first results regarding patterns of reproductive isolation within a species complex exhibiting cases of introgression and mitonuclear discordance.

Reproductive strategy within the L2/L3 lineages of the Allolobophora chlorotica agg

The parentage analysis revealed that, over a 4‐month period, A. chlorotica adults can have offspring from multiple partners. Although multiple mating was known in several earthworm species, e.g. in Eisenia fetida (Monroy et al., 2003), Lumbricus terrestris (Butt & Nuutinen, 1998; Michiels et al., 2001) and Hormogaster elisae (Novo et al., 2013), to the best of our knowledge, only Novo et al. (2013) used similar molecular techniques to examine the fate of sperm after its transfer to a mate. Using microsatellite markers, Novo et al. (2013) found cases of multiple paternity in Hormogaster elisae and showed that paternity is influenced by the order of copulation. Here, we showed that A. chlorotica adults that mated with both their potential mates in the trio produced significantly more juveniles. This result is consistent with the hypothesis of Porto et al. (2012) who investigated the effect of multiple mating on female reproduction in the earthworm Eisenia andrei. These authors found that multiple mating was beneficial for female reproduction, with increased hatching success of the cocoons, and they suggested that this may result from an increase in sperm quantity and/or diversity in the spermathecae of multiple‐mated earthworms.

Further, we found that the sex function used in the L2–L3 crosses was specific to each lineage, with the L2 parents that preferentially used their male function and L3 parents their female function. These differences might be the result of a conflict of interest between the mating partners (i.e. the sperm donor and the sperm recipient, Schärer, 2009). It was suggested that sperm donors develop traits that either directly boost the female function of the recipient or disrupt its male function (Schärer, 2009). For instance, in the earthworm Lumbricus terrestris, the injection of a substance from its setal glands through copulatory setae during mating increases sperm uptake and delays re‐mating (Koene et al., 2005). By analogy, it could be that during reciprocal copulations and sperm transfers, the L2 individual could transfer a substance that would increase female allocation in the L3 recipient. Understanding, on a general basis, these postcopulatory mechanisms of resource allocation between sexual functions and how they may interfere with multiple mating could extend further our understanding of the sexual selection mechanisms at play in this species. Multiple mating with different partners as well as sperm storage in spermathecae (three pairs of spermathecae have been described) reported in the morphospecies A. chlorotica by Sims and Gerard (1999) further advocates for interferences between multiple mating and post‐copulatory sexual selection.

It is noteworthy that the L2 mt‐lineage was more frequent than that of the L3 (50.6% of L2 vs. 35.7% of L3) in the prospected farmer's field when we could expect the opposite since in L2/L3 crosses, it was the L3 parents that most frequently used their female function, thereby transmitting their mt genome to the offspring. An explanation could be that L2/L2 crosses are more frequent that L2/L3 crosses in the field. An alternative hypothesis would be that L2 individuals have an adaptive advantage that makes them disproportionately more abundant.

Reproductive isolation between L1 and L2/L3 lineages of the Allolobophora chlorotica agg

Our results support the observations from previous breeding experiments, according to which post‐zygotic reproductive isolation processes are occurring between L1 and L2/L3 lineages (Lowe & Butt, 2008). Indeed, we found that the high production of cocoons in the L1/L2/L2 and L1/L3/L3 trios was counterbalanced by a high mortality, as evidenced by their low hatching rate (Table 1). Likewise, Lowe and Butt (2008) found that crosses between adults of the green (L2 and L3) and pink (L1 and L4) morphs were able to reproduce, but the viability of their cocoons was substantially reduced, with a hatching rate of 6% and 59% for the green and pink morphs, respectively. It would have been particularly interesting to determine the parents of the cocoons which did not hatch by genotyping. This was attempted, but we were unable to distinguish between the maternal and juvenile DNA in the DNA extracts from the cocoons. Consistently with the observed reproductive isolation, 63% of the L1 adults produced no juveniles and, comparatively, the L2/L3 adults produced significantly more juveniles when reproducing with an adult of the same lineage that with an adult of L1.

Despite evidence of reproductive isolation, we still obtained a few juvenile hybrids of L2/L3 and L1 (in trio O, P and T) from our cross‐breeding experiment. Hybrids were also recorded, anecdotally, amongst individuals from the field population (Figure 2), suggesting that some can reach adulthood. One of these hybrids, with an L3 mitochondrial lineage but assigned at 60% to cluster B and at 40% to cluster A, was used in trio B. The parentage analysis revealed that this adult had produced 10 juveniles using only its female function. Sterility of the male function in hybrids has already been proposed from the breeding experiment by Lowe and Butt (2008) and deserves to be further investigated in order to test for Haldane's rule (i.e. disproportionate hybrid dysfunction in the male function, Haldane, 1922) in this simultaneous hermaphrodite.

Hybrid male sterility would allow an explanation for the existence of introgressed individuals that presented an L2 or L3 mt lineage but that were assigned to the nu cluster A (Figure 2, Table 3, Dupont et al., 2016). The case of such introgressed individuals, collected from field populations and used in three different trios (D, E and O, Table 3), was the most interesting. They were all from the L2 mt lineage but from the nu cluster A. In the trios D and E, composed of L2 and L3 adults with a majority of L2, one of the introgressed individuals produced one juvenile using its female function (in trio D), the other produced seven juveniles using both functions (in trio E: four and three juveniles produced using, respectively, its male and female function). In trio O, composed of one adult of L1 with two adults of L2, the introgressed individual reproduced using both functions (two and four juveniles produced using, respectively, its male and female function, Supporting information) but only with the L1 adult of the same nu cluster (A). Such a small number of juveniles produced from an introgressed parent do not allow any strong conclusion but still suggests that male function could be restored in introgressed individuals, whilst it is not active in hybrids. This advocates that, in hybrids, male‐specific genes are likely to be non‐ or misexpressed, whereas, in introgressed individuals, in which a full nu genome of one species has been reconstructed after successive backcrosses, the expression of these genes could be restored.

CONCLUSION

Overall, our results revealed different degrees of reproductive isolation within a complex of simultaneous hermaphrodite species that have reached different stages in the speciation process, resulting in opposite patterns of mitonuclear congruence. First, we confirmed post‐zygotic reproductive isolation processes and hybrid male sterility for the crossing between lineages that diverged both at the mt and nu levels (i.e mitonuclear congruence). Furthermore, the existence of two mt lineages that cannot be distinguished with nu markers (i.e. mitonuclear discordance) could be, in the first place, explained by a lack of reproductive isolation mechanisms that have evolved during the allopatric phase, which was likely at the origin of the mt divergence. However, despite the apparent gene flow amongst these two mt lineages, we revealed that they present different reproductive strategies, with one lineage using more frequently its female function in cross‐breeding. Basic sex allocation models for simultaneous hermaphrodites generally assume a linear trade‐off between the allocation to male and female functions, so that higher allocation to one function leads to a proportional decrease in the allocation to the other function (Schärer, 2009). It is expected that individuals show a preference for adopting the sex role that tends to offer the higher potential fitness gain per mating (Anthes et al., 2006). Thus, the difference in the use of sex function could result from differential sex allocation mechanisms that have evolved during the evolutionary history of this species complex.

AUTHOR CONTRIBUTIONS

L.D., D.P. and K.R.B. conceived and designed the experiment. K.R.B. and L.D. performed the experiments. L.D. and H.A. conducted the analyses and wrote the first draft of the manuscript. All authors contributed to revising and editing the manuscript.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.14017.

Supplementary Material

Click here for additional data file.