Anna Anastaseni

Antonio Romano

A perceptual experiment testing the Italian /ʎ/ -/j/ contrast

Keywords: Auditory phonetics, identification task, palatal lateral, palatal approximant, duration

Due to the articulatory proximity, Italian /ʎ/ and /j/ are confused by some speakers. The neutralisation of the /ʎ/-/j/ contrast lets prevail [j]-like sounds; /lj/ and /llj/ are sometime realised with a [ʎ]-like sounds.

Two classification tasks were proposed to a sample of 113 native Italian listeners to assess the incidence of the duration factor in distinguishing intervocalic segments [ʎ], [j] and [lj], manipulating the duration of the sounds within the same contexts.

It emerged that participants were always able to tell apart [lj]-like stimuli, while they showed more variation in the classification of /ʎ/ and /j/ categories. The duration plays a significant role in the identification of /ʎ/ and /j/: with less than 140-160 ms, the perceived stimulus originated from a [ʎ]-sound is generally classified as /j/, whereas stimuli with increasing durations of an original [j]-sound need an upper threshold (~200 ms) to be mostly attributed to /ʎ/.

INTRODUCTION

The phonological inventory of Italian distinguishes between the lateral palatal /ʎ(ʎ)/ and the approximant palatal /j/; however, this distinction has low-functional load. In neutral pronunciation the palatal /ʎ/ is always geminate in postvocalic position [1]. The alternation as short and long is governed by position: is not phonemic. Its duration exceeds 100 ms in standard Italian and in several regional varieties [2], but it does not reach the length of the corresponding non-palatal lexical geminate [lː] and, apparently, it does not induce the preceding vowel to shorten [START_REF] Endo | Caratteristiche prosodiche delle così dette 'rafforzate' italiane[END_REF]. On the other hand [j] can be considered as inherently short.

The lateral one has a complex articulatory mechanism ([1], [START_REF] Bladon | Lateral consonants in Italian[END_REF], [START_REF] Recasens | Vowel-to-vowel coarticulation in Catalan VCV sequences: an articulatory and acoustical study[END_REF]). It appears late in the phonological inventory of Italian-speaking children [START_REF] Zanobini | Phonology and Language Development in Italian Children: An Analysis of Production and Accuracy[END_REF] and it is often simplified ( [START_REF] Galatà | Primi dati sull'acquisizione fonetico-fonologica dell'italiano L2 in prescolari rumeni[END_REF], [START_REF] Zampaulo | The evolution of the (alveolo) palatal lateral consonant in Spanish and Portuguese[END_REF]). The phoneme /ʎ/ is rare in majority of the world's languages; it is present in only 20 of the 451 languages considered in UPSID (UCLA Phonological Segment Inventory Database). Moreover, it occurs in only 16 languages whose inventory also includes the palatal approximant /j/ (by contrast, in 378 languages in the database).

The [ʎ]-sound and [j]-sound can be defined as approximants: they have a periodic structure, but a reduced intensity of about 10 dB compared to vowels [START_REF] Pettorino | Indagine acustica sulle approssimanti dell'italiano[END_REF]. The neutralisation of the opposition between /ʎ/ and /j/ is widespread and is favoured by acoustic similarity, as shown by the evolution of phonological systems in other romance languages ( [START_REF] Zampaulo | The evolution of the (alveolo) palatal lateral consonant in Spanish and Portuguese[END_REF], [START_REF] Porras | Spanish Yeísmo: A Cognitive Linguistic Approach to Phonological Change[END_REF], [START_REF] Rost Bagudanch | Una panorámica del yeísmo:?` un proceso acabado o en construcción? Una panorámica del yeísmo:?[END_REF]). In Italy, the delateralisation of /ʎ/ (/ʎ/ → [j]) is a widespread phenomenon in the Centre-South (with considerable micro-areal variation, [START_REF] Romano | An MRI Study on the Articulatory Properties of Italian Consonants[END_REF]- [START_REF] Marotta | La liquida laterale nel livornese[END_REF]), with the exception of Tuscany and Sardinia [1]. The neutralisation of the /ʎ/ -/j/ contrast allows [j]-like sounds to prevail. On the other hand, in the same regions a lengthening of /j/ in intervocalic position is observed ([1], [START_REF] Marotta | La liquida laterale nel livornese[END_REF]). Furthermore, the phonological status of /ʎ/, mainly in Northern regions, has to be assessed taking into account the extension of [ʎ]-like sounds to the realisation of /lj/ and /llj/ [1], [START_REF] Romano | Argomenti scelti di glottologia e linguistica[END_REF], [START_REF] Bertinetto | The sound pattern of Standard Italian, as compared with the varieties spoken in Florence, Milan and Rome[END_REF].

Moreover, in pronunciation, among speakers of younger generations especially in the North, regional features are weak and often phonetic features coming from different regional Italian varieties are co-present [START_REF] Cerruti | Varietà dell'italiano[END_REF]. As the three corpora ParlaTO, KIP and CLIPS show, in the variety of Italian spoken in Turin, [ʎj] and [jː] are also attested, not only the expected variants [ʎ] and [lj] [START_REF] Fontanot | Aspetti di italiano "composito" nella pronuncia di giovani torinesi: un'indagine su corpora[END_REF].

Within the community we have been observing in the last two years by means of acoustical analysis, the loss of the lateral feature seems not to give rise to a real neutralisation between /ʎ/ and /j/ since their common realisations are respectively [jː] and [j] and a length contrast is often preserved.

A perceptual experiment has been set up in order to verify the hypothesis that a duration feature is sufficient to avoid a full merger. It is important to check the ability of Italian speakers to classify different phonemes. In this pilot study, we have prepared two identification tasks.

Firstly, stimuli with /ʎ/, /j/, /lj/ and /llj/ were presented to the listeners. The results showed that participants had more difficulty distinguishing the palatal lateral and palatal approximant. For this reason, we decided to focus our attention on /ʎ/ and /j/, with an exploratory reanalysis of the first experiment.

METHOD

In both experiments, an identification task was proposed to a group of Italian speakers.

Participants had to listen to a list of pseudo-words and had to associate each pseudo-word heard to an orthographic transcription within a closed list of possible labels. The orthographic labels were spelled with ‹gl›, which in Italian is usually associated with /ʎ/, ‹i›, associated with /j/, ‹li›, associated with /lj/, and ‹lli› associated with /llj/ in the first experiment; and in the second experiment with ‹gl› and ‹i›.

The tasks were disseminated online via the Folerpa platform [START_REF] Fernández Rei | FOLERPA: Ferramenta On-Line para ExpeRimentación PerceptivA , Santiago de Compostela[END_REF], an online tool developed at the Instituto da Lingua Galega (ILGA) of the University of Santiago de Compostela (Galicia, Spain).

In the first task 32 pseudo-words were presented to each participant, in the second task 40; each item was tested twice (total stimuli: 64 for the first test; 80 for the second). To avoid list effects the stimuli were presented at each replay in random order. In the middle of the task there was a pause.

The response times were calculated in milliseconds, between the end of the repetition of each stimulus and the listener's response.

Participants

The first task was proposed to 67 university students. Age ranged from 18 years to 29 years (mean: 21.20; SD: 2.52), the participants were Italian native speakers. 62% of the subjects indicated Piedmont as their region of origin.

On the other hand, the second experiment was carried out with 46 students of the University of Turin aged between 21 and 31 (mean 23.93; SD: 2.07). 57% of the participants indicated Piedmont as their region of origin.

Participants were asked to indicate all the languages in which they have active and passive competence. The mother tongue and the L2, in fact, can modify the performance of identity task. Comparative studies conducted on Catalan and Spanish show that Catalan speakers have better performance in distinguishing [ʎː] and [j] than Castilian speakers [START_REF] Bagudanch | La percepción de/ʎ/ y /j/en catalán y en español. Implicaciones en la explicación del yeísmo[END_REF], [START_REF] Rost Bagudanch | Bilingualism and sound change: perception in the /ʎ/-/j/ merger process in Majorcan Spanish[END_REF].

All participants had no hearing problems, nor had they ever suspected hearing loss.

Finally, it was controlled that no subject had given more than 85% of identical answers and that there weren't any technical problems.

Each student was asked to complete the Vinegrad Plus Questionnaire [START_REF] Montesano | SLD: promoting continuity between School and University education. Results of a screening program in Calabria. Form@ re-Open[END_REF], according to the modified version of Montesano and Valenti [23]. This test requests answers to a variety of questions designed to identify common traits of dyslexia. 33% of the participants were in the risk range, while 28% said they suspected to be dyslexic. These data are particularly important because phonological deficits are linked with dyslexia. In fact, many studies have shown the presence of phonological difficulties even in dyslexic adults, especially in phonological tasks that require access to phonology from visual inputs [START_REF] Svensson | How persistent are phonological difficulties? A longitudinal study of reading retarded children[END_REF][START_REF] Danelli | How many deficits in the same dyslexic brains? A behavioural and fMRI assessment of comorbidity in adult dyslexics[END_REF][START_REF] Mccrory | More than words: a common neural basis for reading and naming deficits in developmental dyslexia?[END_REF].

This fact is relevant since the identification task in the proposed modalities has an important visual component. Participants, in fact, had to make an effort of grapho-phonematic conversion to associate the heard stimuli with the proposed orthographic labels ([27]- [START_REF] Stoehr | The impact of orthographic forms on speech production and perception: An artificial vowel-learning study[END_REF]).

Stimuli

The stimuli consisted of four groups of pseudowords with [ʎː], [j] [lj] and [lːj] (Table 1).

All stimuli had the same phonotactic structure. They were disyllabic and paroxytone pseudo-words, with the target sound in intervocalic position. The vowels preceding it were stressed.

Pseudo-words were chosen to avoid the influence of lexical variables of existing words on the identification test. Once the pseudo-words had been selected for the test, a professional speaker was recorded in a silent room, he repeated the stimuli in isolation. Afterwards, the durations of the vowels preceding the target phonics [ʎː], [lːj], on the one hand, and [j],

[lj], on the other hand, were equalised to 145 ms (± 5ms) and 160 (± 5 ms) using GoldWave audio editing software.

The mean of the duration of [ʎː] in the four pseudo-words (ciaglia, naglio, roglio, toglia) was 158 ms. All stimuli were standardised to this value (± 2 ms), and then, keeping the previous vowel identical, the duration of the target sound was varied by reducing it by 40 ms each time and creating four new stimuli with [ʎ] durations of 220 ms, 180 ms, 140 ms and 100 ms.

During the first task, all stimuli were used, whereas for the second task, only stimuli with [ʎː] and [j] were used.

Due to the impact of the duration factor for the classification of [ʎ] ( §2.3.1), on the second task we focused on the rule of duration on [ʎ]-[j] distinction. The same manipulation was done with the stimuli presenting [j]-sound, in the original audio of 180 ms duration, manipulated to reproduce the same continuum as [ʎ] (260 ms, 220 ms, 180 ms, 140 ms, 100 ms).

Data Analysis

Accuracy of participants' responses in the dyslexia risk range

As mentioned, the Vinegrad Plus was administered to the participants of the two experiments in order to check whether there was a percentage of people in the sample who were in the dyslexia risk group. It emerged that 33% were actually in this group. Given the possible correlation between dyslexia and the outcome of the tasks, the data on the accuracy of the answers were cross-referenced with the presence or absence of dyslexia risk.

To verify this hypothesis, the Chi Square test was applied and the χ² value and Cramer's V were calculated. This analysis did not yield significant results (χ² = 0.626; V = 0.013). This fact suggests that difficulties in the classification of /ʎ/ and /j/ are not characteristic only of those with specific difficulties, but more widespread in the population.

Because of dyslexia is not affecting the results of the tasks we proposed, we did not delete the data of any participant for analysis.

Experiment 1

The first task showed that the participants were able to correctly associate /lj/ and /llj/ to the corresponding orthographic label ‹li› and ‹lli› (see The data were processed statistically with SPSS for Analytics (by IBM), v. 28. As we were working with categorical data, crosstabs were used to check the relationship between two variables (type of stimulus and listener categorisation). We conducted Chi-square test and Cramer V. The significance value was set at 0.05 (α=0.5).

As a result, the degree of association of /lj/ and /llj/ with the respective graphic labels ‹li› and ‹lli› is significantly higher (χ²= 606.696; V=0.871) than for /ʎ/ and /j/ (χ²= 191.701; V= 0.490).

In addition, we assessed the participants' sensitivity to the duration of [ʎ] by evaluating whether it affects classification. As discussed, in fact, duration is one of the variables of the realisation of these phonemes. For this reason, each pseudo-word was proposed with the target sound [ʎ] having durations of 260 ms, 220 ms, 180 ms, 140 ms and 100 ms.

It was found that the duration factor played a significant role in the correct recognition of the sound (Pearson's r=0.947**, p=0.01). As the duration of [ʎ] decreased, the number of target responses ‹gl› dropped from 87% at 260 ms to just 9% at 100 ms (Figure 1). Given these initial results, it was decided to investigate this further by proposing a task in which not only the length of [ʎ] but also that of [j] was manipulated (Experiment 2).

Experiment 2

Our results of the first experiment demonstrate that the duration of [ʎ] is decisive in the subjects' choices in the identification task. The data collected in the first task was then confirmed by the second experiment and are presented below. In this experiment, listeners had to associate manipulated stimuli with /ʎ/ and /j/ (100-260 ms) and orthographic labels with ‹gl› and ‹i›.

The results show that in more than 90% of the responses the sound [ʎ]≥220 ms it is solidly associated with the target transcription ‹gl› and the stimuli with [j] at 100 ms is associated with ‹i›.

There is a directly proportional relationship between the duration of [ʎ] and the target association ‹gl› (Pearson's r = 0.935**; p = 0.02): the longer the duration of the sound, the better its recognition.

Whereas, surprisingly, an inversely proportional relationship appears in the recognition of [j] (Pearson's r = -0.988; p = 0.002; Figure 22). The [ʎ]-sound below the 180 ms threshold is no longer identified with certainty. On the other hand, [j] is almost unambiguously associated with ‹i› when it has durations of 100 ms and 140 ms. Significantly, however, the curves describing the classification of [ʎ] and [j] as ‹gl› and ‹i› according to duration do not show an entirely symmetrical profile.

The most significant finding concerns pseudowords with [j] of longer duration, which are alternately associated with ‹gl› and ‹i›.

It appears that different parameters come into conflict: the quality of the sound and its duration. It follows that [j] of duration superior than 180 ms is clearly classifiable neither as /ʎ/ nor as /j/.

This finding suggests that the duration feature is regarded as a significant index by the listeners.

Indeed, in Italian, duration characterises and distinguishes the two phonemes.

CONCLUSIONS AND PERSPECTIVES

Two perceptual identification tasks were submitted to a total of 113 Italian-speaking university students.

The results of these tasks showed that participants were always able to tell apart [lj]-like stimuli, while they showed more variation in the way they classified the stimuli in the /ʎ/ and /j/ categories. In this case, it is confirmed that the duration plays a significant role in the identification of /ʎ/: with less than 140-160 ms. The perceived stimulus originated from a [ʎ]-sound is generally classified as /j/, whereas stimuli with increasing durations of an original [j]-sound need an upper threshold (around 200 ms) to be mostly attributed to /ʎ/.

The [ʎ]-sound is associated with longer duration than [j]. What has been observed has important consequences for spelling. Indeed, it must be taken into account that if at a perceptual level the distinction between phonemes is not clear, in writing the indirect route cannot be used or in any case, it has an increased chance of not being reliable.

Moreover, it is necessary to complete further research in order to reach more solid conclusions. It would therefore be necessary to verify the discrimination and identification capacity of the two phonemes in different age groups, as well as to extend the research geographically. This is the reason why we decided to collect data from a more controlled sample. Therefore, it was decided to propose a reading task, a pseudo-words dictation task and an auditory discrimination task to 60 Piedmontese speakers, divided into age groups [START_REF] Romano | Argomenti scelti di glottologia e linguistica[END_REF][START_REF] Bertinetto | The sound pattern of Standard Italian, as compared with the varieties spoken in Florence, Milan and Rome[END_REF][START_REF] Cerruti | Varietà dell'italiano[END_REF][START_REF] Fontanot | Aspetti di italiano "composito" nella pronuncia di giovani torinesi: un'indagine su corpora[END_REF][START_REF] Fernández Rei | FOLERPA: Ferramenta On-Line para ExpeRimentación PerceptivA , Santiago de Compostela[END_REF][START_REF] Bagudanch | La percepción de/ʎ/ y /j/en catalán y en español. Implicaciones en la explicación del yeísmo[END_REF](20-25](25-30). All data were collected in person. Data analysis is still in progress. Analisi del processo di retroflessione consonantica in area romanza, con dati sperimentali dal còrso e dal siciliano. Pisa: Scuola Normale Superiore.

Figure 1 :

 1 Figure 1: Classification of [ʎ] as a function of sound duration (in ms). Experiment 1.

[ 1 ]

 1 Canepari, L. 1999. Manuale di pronuncia italiana. Zanichelli. Bologna. [2] Celata, C. 2006 6 -PhD unpublished dissertation.

Figure 2 :

 2 Figure 2: Experiment 2. Classification of [ʎ] and [j] as a function of duration (in ms).

Table 1 :

 1 Selected stimuli for the identification tasks (task 1 and 2).

		Task 1, 2	Task 1
		/ʎ/	/j/	/lj/	/llj/
		‹gl›	‹i›	‹li›	‹lli›
	[ˈa]_[a]	/ˈtʃaʎa/	/ˈtʃaja/	/ˈtʃalja/	/ˈtʃallja/
		‹ciaglia›	‹ciaia›	‹cialia›	‹ciallia›
	[ˈa]_[o]	/ˈnaʎo/	/ˈnajo/	/ˈnaljo/	/ˈnalljo/
		‹naglio›	‹naio›	‹nalio›	‹nallio›
	[ˈɔ]_[o]	/ˈrɔʎo/	/ˈrɔjo/	/ˈrɔljo/	/ˈrɔlljo/
		‹roglio›	‹roio›	‹rolio›	‹rollio›
	[ˈɔ]_[a]	/ˈtɔʎa/	/ˈtɔja/	/ˈtɔlja/	/ˈtɔllja/
		‹toglia›	‹toia›	‹tolia›	‹tollia›

Table 2 )

 2 .

		‹gl›	‹i›	‹li›	‹lli›
	/ʎ/*	69%	22% 6%	4%
	/j/	23%	70% 3%	4%
	/lj/	2%	1%	91% 7%
	/llj/	3%	1%	5%	92%

Table 2 :

 2 Summary Here are the reported results of stimuli with a [ʎ] duration of 180 ms, the average value among those tested. By entering the duration values of [ʎ] =180 ms here, [ʎ] and [j] with the same duration are compared.

of identification task's results. Experiment 1. Classification of /ʎ/, /j/, /lj/, /llj/. *


Pierre Andre Eyer
email: pieyer@live.fr

Keywords: Population genetics, Reproductive isolation, Speciation continuum, Conditional parthenogenesis, Cuticular hydrocarbons, Sexual selection, Sibling species, Social insects

Integrative characterization of genetic and phenotypic differentiation in an ant species complex with strong hierarchical population structure and low dispersal abilities

Integrative characterisation of genetic and phenotypic differentiation in an ant species complex with strong hierarchical population structure and low dispersal abilities.

INTRODUCTION

Speciation is usually a long evolutionary process through which a population splits into different evolving lineages that ultimately give rise to sister species [START_REF] Queiroz | Species concepts and species delimitation[END_REF]. This segregation is usually triggered by selection and/or geographic isolation which, over time, leads to the accumulation of barriers to gene flow resulting in reproductive isolation [START_REF] Mayr | Systematics and the Origin of Species[END_REF][START_REF] Endler | Gene Flow and Population Differentiation: Studies of clines suggest that differentiation along environmental gradients may be independent of gene flow[END_REF][START_REF] Coyne | Speciation[END_REF][START_REF] Wang | Quantifying the roles of ecology and geography in spatial genetic divergence (F Adler, Ed.)[END_REF][START_REF] Sexton | Genetic isolation by environment or distance: which pattern of gene flow is most common?[END_REF][START_REF] Richardson | Microgeographic adaptation and the spatial scale of evolution[END_REF]. Before reaching complete isolation, divergent populations can harbour various levels of reproductive isolation that can be positioned on a speciation continuum [START_REF] Stankowski | Defining the speciation continuum[END_REF] and are often coupled with different levels of phenotypic and genetic differentiation [START_REF] Queiroz | Species concepts and species delimitation[END_REF][START_REF] Hendry | Along the speciation continuum in sticklebacks[END_REF]). Defining species boundaries along a speciation continuum is therefore complex and not always possible as it strongly relies on the criteria used to distinguish species. Integrative approaches relying on genetic and phenotypic differentiation associated to mating experiments directly assessing reproductive isolation should help distinguishing within-species structure from species boundaries [START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF].

By affecting gene flow, distinct life history traits, such as dispersal ability and reproductive systems, may profoundly affect the speciation process (e.g., [START_REF] Bolnick | The magnitude of local adaptation under genotype-dependent dispersal[END_REF][START_REF] Servedio | Geography, assortative mating, and the effects of sexual selection on speciation with gene flow[END_REF][START_REF] Slatkin | Isolation by distance in equilibrium and non-equilibrium populations[END_REF]. Low dispersal hence favors speciation events over restricted spatial scale [START_REF] Kisel | Speciation has a spatial scale that depends on levels of gene flow[END_REF] and enhances the establishment and persistence of within-species lineages that can ultimately lead to new species [START_REF] Dynesius | Persistence of within-species lineage: a neglected control of speciation rates[END_REF]. Species with weak dispersal abilities might therefore be particularly prone to the evolution of reproductive isolation over small geographic scale [START_REF] Pante | Species are hypotheses: avoid connectivity assessments based on pillars of sand[END_REF]. Similarly, the transition to asexual reproduction enhances speciation when asexual lineages become completely isolated from the sexual ancestors [START_REF] Neumann | Reproductive division of labour and thelytoky result in sympatric barriers to gene flow in honeybees (Apis mellifera L.)[END_REF][START_REF] Gebiola | Integration of molecular, ecological, morphological and endosymbiont data for species delimitation within the Pnigalio soemius complex (Hymenoptera: Eulophidae)[END_REF][START_REF] Confalonieri | Game of clones: Is Wolbachia inducing speciation in a weevil with a mixed reproductive mode?[END_REF]. Even when facultative, asexual reproduction might strongly impede gene flow. Species combining these two characteristics, limited dispersal and asexual reproduction, usually exhibit strong population structure. The use of multidisciplinary approaches is therefore particularly relevant to understand in which stage of the speciation continuum they are, as genetic data alone might not convey enough information.

In this study, we investigate genetic and phenotypic differentiation between and within two sibling species of the thermophilic ant genus Cataglyphis, known to have low dispersal abilities, and further coupled these analyses with behavioral experiments. The species status of C. cursor and C. piliscapa have long been debated (note that most of the publications based on localities of C. piliscapa wrongly used the species name C. cursor). Their identification is primary based on the number of hairs on the first tergite of the gaster [START_REF] Bondroit | Les fourmis de France et de Belgique[END_REF][START_REF] Agosti | A provisional list of the Balkan with a list to the worker caste. II. Key to the worker caste, including the European species without the Iberian[END_REF] and on cuticular hydrocarbon analyses [START_REF] Nowbahari | Individual, geographical and experimental variation of cuticular hydrocarbons of the ant Cataglyphis cursor (Hymenoptera: Formicidae): Their use in nest and subspecies recognition[END_REF]), but a robust quantitative assessment is still lacking. Interestingly, these two allopatric species are separated by the Rhône River in the Southern France, with C. cursor occurring on the East side and C. piliscapa on the West side [START_REF] Nowbahari | Individual, geographical and experimental variation of cuticular hydrocarbons of the ant Cataglyphis cursor (Hymenoptera: Formicidae): Their use in nest and subspecies recognition[END_REF]). In both species, female dispersal abilities are particularly weak, consistently with the low number of colonies produced each year and the widespread adoption of the fission process [START_REF] Lenoir | Reproduction and dispersal in the ant Cataglyphis Cursor (Hymenoptera, Formicidae)[END_REF][START_REF] Cronin | Recurrent evolution of dependent colony foundation across eusocial insects[END_REF]. During colony fission, young queens of both species leave the natal nest by foot with nestmate workers to found few new colonies nearby (four new nests at 7 m of distance on average) [START_REF] Lenoir | Reproduction and dispersal in the ant Cataglyphis Cursor (Hymenoptera, Formicidae)[END_REF][START_REF] Chéron | Unequal resource allocation among colonies produced by fission in the ant Cataglyphis cursor[END_REF]. Colonies of both species are headed by a single queen mated with several males [START_REF] Lenoir | Reproduction and dispersal in the ant Cataglyphis Cursor (Hymenoptera, Formicidae)[END_REF][START_REF] Pearcy | Conditional use of sex and parthenogenesis for worker and queen production in ants[END_REF]). Four main factors limit male gene flow in C. piliscapa : i) males are not good flyers and their effective dispersal distance was estimated to be less than 250 m [START_REF] Hardy | Small-scale spatial genetic structure in an ant species with sex-biased dispersal[END_REF]; ii) male gene flow is restricted to the pre-mating process since queen dispersal occurring after mating is limited to the worker's walking distance; iii) workers aggress males during mating, acting as a potential filtering process and indirect source of sexual selection [START_REF] Helft | First evidence of inclusive sexual selection in the ant Cataglyphis cursor : worker aggressions differentially affect male access to virgin queens (T Tregenza[END_REF]. By aggressing differentially males from different localities of origins, workers might therefore also play a role in limiting the mating efficiency of immigrant males. iv) Queens mate with males from different colonies to produce workers but use thelytokous parthenogenesis (originating daughters by asexual reproduction) to produce part, if not all, of the new queens [START_REF] Doums | Facultative use of thelytokous parthenogenesis for queen production in the polyandrous ant Cataglyphis cursor[END_REF][START_REF] Pearcy | Conditional use of sex and parthenogenesis for worker and queen production in ants[END_REF]). The remaining fraction of the new queens is produced by sexual reproduction. Queens also reproduce through arrhenotokous parthenogenesis (originating males by asexual reproduction) to produce new haploid males, as it is classically the case in Hymenoptera [START_REF] Trivers | Haploidploidy and the evolution of the social insect[END_REF]. This unusual reproductive system reduces male gene flow and limits genetic mixing in reproductive individuals, as most of the new queens and all males are asexually produced. Males contribute to gene flow through the small fraction of sexually produced queens (Figure 1A). Occasionally, males can also contribute to gene flow in orphaned colonies, as workers can asexually produce new males and queens in the absence of the mother queen (Figure 1A). These reproductive characteristics result in a strong genetic structure at fine spatial scale (from 100 m up to 120 Km; [START_REF] Clémencet | Hierarchical analysis of population genetic structure in the monogynous ant Cataglyphis cursor using microsatellite and mitochondrial DNA markers[END_REF][START_REF] Hardy | Small-scale spatial genetic structure in an ant species with sex-biased dispersal[END_REF], and may therefore favor evolution of reproductive isolation over small geographical scale. It is not yet clear whether C. cursor shares the same reproductive system and its population genetic structure is not yet known.

Here, we used an integrative approach to investigate the level of genetic and phenotypic differentiation between localities both within and between the two defined species C. piliscapa and C. cursor. We aimed to detect the highest possible within-species differentiation by sampling localities either separated by a geographic barrier (the Durance River for C. cursor) or at the extreme ranges of the geographic distribution (for C. piliscapa).

To this end, we compared the level of phenotypic and genetic differentiations within vs. between species by: i) considering three types of genetic markers (nuclear and mitochondrial sequences, and microsatellite markers); ii) measuring morphological traits (number of hairs, tibia length and morphology of male genitalia); and iii) performing a chemical analysis of the cuticular hydrocarbons. Finally, we further investigated mechanisms of reproductive isolation performing behavioral experiments to test for the mating success of conspecific and heterospecific males as well as the aggressive behavior of workers towards them.

MATERIALS & METHODS

Colony sampling

Between 2015 and 2017, workers from 369 distinct colonies covering 26 localities were sampled in Southern France and Northern Spain (Figure 1B; Table S1). For C. cursor, we sampled 15 localities on each side of the Durance River (213 workers), whereas for C. piliscapa, we sampled 11 localities at the two extreme ranges of its geographic distribution (156 workers). For all subsequent analyses, individuals were assigned to the species according to their geographic position, but we did not use any other a priori to assign individuals to a given species. Given the presence of a putative barrier to gene flow for C. cursor (the Durance River) and the sampling gap between the Northern and Southern distribution of C. piliscapa (Figure 1B), we refer for all the following analyses to four sampled regions, two for each species (Figure 1B). For a subset of these colonies (n = 241), one individual per colony was cooled down in an icebox immediately upon collection and frozen at the end of the day for later use for chemical analyses (Table S1). The remaining samples were stored in 95% ethanol (5%TE; Tris 10 mM: EDTA 1 mM) until DNA extraction. One worker per colony (n = 369) was genotyped at 13 microsatellite markers. Due to expected low within-locality polymorphism [START_REF] Clémencet | Hierarchical analysis of population genetic structure in the monogynous ant Cataglyphis cursor using microsatellite and mitochondrial DNA markers[END_REF], we only sequenced two workers per locality (n = 52 individuals) for the COI mitochondrial marker and 12 EPIC nuclear genes. We chose individuals from the most distant colonies in order to capture most of the total genetic variability [START_REF] Clémencet | Hierarchical analysis of population genetic structure in the monogynous ant Cataglyphis cursor using microsatellite and mitochondrial DNA markers[END_REF]. All genotyped workers were measured for the number of hairs on the first tergite of the gaster and the length of the tibia (n = 369). An additional dataset of 92 males was used to assess morphological differences in male genitalia. These males were obtained from 21 colonies sampled in seven localities (4 localities of the genetic sampling and 3 news ones) distributed within the four regions (Figure 1B, Table S1).

Genetic analysis

Total genomic DNA was extracted from half of the head using the extraction kit Qiagen (DNA Easy blood and tissue kit). One mitochondrial locus (Cytochrome Oxidase I, COI) was amplified following the protocol of [START_REF] Clémencet | Hierarchical analysis of population genetic structure in the monogynous ant Cataglyphis cursor using microsatellite and mitochondrial DNA markers[END_REF]. Twelve Exon-Primed-Intron-Crossing nuclear loci (EPIC) designed by [START_REF] Ströher | Exon-primed intron-crossing (EPIC) markers as a tool for ant phylogeography[END_REF] were sequenced (see loci characteristics and sequencing conditions in Table S2). Two individuals from C. italica (from Bari, Italia) were also sequenced and used as outgroup. PCR products of nuclear and mitochondrial loci were sent to GENOSCREEN (Lille, France) for Sanger sequencing. Base calling and sequence alignment were performed using CodonCode Aligner (CodonCode Corporation, Dedham, MA, USA). Phase reconstruction in nuclear loci is detailed in Table S2. All DNA sequences were deposited in GenBank database (accession numbers in Table S2). In addition, thirteen microsatellite loci were amplified for genetic analyses (see characteristics and PCR conditions in Table S3).

All analyses (genetic, morphological and chemical) were conducted on R v.3.2.2 (R Core Team 2016) except when otherwise stated.

EPIC loci. The Haploweb approach was used to delimit genetic clusters using phased haplotypes (alleles) of EPIC loci that co-occur within heterozygous individuals [START_REF] Flot | Haplowebs as a graphical tool for delimiting species: a revival of Doyle's 'field for recombination' approach and its application to the coral genus Pocillopora in Clipperton[END_REF]. Haploweb uses haplotype networks to identify connections between haplotypes co-occurring in heterozygous individuals to visualize and delineate single-locus 'fields for recombination' (FFRs). Haplotype networks were built with NETWORK v.10.2 using the median-joining method [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF]. The congruence of FFRs obtained from several independent loci is used to estimate clustering support. The congruence between nuclear markers is reported on a haploweb matrix. This matrix represents, for each pair of individuals, a score equal to the number of EPIC markers supporting the clustering of the two individuals (i.e., for which the two individuals belong to the same FFR). The low genetic diversity observed for the nuclear sequences did not allow inferring summary statistics and testing the correlation between genetic and geographical distances. mtDNA. Bayesian Inference (BI) was used to reconstruct phylogenetic relationships among haplotypes using the program MrBayes v.3.1.2 [START_REF] Ronquist | MrBayes 3: Bayesian phylogenetic inference under mixed models[END_REF]. The two samples of the sister species C. italica from Bari were used as outgroup. Tree topology and nodes support were estimated using a Metropolis-coupled Markov chain Monte Carlo sampling approach, with two runs implementing four chains (three hot, one cold) running for 1 x 10 7 generations, saving one tree every 100 generations. Convergence was checked by examining the standard deviation of the split frequencies between the two runs (< 0.01) and the Potential Scale Reduction Factor diagnostic. The first 25,000 trees of each run were discarded as burnin, and a majority-rule consensus tree was generated from the 75,000 remaining trees.

The genetic distance between localities was computed as the average proportion of nucleotide differences (p-distance). Visualization of the pairwise distance matrix was obtained by means of a non-parametric multidimensional scaling plot (NMDS) as implemented in the MASS library. Correlation between geographic and genetic distances was assessed by a Mantel test [START_REF] Mantel | The detection of disease clustering and a generalized regression approach[END_REF]) computed with the VEGAN library. We further tested the influence of the Rhône River on genetic differentiation using a partial Mantel test. To this end, we built a binary pairwise matrix with values of 0 corresponding to pairs of localities separated the river (belonging to different putative species) and 1 for pairs of localities on the same side of the river (belonging to the same species).

Microsatellite data.

We tested for linkage disequilibrium between loci and for Hardy-Weinberg equilibrium within each locality using exact tests within the Genepop Package [START_REF] Rousset | genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF]. Allelic richness was estimated using the rarefaction method (it standardizes to the smallest number of individuals sampled across all combinations of locus and locality multiplied by the ploidy level of the species) implemented in the package PopGenReport, allowing to control for different sample size [START_REF] El Mousadik | High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco[END_REF]. We used a mixed effect model with locus as random factor to test for a difference in allelic richness and gene diversity between the two species (considered as fixed factor in the model) using the package nlme. Pairwise GST [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF], which corrects for the FST maximal value given the observed within locality diversity, was estimated with the PopGenReport package. Similarly to mtDNA, GST distance matrix was visualized by an NMDS and the effects of the Rhône River and geographical distance on GST were tested using the same partial Mantel test. We run the Bayesian clustering approach implemented in STRUCTURE v.2.3 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF], varying K (number of clusters) from 1 to 26 (i.e., number of localities). We run the algorithm 10 times for each K with a combination of an admixture and a correlated-allele frequencies model. Each run included a 5 × 10 4 burn-in period followed by 1 × 10 5 iterations of the MCMC. Data convergence was assessed by analyzing Alpha equilibrium across iterations. The most likely number of genetic clusters was evaluated using the ΔK statistics [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF]) and the log-likelihood value generated by Structure Harvester v.0.6.8 [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF].

Chemical analyses of cuticular hydrocarbons

For each colony sampled in 2015 (the South region of C. piliscapa distribution was missing), chemical extractions were carried out a few weeks after sampling, with each worker individually retrieved from the freezer immediately before its chemical extraction. After removal of the gaster, the individual's body was immersed in 200 µl of pentane for 5 min to extract the cuticular hydrocarbons. Samples were analyzed by gas chromatography-mass spectrometry (GC-MS) following the protocol described in Table S4. Peak areas were integrated with MSD ChemStation software v.E.02.01.1177 (Agilent technologies). Chemical compounds were identified by their retention times and fragmentation patterns, and were confirmed through comparisons with previous results on C. piliscapa [START_REF] Monnin | Chemical characterization of young virgin queens and mated egg-laying queens in the ant Cataglyphis cursor: random forest classification analysis for multivariate datasets[END_REF]; Table S4).

To assess the chemical differentiation between the two species and between the two regions within C. cursor, we used a Random Forest analysis (RF) with the R package randomForest [START_REF] Liaw | Classification and regression by random forest[END_REF]. RF is a tree-based machine learning tool increasingly used in chemical studies, including the chemical ecology of social insects (i.e., [START_REF] Mitra | Nestmate discrimination in the social wasp Ropalidia marginata: chemical cues and chemosensory mechanism[END_REF][START_REF] Monnin | Chemical characterization of young virgin queens and mated egg-laying queens in the ant Cataglyphis cursor: random forest classification analysis for multivariate datasets[END_REF] as it deals well with multicollinearity, assumes no specific distribution of the variables and does not require transformation of constrained data (such as chemical compounds) [START_REF] Ranganathan | Reducing the babel in plant volatile communication: using the forest to see the trees[END_REF]. We used the relative abundances of peaks (i.e., the proportion each peak represented relative to the sum of all peaks) computed from the raw dataset, including all the zero values. We performed RF using 10,000 trees to determine the importance of each peak in classifying the individuals into the two species. We used the default values for the number of input variables randomly selected to build each node of the tree, and for the number of observations not used for building the current tree (i.e., the out of bag sample -OOB). This OOB sample (33% of individuals) was used to generate the confusion matrix and the OOB error rate (i.e., rate of unsuccessful assignment of the OOB samples to the correct species). A low OOB error rate indicates a high ability of the indicator variables to predict the species of origin of the individual. The pairwise proximity matrix between individuals estimated by the RF was visualized using a metric multidimensional scaling plot and was used to perform a Mantel test for testing the correlation between chemical and genetic distances (using GST). The individual chemical distances were averaged to obtain mean chemical distances between localities before performing the Mantel test.

Morphometric analyses of workers

The left foreleg of each genotyped worker (n = 369) was detached from the body, glued on a double-faced tape, and the tibia was measured to the nearest 0.01mm using a stereomicroscope (Leica Microsystems, Wetzlar, Germany) and the software ImageJ v.1.8 (NIH, Bethesda, MD, USA; http://imagej.nih.gov/ij/). In addition, we counted the number of hairs on the first tergite of the gaster. This trait is the main criterion used to discriminate C. cursor and C. piliscapa, the former having more than five hairs whereas the latter has less than five [START_REF] Agosti | A provisional list of the Balkan with a list to the worker caste. II. Key to the worker caste, including the European species without the Iberian[END_REF].

We first used the R library pdfCluster [START_REF] Azzalini | Clustering via Nonparametric Density Estimation: The R Package pdfCluster[END_REF] to cluster individuals based on the two morphological variables. This package performed clustering based on nonparametric density estimation. The first phase consists in determining space regions (cluster cores) containing objects with low dissimilarities and separated from each other by valley of high dissimilarities. The second phase consists in assigning the data points still unlabeled (those not falling into a determined cluster core) among the cluster cores. We used the default arguments proposed by the function pdfCluster. We performed the clustering based on both variables.

Second, we investigated the differences between the two species using a mixed effect model with the package nlme [START_REF] Pinheiro | nlme: Linear and Nonlinear Mixed Effects Models[END_REF], considering the species as a fixed effect whereas localities were considered as a random effect. We allowed heterogeneous variance between species for the number of hairs using the weight function (varIdent). We checked visually that residuals follow model assumptions for both dependent variables. The fixed effect was tested by comparing the model fitted with the maximum likelihood method with and without the fixed factor using a log-likelihood ratio test.

Geometric morphometrics of male genitalia

The 92 male genitalia were dissected under a binocular microscope after a 4-hour bath in KOH 10%, and only the left side of the genitalia was used for measurement. Four sclerified elements of the genitalia were separated (sternum IX, paramere, volsellae and penisvalvae; Boudinot 2013) and digitalized using a Leica M165C binocular microscope connected to a Leica DFC425 camera. Anatomical landmarks and sliding semi-landmarks were defined for each element using the software tpsDig2 [START_REF] Gunz | Semilandmarks: a method for quantifying curves and surfaces[END_REF][START_REF] Rohlf | The tps series of software[END_REF]. Five anatomical landmarks and 18 sliding semi-landmarks were placed on each sternum IX (paramere: 6 landmarks /44 semi-landmarks; volsella: 5 landmarks /42 semi-landmarks and penisvalvae: 4 landmarks /35 semi-landmarks). These landmarks were used to estimate the shape of each element for every individual. The sliding phase was performed by minimizing the bending energy between each specimen and the consensus [START_REF] Gunz | Semilandmarks: a method for quantifying curves and surfaces[END_REF]. The coordinates were then superimposed using a Generalized Procrustes Analysis [START_REF] Rohlf | Extensions of the procrustes method for the optimal superimposition of landmarks[END_REF] with the R library geomorph v.3.0.6 [START_REF] Adams | Geomorph: Software for geometric morphometric analyses[END_REF]). This analysis removes the effect of size, location and orientation. The centroid size (the square root of the sum of the squares of the distance from each point to the centroid of the set of landmarks), a measure of size, was also taken into account for the following analyses to control shape variation for size.

For each element of the male genitalia, variation in conformation between individuals was analyzed with a Principal Components Analysis (PCA). The first components accounting for at least 90% of the PCA total variation were used as new variables to define the shape of each element (shape variables) [START_REF] Baylac | Fourier descriptors, procrustes superimposition, and data dimensionality: an example of cranial shape analysis in modern human populations[END_REF]. This procedure reduces the dimensionality of the data set and works on independent variables [START_REF] Baylac | Fourier descriptors, procrustes superimposition, and data dimensionality: an example of cranial shape analysis in modern human populations[END_REF]. In order to reach 90% of the variation, five components were used as shape variables for sternum IX, eight for paramere, nine for volsellae and seven for penisvalvae. Multivariate analyses of variance (MANOVA) were performed to test for statistical difference in shape for each genitalia element. Three MANOVAs were performed. The first one was used to test for difference between species and the other two were used to test for difference between regions within each species. We used centroid size, the square root of the sum of the variances of the landmarks around the centroid in x and y directions to control for size in our analysis, as in the absence of allometry the centroid size is supposed to be uncorrelated with any shape variable [START_REF] Bookstein | Size and shape spaces for landmark data in two dimensions: comment[END_REF][START_REF] Jungers | Shape, relative size, and size-adjustments in morphometrics[END_REF]. Shape variations between species and regions were visualized using a canonical variate analysis (CVA) for each of the four elements.

Behavioral experiments

We used seven colonies from five localities that were collected either during the sampling described in this paper or in a previous sampling (Table S1). Colonies hibernated in the laboratory at about 14°C, and sexuals were produced once colonies were removed from hibernation. Males were removed from their colonies twice a week to prevent intra-nest mating and ensure that young queens remained virgins (gynes). Gynes were marked every three to four days to record their age. A total of 11 gynes (six C. piliscapa and five for C. cursor) and 22 males (11 for each species) were used (Table S5).

Each behavioral experiment was divided into three successive trials of 10 min separated by 15 min of resting time between two trials and was conducted with the same individuals. Each behavioral experiment took place in a mating box (10.5*9.5 cm) under natural light complemented with a heating lamp to raise temperature to ~35°C. This mimics mating conditions in the field, which occurs in spring during the hottest hours of the days. In Trial 1, we placed one gyne and two males (one per species) in the mating box, to test for unimpeded mate choice by the gyne. After having removed the males, we added 15 workers originating from the same colony of the gyne (i.e., nestmates) in the mating box and let them 15 min for acclimation. Trial 2 began when the two previous males were placed again in the mating box for 10 minutes, to test for the behaviors of workers towards males. Workers may differentially aggress the conspecific and heterospecific males and thus affect their mating success, as previously suggested within C. piliscapa [START_REF] Helft | First evidence of inclusive sexual selection in the ant Cataglyphis cursor : worker aggressions differentially affect male access to virgin queens (T Tregenza[END_REF]. After having removed the males, we let again the gyne and workers rest for 15 min. Trials 3 began when we added two foreign workers, each one being a sister of the males previously tested in Trials 1 and 2 (i.e., originating from the same colonies). The idea was to test whether workers behave more aggressively towards heterospecific than conspecific workers, and to test whether this effect was stronger when targeting males (Trial 2) or workers (Trial 3). In all trials, we paid attention to select males that are dark pigmented, active, moving around, and excited with their wings open to ensure their readiness to mate. Note that workers of C. piliscapa and C. cursor are not territorial and are generally peaceful with workers from colonies of the same locality, which can easily be adopted in natural colonies [START_REF] Nowbahari | Age-related changes in aggression in ant Cataglyphis cursor (Hymenoptera, Formicidae): influence on intercolonial relationships[END_REF] and are found at low frequency in field colonies [START_REF] Doums | Worker thelytoky allows requeening of orphaned colonies but increases susceptibility to reproductive cheating in an ant[END_REF]. We performed 11 behavioral experiments, five with a C. cursor gyne (all originating from a single colony) and six with a C. piliscapa gyne (five gynes from one colony and one gyne from another one) (see Table S5). The conspecific foreign male and workers came from a different locality than the gyne. They were from the same region in five experiments (with C. piliscapa gynes) and from a different region in six experiments (five experiments with C. cursor gynes and one experiment with C. piliscapa gyne) (Table S5). Unfortunately, the low number of sexuals produced did not allow us to robustly test for the occurrence of indirect sexual selection within each species, calling for a cautious interpretation of the results.

Each trial of 10 min was videotaped with a Sony Handycam HD and the behaviors were recorded using the free software BORIS v.7.4.10 (Friard and Gamba 2016). Before introduction into the mating box, the two males and two foreign workers were marked with a color dot on the thorax using UniPaint Marker. In Trials 1 and 2, we recorded the number of mating attempts (as described in [START_REF] Helft | First evidence of inclusive sexual selection in the ant Cataglyphis cursor : worker aggressions differentially affect male access to virgin queens (T Tregenza[END_REF] and the number and duration of successful matings (male and gyne locked into mating). In Trials 2 and 3, we recorded the duration of biting of workers towards foreign males or workers. Workers sometimes hold biting for a prolonged period of time so that the number of biting alone does not appropriately reflect the aggression intensity.

Using data from Trials 1 and 2, in which males were present, we tested whether the number of mating attempts differed as a function of the origin of males (conspecific versus heterospecific), of the presence/absence of workers, and of the gyne species (C. cursor or C. pilliscapa) using a mixed effect model with the behavioral experiment as a random factor. Each main fixed effect was tested by comparing the full model with a reduced model without the tested main effect using a log likelihood ratio test. Homoscedasticity and normality of residuals were checked visually. Using data from Trials 2 and 3, in which workers that are sisters of the gyne were present, we tested whether the duration of aggressive behaviors differed as a function of the origin of individuals (conspecific versus heterospecific), of the type of foreign individuals (males or workers) and the gyne species (C. cursor or C. pilliscapa) using the same model than above. In this model, we included the interaction between the type of individuals introduced (males or workers) and their origin, as the behavior of worker sisters of the gyne towards introduced individuals is expected to differ according to their type. When the interaction was significant, mean comparisons were performed using the package emmeans.

RESULTS

Genetic analysis

Nuclear sequences. Nine out of the twelve nuclear EPIC markers were successfully sequenced for all individuals analyzed (n = 52, Table S2). After reconstructing the phase for the nine markers, we found an average of 8.8 haplotypes and 9.4 SNPs per marker (Table S2). We found a single SNP (EPIC 1503) discriminating the two species, (i.e., without shared alleles between the two). Haploweb uncovered three main genetic clusters (Figure 2): one cluster grouping all C. cursor localities, and the two others grouping all C. piliscapa localities (the haploweb for each EPIC marker are given in Figure S1). The outgroup C. italica was well differentiated. The two distinct clusters within C. piliscapa (dark and light blue color, Figure 2B) match the geographical regions except for the localities of VIL and FRO. However, although the Northern localities of C. piliscapa form a homogenous cluster, the result is less evident for the South region of this species. Within C. cursor (dark and light red color, Figure 2B), two clusters could also be distinguished (even though less clearly) and differentiate the two sides of the Durance River, except for the localities CAD and ROV.

MtDNA: the mitochondrial COI marker was successfully sequenced for all individuals analyzed (n = 52). We found 23 unique haplotypes defined by 93 SNPs, representing 16% of the total sequence. In 21 out of the 26 localities, the two individuals sequenced shared the same haplotype. Five SNPs were fixed for a different allele in each species, potentially representing markers for species identification. The phylogenetic tree (Figure 2A) identifies two main clusters corresponding to the two species. The mean p-distance between regions of C. piliscapa was 0.061 and of 0.029 between regions of C. cursor, while the p-distance between species was 0.064. This was reflected in the pattern shown by the NMDS plot (Figure 3A), as it separated only three main clusters that were more or less genetically equidistant. This analysis identified the two clusters previously uncovered within C. piliscapa (except for localities FRO and VIL). Within C. cursor, the two sides of the Durance River were only moderately discriminated, with discrepancies for the localities CAD and ROV. We found a significant positive association between the mean p-distance and geographical distances after controlling for the presence of the Rhône River with the partial Mantel test (r = 0.40; P = 0.001, Figure 4A). Similarly, the presence of the Rhône River separating pairs of localities also significantly increased genetic differentiation after controlling for the geographical distances (partial mantel test, r = 0.61; P = 0.001), showing that for a given geographical distance, localities from the two species were more genetically differentiated than within species. This was always true for C. cursor but not for C. piliscapa (Figure 4A), where some pairs of localities showed the same level of genetic distance than pairs of localities from the two species (the matrix of pairwise genetic distances is provided in Table S6). These pairs of localities were those including one locality from each of the two genetic clusters identified by the NMDS plot within C. piliscapa (i.e., separating FRO, VIL, EMP, TET, SAC and TRO from the rest of the localities of this species).

Microsatellite data:

We did not detect linkage disequilibrium between microsatellite loci (Table S7) and all except one HW equilibrium tests were not significant after Bonferroni's correction (Table S8). All thirteen microsatellite loci used in this study were polymorphic with the allelic richness per locality and locus ranging from 1 to 7.37 (mean± SD = 3.65 ± 1.49) and gene diversity from 0 to 0.95 (mean± SD = 0.61 ± 0.24) (Table S9 andS10). Interestingly, C. piliscapa presented a significantly higher allelic richness (4.22 ± 0.93) than C. cursor (3.23 ± 0.32) for all but one locus (mixed effect model, L-ratio = 55.23, P < 0.0001, Table S9). Gene diversity showed a similar pattern (C. piliscapa: 0.68 ± 0.12; C. cursor: 0.56 ± 0.06; mixed effect model: L-ratio = 28.18, P < 0.0001, Table S10).

Highly significant population differentiation was found among all but three comparisons (these three non-significant comparisons were found between neighboring localities of C. cursor in the South region, Table S11). The global Hedrick's GST was high (0.73 ± 0.20) (pairwise values of GST and FST are given Table S11). The genetic differentiation was always high between species (mean±sd, GST = 0.81±0.07) as well as between regions for both species (mean±sd C. cursor: GST = 0.81±0.09; C. piliscapa GST = 0.80±0.08) whereas it was more variable, with a lower mean when comparing localities within region for both species (mean±sd C. cursor: GST = 0.43±0.24; C. piliscapa GST = 0.59±0.27) (Table S11). Localities from the same geographic regions clustered together (with two exceptions, FRO and EMP, see Figure 3B) but the four clusters were roughly equidistant. The microsatellites better discriminated the two regions within C. cursor than the mtDNA (Figure 3). The proportion of private alleles between species was 39.5% (109 out of 276 alleles) with 48% of them found in a single locality and all except two occurring at low frequency (less than 13%). Interestingly, these two private alleles both found in C. piliscapa at high frequency (45% for one and 27% for the other) also discriminated the two regions of C. piliscapa, with the exception of VIL. For both species, the proportion of private alleles between regions was higher than between the two species (within C. piliscapa: 50% (131/263) and within C. cursor: 50% (90/180)).

In agreement with the large level of pairwise genetic differentiation between localities, STRUCTURE analyses detected a large number of genetic clusters with the optimal k being 10; K = 4 being the second best solution (Figure S2). At K = 10, the clusters identified by STRUCTURE were well delimited and consistent among runs. In all individuals, except one of C. cursor (from LAC), the highest level of assignment matched their region of origin and therefore species (Figure S2). Only one locality of C. cursor (CAD), harbored a high number of individuals (12/15) where the largest assignment value was below 80%. Interestingly, this locality also showed discrepancy for the EPIC and mtDNA markers suggesting ongoing gene flow across the Durance River. Only two percent of individuals (7/369; five in C. piliscapa and two in C. cursor) had more than 20% assignment corresponding to populations of the other species. When forcing to K = 4, the four clusters defined by STRUCTURE were consistent among runs and corresponded to the four geographical regions (except for ROV). When forcing to K = 2, the results were not consistent among runs, with no clear biological interpretation of the results (Figure S2). Overall, these results reflect the strong genetic differentiation between regions within species and among localities within a given region with low level of admixture (Figure S2).

The partial Mantel test showed a significant positive association between pairwise GST and geographical distances after controlling for the presence of the Rhône River (r = 0.32; P = 0.001). Similarly, including the presence of the Rhône River separating pairs of localities significantly increased genetic differentiation after controlling for the geographical distances (partial Mantel test, r = 0.18; P = 0.006). This is graphically represented in Figure 4B where inter-specific comparisons always have high values of genetic differentiation (mean±sd GST = 0.79±0.06) whereas intra-specific comparisons for the same geographical distances (between 20 and 50 km) vary from low to high values of genetic differentiation within both species (mean±sd C. cursor: GST = 0.65±0.25; C. piliscapa GST = 0.53±0.25).

Chemical analyses

Our chemical analyses identified 44 peaks containing 74 compounds (Table S4). When training the RF analysis for discriminating the two species, the rate of individual misclassification to their group (OOB error rate) was very low (0.8%), showing that the species could confidently be identified from their cuticular hydrocarbons. Only two misclassifications (out of 239 individuals) were observed (Figure 5), from the localities of StM and LUN. Training the RF to differentiate the two regions within C. cursor led to an OOB error rate of 5% (8 misclassifications out of 145 individuals; Table S12). The chemical analysis further segregated two clusters within C. cursor, which roughly correspond to the North and South regions identified using genetic markers. Note that our chemical analyses did not include individuals from the South region of C. piliscapa (see Materials and Methods). We observed a significant correlation between chemical and genetic distances (partial Mantel test, r = 0.42; P = 0.005) only due to the similar high genetic (based on GST) and chemical differentiation observed between localities of the two species (Figure S3).

Morphometric analyses of workers

The clustering analysis revealed four different clusters based on both tibia size and the number of hairs (Figure 6). Ninety-eight percent (147/150) of the individuals from clusters 1 and 2 belonged to C. piliscapa whereas 96% (210/219) of individuals from clusters 3 and 4 belonged to C. cursor (overall, twelve individuals were misclassified, Table S12). The clustering analysis confirmed that the two species can be discriminated with relatively high accuracy based on the number of hairs, with the value of five representing the threshold between the two species. Considering only the number of hairs, 50% of individuals with five hairs (8/16) were wrongly assigned. Not surprisingly, in C. cursor, the three wrongly assigned individuals had five hairs and in C. piliscapa, the nine wrongly assigned individuals had either five, six or nine hairs. However, by averaging the number of hairs by locality, all localities were properly assigned (see Table S1 for the mean number of hairs per locality) suggesting that sampling many individuals from a given colony would allow a proper determination at the colony level. Within C. piliscapa, the two other clusters (1 and 2) determined by the analysis were discriminated by the tibia size with cluster 2 containing larger workers. However, they did not correspond to the two geographical regions since 52% and 63% of individuals from respectively clusters 1 and 2 belonged to the South region of C. piliscapa. Within C. cursor, the two clusters (3 and 4) were differentiated mainly by the number of hairs with cluster 4 having more hairs (Figure 6). Interestingly, 82% of the individuals from cluster 4 belonged to the North region of C. cursor, whereas 74% of the individuals from cluster 3 belonged to the South region.

In agreement with the clustering analysis, the number of hairs significantly differed between the two species (L-ratio = 49.09, P < 0.0001) with an average of 0.97 (± 1.77) hairs for C. piliscapa and of 12.95 (± 5.16) for C. cursor. Although the differences are small, the two species also differed significantly by their tibia size (L-ratio = 9.27, P = 0.002), with C. piliscapa having longer legs (2.14 ± 0.20) than C. cursor (2.06 ± 0.17). Within each species, the number of hairs significantly differed between the two regions (C. piliscapa: L-ratio = 10.12, P = 0.002; C. cursor: L-ratio = 25.77, P < 0.0001) whereas the tibia size did not significantly differ (C. piliscapa: L-ratio = 0.21, P = 0.64; C. cursor: L-ratio = 0.02, P = 0.96).

Morphometric analyses of male genitalia

The shape of the four elements of the genitalia significantly differed between the two species (MANOVA: P < 0.001 for all the four elements; Figure 7). This was especially clear for paramere and volsellae as the two species scattered along the two axes of the CVA (Figure 7). The North and South regions of C. cursor differed in shape for volsellae (P < 0.001) and penisvalvae (P < 0.001), but not for paramere (P = 0.51) and sternum IX (P = 0.529). The North and South regions of C. piliscapa differed in shape for all the four elements (Sternum IX: P = 0.003; paramere: P < 0.001; volsellae: P < 0.001; penisvalvae: P < 0.001).

Behavioral experiments

The number of mating attempts was about twice higher with conspecific males than with heterospecific ones (L-ratio = 28.09, df = 1, P < 0.001), and in the absence of workers than in their presence (L-ratio = 6.77, df = 1, P = 0.009) (Figure 8A,8B). The interaction between the type of males and the presence of workers was not significant (L-ratio = 0.87, df = 1, P = 0.35). Overall, the number of mating attempts did not significantly differ between the gynes of the two species (L-ratio = 0.94, df = 1, P = 0.33). Interestingly, C. cursor and C. piliscapa males behaved differently towards heterospecific gynes. None of the six C. cursor males tried to mate with a C. piliscapa gyne, whereas three out of five C. piliscapa males tried to mate with a C. cursor gyne (Figure 8A, 8B, Table S5). As expected, males of both species attempted to mate with conspecific gynes (at least in the absence of workers). All except one C. piliscapa gynes mated with the conspecific male whereas no C. cursor gyne mated whatever the males (Figure 8A, 8B, Table S5), even in absence of workers. Although this can result from the C. cursor gynes being not yet sexually mature (even though males tried to mate, which usually indicates that gynes were sexually mature too; C. Doums & T. Monnin, pers. obs.), it can also suggest that the conspecific males were too genetically dissimilar (they all originated from a different region whereas males were from the same region for C. piliscapa gynes).

Workers were more aggressive towards heterospecific individuals (males or workers) than towards conspecifics. However, the extent of the difference varied according to the type of individuals (males or workers) (significant interaction terms: L-ratio = 4.27, df = 1, P = 0.04, Figure 8C,8D). Even though both heterospecific males and heterospecific workers were significantly more aggressed than their conspecific counterparts (males: contrast estimate = 27.3, t = 3.26, df = 36.3, P = 0.0024; workers: contrast estimate = 51.4, t = 6.13, df = 36.3, P < 0.001), the extent of the difference was twice higher for workers (Figure 8D) than males (Figure 8C). The reason for this is that while conspecific foreign workers (Figure 8D) were rarely aggressed, conspecific males were often victims of attacks (Figure 8C), as previously observed [START_REF] Helft | First evidence of inclusive sexual selection in the ant Cataglyphis cursor : worker aggressions differentially affect male access to virgin queens (T Tregenza[END_REF]. The only exception to this rule concerned one outcome where the conspecific workers suffered the highest level of aggression (Figure 8D). Interestingly, this was the only test using a C. piliscapa gyne and a conspecific male originating from a different region and resident workers were very aggressive towards the two foreign workers, suggesting high aggression between regions in C. piliscapa. The workers of the two species did not significantly differ in their duration of aggression (L-ratio = 1.00, df = 1, P = 0.32).

DISCUSSION

We combined genetic, morphological, chemical and behavioral analyses to assess degree of differentiation between the sister ant species, C. cursor and C. piliscapa [START_REF] Bondroit | Les fourmis de France et de Belgique[END_REF][START_REF] Agosti | A provisional list of the Balkan with a list to the worker caste. II. Key to the worker caste, including the European species without the Iberian[END_REF]) could be well differentiated. These two species were supposed to have a non-overlapping distribution on each side of the Rhône River [START_REF] Nowbahari | Individual, geographical and experimental variation of cuticular hydrocarbons of the ant Cataglyphis cursor (Hymenoptera: Formicidae): Their use in nest and subspecies recognition[END_REF]). We therefore specifically compared localities on each side of the river, assigning individuals to species according to their geographical origin in order to test independently each marker (even the taxonomic character). Our sampling was designed to maximize within-species differentiation by sampling localities either at the extreme ranges of the distribution (C. piliscapa) or on each side of a putative geographic barrier (the Durance River for C. cursor). Our results revealed that the between-species differentiation is only slightly higher than withinspecies differentiation, especially because of the high differentiation found between regions within species. Three markers, the mitochondrial DNA, the cuticular hydrocarbons and the number of hairs, were able independently distinguish the two species with no or low error rates whereas they were less successful in distinguishing the two regions within species (Table S12). Interestingly, we observed higher levels of aggression of workers towards heterospecific males that could lead to reproductive isolation between species by limiting mating attempts. No mating of heterospecific males was observed between species, even in absence of workers in our laboratory setting. Even though this result should be confirmed on a larger sample size, they suggest reproductive isolation between the two species. Altogether, our results support the taxonomic separation of the two species even though the genetic integrity of the two species is still far from clear.

The three genetic markers do not offer a similar pattern even though all of them showed a slightly higher genetic differentiation between species than within species. MtDNA phylogenetic tree separated the two species; this marker therefore appears to be a promising marker for species identification, with five SNPs being species-diagnostic (Table S12), though additional sampling is needed to confirm that these SNPs are discriminant over all the distribution area. However, the genetic differentiation between the two species did not reach the level observed between them and the most closely related species C. italica (Figure 2). Moreover, the level of genetic differentiation was only slightly higher between species (pdistance = 0.064) than between regions, at least for C. piliscapa (p-distance = 0.061). Among the nine nuclear EPIC markers studied, only a single SNP was found to discriminate the two species. These markers revealed three main clusters (Figure 2), identifying two regions of C. piliscapa and a single cluster for C. cursor. Most EPIC markers shared at least one haplotype between the two species, which could reflect a lack of complete reproductive isolation or may stem from incomplete lineage sorting, whereby the low substitution rate of these genetic markers did not allow to accumulate enough divergence to fully segregate the two species. The microsatellite markers showed a more complex picture, where higher between-species genetic differentiation was only observed over short geographical distances (between 20 and 50 km, Figure 4). In both species, the genetic differentiation was high between localities, even within regions, and could be as high as between species. Accordingly, ten well delimited groups were detected by STRUCTURE and the two species could not be consistently distinguished for K = 2 (Figure S2). For K = 10, few individuals in both species (2%) had more than 20% assignment corresponding to populations of the other species. Overall, the genetic markers therefore suggest that speciation is still ongoing between the two species and that the various levels of population genetic differentiation, even including those observed within each species, could represent pairs of populations at varying stages of the speciation continuum [START_REF] Feder | The genomics of speciation-with-gene-flow[END_REF][START_REF] Shaw | Speciation Continuum[END_REF]Stankowski andRavinet, 2021a, 2021b). However, our genetic markers were not enough powerful to confirm hybridization, and to disentangle current gene flow from shared ancestral polymorphism, as was done for instance by [START_REF] Cordonnier | From hybridization to introgression between two closely related sympatric ant species[END_REF] in the ant Tetramorium. Hybridization is a common phenomenon in ants that could promote adaptive genetic changes [START_REF] Kulmuni | Introgression in hybrid ants is favored in females but selected against in males[END_REF][START_REF] Martin-Roy | Differences in thermal tolerance between parental species could fuel thermal adaptation in hybrid wood ants[END_REF] and is often associated with major changes in social organization and reproductive systems (i.e. [START_REF] Schwander | Characterization and distribution of Pogonomyrmex harvester ant lineages with genetic caste determination[END_REF][START_REF] Eyer | Combined hybridization and mitochondrial capture shape complex phylogeographic patterns in hybridogenetic Cataglyphis desert ants[END_REF][START_REF] Stolle | Recurring adaptive introgression of a supergene variant that determines social organization[END_REF]. Whole genome analyses would be very promising in our study system to investigate these questions. They would also allow to identify putative genomic regions directly involved in reproductive isolation.

Morphological (number of hairs, size of workers and male genitalia) and chemical (cuticular hydrocarbons) characters also showed a higher between-species than withinspecies differentiation. The taxonomic morphological criterion (number of hairs) can easily be used to identify the two species with low error rate (3%). Interestingly, our quantitative assessment of the threshold value (5 hairs) was the same than the one stated in the determination key by [START_REF] Agosti | A provisional list of the Balkan with a list to the worker caste. II. Key to the worker caste, including the European species without the Iberian[END_REF]. The chemical assessment was also successful in identifying the two species with a very low error rate at the individual level (0.8%). Interestingly, the misclassified individuals are not the same for the different criteria used, suggesting that the misclassification errors likely stem from variations within each criterion or are old generation backcrosses. Cuticular hydrocarbon (CHC) blends are used by social insects as cues to discriminate nestmates from non-nestmates. They are influenced by both genetic and environmental factors [START_REF] Walsh | Ant cuticular hydrocarbons are heritable and associated with variation in colony productivity[END_REF] and are usually variable between and within species [START_REF] Bos | Recognition of social identity in ants[END_REF][START_REF] Sturgis | Nestmate recognition in ants (Hymenoptera: Formicidae): a review[END_REF][START_REF] Eyer | Divide and conquer: Multicolonial structure, nestmate recognition, and antagonistic behaviors in dense populations of the invasive ant Brachymyrmex patagonicus[END_REF][START_REF] Blumenfeld | Consistent signatures of urban adaptation in a native, urban invader ant Tapinoma sessile[END_REF]. Our results agree with other studies showing that chemical profiles can be used as relevant taxonomic tools in ants [START_REF] Blomquist | Insect Hydrocarbons: Biology, Biochemistry, and Chemical Ecology[END_REF][START_REF] Seppä | The many ways to delimit species: Hairs, genes and surface chemistry[END_REF][START_REF] Kather | Cuticular hydrocarbon profiles as a taxonomic tool: advantages, limitations and technical aspects[END_REF][START_REF] Villalta | Origin and distribution of desert ants across the Gibraltar Straits[END_REF]. The morphology of male genitalia is also frequently used as a taxonomic tool due to its high rate of evolution [START_REF] Hosken | Sexual selection and genital evolution[END_REF][START_REF] Simmons | Sexual selection and genital evolution[END_REF][START_REF] Kjer | Progress, pitfalls and parallel universes: a history of insect phylogenetics[END_REF]. However, although significant differentiation was detected between species (especially for the paramere), species distinction only based on this morphological character might be arduous unless using clustering analysis. The absence of major morphological differentiation in male genitalia supports the idea that reproductive isolation during the speciation process could have been mainly driven by pre-copulatory mechanisms of sexual selection and/or mate recognition based on chemical signatures. This hypothesis deserves further tests in subsequent studies on this system.

Our behavioral experiments, even though based on a limited sample size, suggest that the chemical differences detected between the two species are perceived, potentially with other morphological cues, and used by workers in the context of mating. Workers are usually not aggressive towards workers from other colonies, and foreign workers can be easily adopted if they do not come from geographically distant localities [START_REF] Nowbahari | Age-related changes in aggression in ant Cataglyphis cursor (Hymenoptera, Formicidae): influence on intercolonial relationships[END_REF]. On the contrary, workers belonging to the two distinct species are rapidly attacked and killed (Nowbahari et al. 1990, Figure 8, table S5). Our behavioral experiments using foreign workers as focal individuals confirm these results as they show stronger aggression towards heterospecific workers than conspecific ones. Similarly, heterospecific males were significantly more aggressed than conspecific males. Interestingly, the level of aggression towards conspecific males was higher than towards conspecific workers. This suggests that aggression towards males is specifically linked with the mating process rather than simply associated with worker colonial recognition. Our laboratory setting mimicking natural conditions ensures that these results reflect behavior occurring in natural populations. Moreover, similar aggressive behaviors of workers towards conspecific males from the same locality were already observed in a previous lab study [START_REF] Helft | First evidence of inclusive sexual selection in the ant Cataglyphis cursor : worker aggressions differentially affect male access to virgin queens (T Tregenza[END_REF] and on the field [START_REF] Cronin | Opportunities for mate choice in the fission-performing ant Cataglyphis cursor[END_REF]C. Doums & T. Monnin, pers. obs.). If workers preferentially aggress migrant males (i.e., heterospecific or from a distant locality), this indirect sexual selection could further limit male's contribution to gene flow and play a major role in population genetic differentiation and speciation. Sexual selection is known to be a major cause of reproductive isolation between species [START_REF] Coyne | Speciation[END_REF]). To our knowledge, our results are the first line of evidence that this could be the case in a social insect species. Estimating the strength of pre-mating reproductive isolation and sexual selection [START_REF] Rolán-Alvarez | Estimating sexual selection and sexual isolation effects from mating frequencies[END_REF][START_REF] Sobel | Unification of methods for estimating the strength of reproductive isolation[END_REF] for various pairs of populations differing by their level of genetic and phenotypic differentiation would shed light on the evolution of reproductive isolation at the different stages of speciation.

Our study also allows a better assessment of the population genetic structure of each species. Although this had been previously studied in C. piliscapa, this was only in the southern region [START_REF] Clémencet | Hierarchical analysis of population genetic structure in the monogynous ant Cataglyphis cursor using microsatellite and mitochondrial DNA markers[END_REF], and nothing was known about C. cursor. The strong genetic structure observed in both species was expected given the mode of colony foundation by fission that reduces female gene flow by limiting colony dispersal to ant walking distances (few meters) and by diminishing the number of new colonies founded due to the high energetic cost in both putative species [START_REF] Cronin | Recurrent evolution of dependent colony foundation across eusocial insects[END_REF]. As observed in our study, this mechanism generally yields a pattern of isolation-by-distance among localities and high genetic differentiation between localities, consistent to what has been previously described in many ant species [START_REF] Seppä | Gene flow and population viscosity in Myrmica ants[END_REF][START_REF] Liautard | Restricted effective queen dispersal at a microgeographic scale in polygynous populations of the ant Formica exsecta[END_REF][START_REF] Doums | Population genetic structure and male-biased dispersal in the queenless ant Diacamma cyaneiventre[END_REF][START_REF] Knaden | Fundamental difference in life history traits of two species of Cataglyphis ants[END_REF][START_REF] Leppänen | Comparative phylogeography of the ants Myrmica ruginodis and Myrmica rubra (A Parmakelis, Ed.)[END_REF]. Even though flying males could erode the isolationby-distance pattern induced by low female dispersal, they are not very good flyers [START_REF] Hardy | Small-scale spatial genetic structure in an ant species with sex-biased dispersal[END_REF]. Moreover, their contribution to gene flow can be ineffective if thelytoky is used at high frequency to produce new queens [START_REF] Pearcy | Conditional use of sex and parthenogenesis for worker and queen production in ants[END_REF][START_REF] Doums | Facultative use of thelytokous parthenogenesis for queen production in the polyandrous ant Cataglyphis cursor[END_REF]. The magnitude of genetic differentiation among populations within a species can be an indicator of the potential for speciation, as shown in vertebrates [START_REF] Bernatchez | Comparative phylogeography of Nearctic and Palearctic fishes[END_REF] and insects [START_REF] Ikeda | Erratum: Corrigendum: Loss of flight promotes beetle diversification[END_REF][START_REF] Polato | Narrow thermal tolerance and low dispersal drive higher speciation in tropical mountains[END_REF]. For instance, marine species with non-planktotrophic development and therefore with low dispersal show higher speciation rates than do planktotrophic forms [START_REF] Hansen | Modes of larval development and rates of speciation in early tertiary neogastropods[END_REF]. Whether colony fission and conditional parthenogenesis increase speciation rates in ants would merit further investigation.

The level of genetic diversity was high in both species even if slightly lower in C. cursor (Table S9 andS10). Whether this reflects different reproductive strategies (use or not of thelytoky) awaits further comparative analysis. The high level of genetic diversity and allelic richness of the microsatellite markers could somehow be surprising under a high level of thelytoky [START_REF] Pearcy | Conditional use of sex and parthenogenesis for worker and queen production in ants[END_REF][START_REF] Doums | Facultative use of thelytokous parthenogenesis for queen production in the polyandrous ant Cataglyphis cursor[END_REF]) and low male dispersal [START_REF] Hardy | Small-scale spatial genetic structure in an ant species with sex-biased dispersal[END_REF]. Indeed, both mechanisms should increase the probability of mating between related males and females and subsequently lead to high level of inbreeding in workers. This discrepancy suggests that thelytoky might be lower than estimated directly based on gyne production in C. piliscapa. This could be due to a higher fitness of sexually produced gynes or to strong variation in the use of thelytoky over space and time not detected by our previous estimates [START_REF] Pearcy | Conditional use of sex and parthenogenesis for worker and queen production in ants[END_REF][START_REF] Doums | Facultative use of thelytokous parthenogenesis for queen production in the polyandrous ant Cataglyphis cursor[END_REF]. Moreover, indirect sexual selection by workers could also play a role in preventing inbreeding and/or outbreeding depression by being more aggressive towards highly related/differentiated males. Experiments are currently being conducted in the laboratory to test these hypotheses, assessing the spatial and temporal variation in the rate of thelytoky in the two species and testing for a cost of inbreeding on queen fitness. The C. cursor species complex is therefore characterized by strong genetic differentiation not only between species but also among populations within species. The differentiation observed between species was only slightly higher than within species. However, it is not clear whether these patterns reflect genetic structure related with divergence time or isolation by distance without necessarily involving stronger reproductive isolation. By using an interdisciplinary approach, we found that the genetic differentiation between species is concomitant to morphological and chemical differentiation, as well as to an absence of mating between species and an elevated aggressiveness of workers towards heterospecific males. Behavioral experiments would be necessary within each species to test whether the strong genetic differentiation found between some pairs of populations, especially between regions, is associated with reproductive isolation. Our study therefore emphasizes that assessing phenotypic traits, antagonism, and mating attempt/success between individuals from distinct putative species is important in species with low dispersal and strong population structure to position the diverging populations along the speciation continuum. Species complex with strong hierarchical genetic structure can provide ideal case studies of taxon pairs near the completion of speciation, which could ultimately help to understand the evolution of reproductive isolation along a speciation continuum [START_REF] Stankowski | Defining the speciation continuum[END_REF]. S1. S1). Localities are aligned in the same order as the phylogenetic tree. This heatmap matrix represents, for each pair of individuals, a score equal to the number of independent EPIC markers supporting the hypothesis that the two compared individuals belong to the same FFR ("fields for recombination", i.e., the same haplotype is co-occurring in heterozygous individuals) indicating by level of grey. Black reflects the highest scores (9 out of 9; the pair of individuals belong to the same FFR for all 9 EPIC markers) whereas the lowest scores (0 out of 9) are shown in white.

Figure 1 :

 1 Figure 1: (A) Schematic representation of the mating system observed in C. piliscapa. Individuals are colored according to their genetic background. Faint arrows represent the small contribution of males to gene flow. (B) Sampling localities of C. cursor and C. piliscapa. Each locality is colored according to its species and region of origin (red for C. cursor South region, pink for C. cursor North region, dark blue for C. piliscapa South region and light blue for C. piliscapa North region). Samples used for genetic, morphological and chemical analysis are coded by a number in a circle whereas the four localities used for genitalia measurements and/or behavioral experiments are coded by a letter in a square. Detailed sample size for each analysis and each locality can be found in TableS1.

Figure 2 :

 2 Figure 2: (A) Bayesian tree based on the mitochondrial CO1 marker. Color boxes at the end of the tip correspond to the color of the geographical region of Figure 1B. The two black boxes correspond to the two samples of the outgroup species C. italica; (B) Conspecific heatmap matrix based on the haploweb of the nine EPIC nuclear markers (FigureS1). Localities are aligned in the same order as the phylogenetic tree. This heatmap matrix represents, for each pair of individuals, a score equal to the number of independent EPIC markers supporting the hypothesis that the two compared individuals belong to the same FFR ("fields for recombination", i.e., the same haplotype is co-occurring in heterozygous individuals) indicating by level of grey. Black reflects the highest scores (9 out of 9; the pair of individuals belong to the same FFR for all 9 EPIC markers) whereas the lowest scores (0 out of 9) are shown in white.

Figure 3 :

 3 Figure 3: Non-parametric multidimensional scaling plot over two dimensions based on (A) the haplotypic distance (p-distance) for the mitochondrial DNA and (B) the Hedricks'Gst for the microsatellites data.Each point is a locality. The species is identified by the shape (triangle for C. piliscapa and square for C. cursor) and the region of origin by the same color than in Figure1B. The locality identity is given for those that do not cluster with the other localities of their region. The stress value of the NMDS was 3.67 for the mtDNA and of 20.92 for the microsatellite data.

Figure 4 :

 4 Figure 4: Patterns of isolation-by-distance for (A) the mitochondrial marker based on the mean pdistance and (B) the microsatellite markers based on the Hedricks'Gst. Both are expressed on the logarithm of geographic distance. Pairs of localities of the same species are colored in red for C. cursor and blue for C. piliscapa whereas pairs comparing localities of the two species are black.

Figure 5 :

 5 Figure 5: Multidimensional scaling plot based on the cuticular hydrocarbon profile of workers. The metric multidimensional scaling plot was based on the pairwise proximity matrix between individuals estimated by the RF. When the classification of the RF is congruent with the origin of the individuals, the triangles (individuals classified as C. piliscapa) should be light blue (individuals originating from C. piliscapa North region), the circles (individuals classified as C. cursor) should be red (individuals originating from C. cursor North region) or pink (individuals originating from C. cursor South region). Note here that only two individuals noted with an arrow were misclassified into their species of origin by the RF.

Figure 6 :

 6 Figure 6: Clustering based on the number of hairs (X axis) and tibia size (Y axis) of workers. Each point corresponds to an individual identified by the clusters it belongs to according to the clustering analysis (labelled from 1 to 4) and by the species (red or pink for C. piliscapa from the South and North region respectively and blue or light blue for C. cursor from the South and North region respectively). The clusters 1 and 2 correspond to C. piliscapa whereas clusters 3 and 4 correspond to C. cursor. The two species can clearly be discriminated based on the number of hairs. Misclassification can be observed for 12 individuals (cluster 1 or 2 classified as C. cursor (blue color) and cluster 3 classified as C. piliscapa (pink color).

Figure 7 :

 7 Figure 7: First two axes of the canonical variate analysis for the four genitalia elements: A) Volsellae B) Paramere C) Sternum IX and D) Penisvalvae. Ellipses represent the 95% confidence intervals of the conformation average for each group.

Figure 8 :

 8 Figure 8:For each of the 11 tests, the number of mating attempts by the conspecific and heterospecific males competing for the tested gyne is given in absence of workers (A) and in their presence (B), as well as the percentage of time the introduced individuals (males (C) or foreign workers (D)) were bitten by resident workers. Each dot represents the performance of the tested individuals, the lines connecting the two males or two foreign workers used in the same test reflect an increasing or decreasing performance from conspecific to heterospecific. The color of the lines and dots indicates the species of origin of the gyne. For panels A and B, the arrows identify males with at least one successful mating (all arrows are blue indicating that all successful males were from C. piliscapa).

ACKNOWLEDGMENTS

This work was funded by ATM grant (Action Thématique de Recherche) from the MNHN (Muséum National d'Histoire Naturelle) and AP grant (Action Ponctuelle) from the EPHE (Ecole Pratique des Hautes Etudes) to CD. We thank Neyla Turak and Tiphaine Boursier (master students) for their help with the morphological and genetic analyses, as well as Arnaud Suwalski for his help in field sampling.

DATA ARCHIVING DNA sequences: GenBank accession numbers MK684255 to 684308 (see Table S2). Microsatellite dataset has been deposited in the Open Science Framework database, https://osf.io (DOI 10.17605/OSF.IO/Q6485). Chemical, morphological and behavioral data will be made accessible via DRYAD after acceptance.

AUTHOR CONTRIBUTIONS

CD designed the studies. CD, TM and PCB collected samples. PCB and PF performed the genetic experiments, BF and RC performed the morphometric experiments, PdE, CL and TM performed the chemical experiments, BF and PCB performed the behavioral experiments. CD, PAE, BF, AK and SM analyzed the data and all co-authors participate in the writing of the paper, with a major contribution of PAE, BF and CD.

COMPETITING INTEREST

The authors declare having no competing financial interests in relation to the work described.