Noppol Kobmoo
email: noppol.kob@biotec.or.th

Suchada Mongkolsamrit

Nuntanat Arnamnart

Janet Jennifer Luangsa-Ard

Tatiana Giraud

Population genomics revealed cryptic species within host-specific zombie-ant fungi (Ophiocordyceps unilateralis)

Keywords:

The identification and delimitation of species boundaries are essential for understanding speciation and adaptation processes and for the management of biodiversity as well as development for applications. Ophiocordyceps unilateralis sensu lato is a complex of fungal pathogens parasitizing Formicine ants, inducing zombie behaviors in their hosts. Previous taxonomic works with limited numbers of samples and markers led to the "one ant-one fungus" paradigm, resulting in the use of ant species as a proxy for fungal identification. Here, a population genomics study with sampling on three ant species across Thailand supported the existence of host-specific species in O. unilateralis s.l. with no footprints of long term introgression despite occasional host shifts and first-generation hybrids. We further detected genetic clusters within the previously delimited fungal species, with each little footprints of recombination, suggesting high levels of inbreeding. The clusters within each of O. camponoti-leonardi and O. camponoti-saundersi were supported by differentiation throughout the genome, suggesting they may constitute further cryptic species parasitizing the same host, challenging the one ant-one fungus paradigm. These genetic clusters had different geographical ranges, supporting different biogeographic influences between the north/center and the south of Thailand, reinforcing the scenario in which Thailand endured compartmentation during the latest Pleistocene glacial cycles.

Introduction

The identification and delimitation of species boundaries are essential for understanding speciation and adaptation processes and for the management of biodiversity and conservation. The existence of cryptic species, i.e. distinct species with hardly distinguishable morphologies, have been increasingly documented in all phyla and constitutes a challenge for the study of speciation, the assessment of biodiversity and the implementation of conservation plans [START_REF] Fišer | Cryptic species as a window into the paradigm shift of the species concept[END_REF][START_REF] Struck | Finding evolutionary processes hidden in cryptic species[END_REF]. Particularly for fungi, the recognition of species can be difficult as they are versatile organisms, with variable morphology and life styles [START_REF] Khonsanit | Cryptic species within Ophiocordyceps myrmecophila complex on formicine ants from Thailand[END_REF][START_REF] Molina | Another example of cryptic diversity in lichen-forming fungi: The new species Parmelia mayi (Ascomycota: Parmeliaceae)[END_REF][START_REF] Mongkolsamrit | Disentangling cryptic species with isaria-like morphs in Cordycipitaceae[END_REF][START_REF] Montarry | Spatio-temporal distribution of Erysiphe necator genetic groups and their relationship with disease levels in vineyards[END_REF][START_REF] Vincenot | Out of Asia: biogeography of fungal populations reveals Asian origin of diversification of the Laccaria amethystina complex, and two new species of violet Laccaria[END_REF]. Fungi are frequently symbionts to other organisms and often evolve host specialization [START_REF] Fournier | Sympatric genetic differentiation of a generalist pathogenic fungus, Botrytis cinerea, on two different host plants, grapevine and bramble[END_REF][START_REF] Gerardo | Ancient host-pathogen associations maintained by specificity of chemotaxis and antibiosis[END_REF][START_REF] Gerardo | Exploiting a mutualism: parasite specialization on cultivars within the fungus-growing ant symbiosis[END_REF], which is a fertile background for host-specific, cryp tic species. The use of the phylogenetic species criterion based on concordance between gene genealogies has been highly useful for delimiting fungal cryptic species [START_REF] Taylor | Phylogenetic species recognition and species concepts in fungi[END_REF]. With the advent of next-generation sequencing, the concordance among gene genealogies can be assessed throughout the genome and new phylogenetic species can be revealed with even more power and reliability [START_REF] Matute | Fungal species boundaries in the genomics era[END_REF][START_REF] Weiss | Genome-wide single-nucleotide polymorphism data reveal cryptic species within cryptic freshwater snail species-the[END_REF] or well-established species can be questioned [START_REF] Lifjeld | When taxonomy meets genomics: lessons from a common songbird[END_REF]. Furthermore, a population genomics approach can identify gene flow and introgression at fine scale along the genome [START_REF] Chan | Species delimitation with gene flow: a methodological comparison and population genomics approach to elucidate cryptic species boundaries in Malaysian Torrent Frogs[END_REF][START_REF] Savary | A population genomics approach shows widespread geographical distribution of cryptic genomic forms of the symbiotic fungus Rhizophagus irregularis[END_REF].

Ophiocordyceps unilateralis sensu lato is a ubiquitous species complex of ant-pathogenic fungi with a pan-tropical distribution [START_REF] Araújo | Zombie-ant fungi across continents: 15 new species and new combinations within Ophiocordyceps I. Myrmecophilous hirsutelloid species[END_REF] and an outstanding biology; these fungi modify the behavior of their ant hosts to promote their own dispersal. Infected ants are called zombie ants because they develop erratic behaviors, excluding themselves from their nest before climbing into vegetation and biting vege tal substrates until death [START_REF] Hughes | Behavioral mechanisms and morphological symptoms of zombie ants dying from fungal infection[END_REF]. The taxon called O. unilateralis has historically been characterized by the presence of fungal stalks arising from dead ants with lateral cushions containing perithecia and whole ascospores. Recent studies based on morphology and molecular markers have shown that O. unilateralis represented a complex of cryptic species, with allopatric species from distinct continents [START_REF] Araújo | Zombie-ant fungi across continents: 15 new species and new combinations within Ophiocordyceps I. Myrmecophilous hirsutelloid species[END_REF] as well as host-specific sympatric species [START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF]. A single and host-specific fungal species appeared to be associated to each ant species, leading to the "one ant-one fungus" paradigm and to the use of ant host identity as proxy for fungal species identification. However, little attention has been paid to intra-specific diversity and to the possibility of differentiated populations or species occurring on a given ant species. Species in the O. unilateralis complex are closely related and many are sympatric, which could lead to introgression and hybridization, although these questions have never been comprehensively studied.

In Thailand, six cryptic species of the O. unilateralis s.l. complex have been identified based on molecular phylogenies and morphological traits [START_REF] Kobmoo | New species of Ophiocordyceps unilateralis, an ubiquitous pathogen of ants from Thailand[END_REF][START_REF] Luangsa-Ard | Ophiocordyceps halabalaensis: A new species of Ophiocordyceps pathogenic to Camponotus gigas in Hala Bala Wildlife Sanctuary, Southern Thailand[END_REF]: O. polyrhachis-furcata, O. camponoti-leonardi, O. camponoti-saundersi, O. septa and O. rami and O. halabalaensis. At the macro-morphological level, only O. rami and O. halabalaensis can be easily distinguished from the others, having one stalk arising from the junction of the ant head and thorax, and two others arising from the junction of the forelegs and thorax; O. rami has sinuous purple stromata while O. halabalaensis has bigger but less sinuous blackish-brown stromata. At the micro-morphological level, O. septa has distinctive swollen and multi-septate ascospores; O. polyrhachis-furcata, O. camponoti-leonardi and O. camponoti-saundersi are the most similar genetically and morphologically and can be found on the ants Polyrhachis furcata, Colobopsis leonardi and C. saundersi, respectively [START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF]. Although morphological traits can distinguish these species, host species is the most convenient proxy for species identification and appears reliable [START_REF] Araújo | Zombie-ant fungi across continents: 15 new species and new combinations within Ophiocordyceps I. Myrmecophilous hirsutelloid species[END_REF][START_REF] Kobmoo | New species of Ophiocordyceps unilateralis, an ubiquitous pathogen of ants from Thailand[END_REF][START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF]. The criterion of congruence among gene genealogies was consistent with the separation of these host-specialized cryptic species in O. unilateralis s.l. [START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF]; however, they were based on only two gene sequences (Elongation factor 1-alpha and Beta tubulin).

Due to frequent expansion and contraction events of the forest during the Pleistocene glacial cycles, Southeast Asia shows strong biogeographic compartmentation (Woodruff, 2010). Genetic differentiation between the north and the peninsular south have been reported in many organisms, including insects such as bees [START_REF] Insuan | Genetic differentiation of the giant honey bee (Apis dorsata) in Thailand analyzed by mitochondrial genes and microsatellites[END_REF][START_REF] Sittipraneed | Genetic differentiation of the honey bee (Apis cerana) in Thailand revealed by polymorphism of a large subunit of mitochondrial ribosomal DNA[END_REF] and flies [START_REF] Pramual | Phylogeography of the black fly Simulium tani (Diptera: Simuliidae) from Thailand as inferred from mtDNA sequences[END_REF], trees (Yu et al., 2019), vertebrates including frogs [START_REF] Chen | Large-scale phylogenetic analyses provide insights into unrecognized diversity and historical biogeography of Asian leaf-litter frogs, genus Leptolalax (Anura: Megophryidae)[END_REF], monkeys [START_REF] Karanth | Molecular phylogeny and biogeography of langurs and leaf monkeys of South Asia (Primates: Colobinae)[END_REF] and bats [START_REF] Tu | Comparative phylogeography of bamboo bats of the genus Tylonycteris (Chiroptera, Vespertilionidae) in Southeast Asia[END_REF]; so far this geographical pattern has not been studied in fungi. Furthermore, an increase in genetic diversity from North to South has been reported in various organisms in the Northern hemisphere, including fungi [START_REF] Miraldo | An Anthropocene map of genetic diversity[END_REF][START_REF] Smith | A latitudinal phylogeographic diversity gradient in birds[END_REF][START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF] although investigations in Southeast Asia are still lacking.

Furthermore, genetic diversity in fungi is also influenced by reproductive modes (sexual vs clonal) and mating systems (selfing vs outcrossing) [START_REF] Nieuwenhuis | The frequency of sex in fungi[END_REF]. The sexual reproductive structures of O. unilateralis has been widely observed in nature; the fungi develop ascomata from diseased ants, which release ascospores into the environment. However, additional transmission via conidia (asexual spores) cannot be ruled out and sexual reproduction can occur between closely related individuals, and even between clonemates in case these fungi are homothallic, which is unknown. Which reproductive mode and mating system are predominant, and whether there are differences between species of O. unilateralis, are still open questions that can have major implications for genetic diversity, reproductive barrier and adaptation in these fungi.

In the present study, we therefore used a population genomics approach to investigate: (i) whether the existence of host-specific cryptic species was supported by genome-wide genetic differentiation in O. unilateralis, (ii) in this case, whether introgression occurred between species, and in which genomic locations, (iii) whether cryptic genetic subdivision could be identified beyond host specific species, e.g. according to geographic regions, (iv) whether genetic diversity increased towards the South within each species, and (v) what the reproductive mode and mating system of the different species were. For these goals, we collected samples of O. polyrhachis-furcata, O. camponoti-leonardi and O. camponoti-saundersi (identified based on host identity as a first proxy) across Thailand, and sequenced their genomes. We investigated their species status using the criterion of concordance between genealogies performed using multiple genomic regions. We also studied the genetic subdivision within each species as well as possible clines in diversity, and the possibility of introgression between species through analyses of Bayesian clustering and allelic spectrum-based models. Finally, we also searched for footprint of recombination to gain insight into the mode of reproduction.

Materials and methods

Sampling, isolation and sequencing

A total of 59 samples of O. unilateralis s.l. were collected in Thailand, 29 on C. leonardi, 16 on C. saundersi and 14 on P. furcata (Fig. 1a and Table 1S); fungal samples were named according to their host species (i.e., O. polyrhachis-furcata, O. camponoti-leonardi and O. camponoti-saundersi, respectively). To conform to the Nagoya protocols on access and benefit-sharing, we obtained an authorization from the Department of National Parks, Wildlife and Plant Conservation (DNP) at the Ministry of Natural Resources and Environment of Thailand for all sample collections. Isolations and extractions of fungal DNA were conducted following [START_REF] Kobmoo | A genome scan of diversifying selection in Ophiocordyceps zombie-ant fungi suggests a role for enterotoxins in co-evolution and host specificity[END_REF]. Basically, ascospores were harvested after being discharged from the sexual reproductive structures (perithecia) of each sample on potato dextrose agar; whenever possible, random single ascospores were picked up with a needle and transferred to a liquid medium for culturing haploid strains. A single monospore strain per sample was retained for analyses. However, initiating cultures from single spores was not consistently successful; some strains were thus obtained by pooling the discharged ascospores together and thus resulted from a mixture of products of meiosis from a diploid individual; the corresponding DNAs thus correspond to diploid genotypes. Ascospores are thought to constitute the principal mode of infection for O. unilateralis [START_REF] Araújo | Diversity of entomopathogenic fungi. Which groups conquered the insect body?[END_REF] but it is still unknown whether infections by two mating-type compatible haploid individuals are required to allow for sexual reproduction or whether these fungi are homothallic.

DNA quality and quantity were checked with gel electrophoresis using 0.8% agarose gel and a Nanodrop (ThermoFisher). Paired-end genomic libraries (150 bp) were constructed and sequenced with an Illumina HiSeq3000 system at the GenoToul platform (Toulouse, France).

Mapping and detection of single nucleotide polymorphisms

The raw reads were trimmed for adaptors and low-quality bases on both ends (q < 20) with Trim Galore! (Krueger, 2015). Reads shorter than 50 bp were discarded. Processed reads were then initially mapped to the latest version of O. polyrhachis-furcata reference [START_REF] Kobmoo | A genome scan of diversifying selection in Ophiocordyceps zombie-ant fungi suggests a role for enterotoxins in co-evolution and host specificity[END_REF], thanks to its higher quality (43.25 Mb, 68 scaffolds, N50 ~2.9 Mb) compared to those of the two other species 531 scaffolds,N50 ~140 Kb; Mb, 1,700 scaffolds, N50 ~102 Kb), using bwa-mem (Li, 2013) -M. SNPs were detected using GATK [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF] through the GATK HaplotypeCaller function. The ploidy parameter was set as haploid for single-spore strains, and as diploid for multiple-spore strains as a bulk of spores from a sample represents a mix of products of multiple meioses from the same diploid individual. As there were no reference variants, the variant recalibration step was dismissed but the variants were filtered a posteriori with total depth between 700× and 5000×, depth for individual samples >10, mapping quality >10, quality normalized by depth >20, genotyping quality >50; various corrections for biases related to the difference of quality between reference and variant alleles were also applied. Any SNP overlapping with regions including repeats identified beforehand with RepeatMasker (Smit et al., http://www.repeatmasker.org) was excluded. Finally, only biallelic SNPs without indels and with at least 80% samples with reliable inferred variants were retained.

For the within-species analyses, the reads from O. camponoti-leonardi and O. camponoti-saundersi were mapped to their respective reference [START_REF] Kobmoo | A genome scan of diversifying selection in Ophiocordyceps zombie-ant fungi suggests a role for enterotoxins in co-evolution and host specificity[END_REF]. The same procedures of mapping and variant calling were applied as in O. polyrhachis-furcata above with slightly different filtering (for O. camponoti-leonardi, final SNPs with total depth between 50× and 3000×, depth for individual samples >20, mapping quality >10, quality normalized by depth >20; for O. camponoti-saundersi, final SNPs with total depth between 50× and 1000×, depth for individual samples >20, mapping quality >10, quality normalized by depth >20). Kobmoo et al. Molecular Phylogenetics and Evolution xxx (xxxx) xxx-xxx

U N C O R R E C T E D P R O O F N.

Phylogenetic reconstruction and genealogical concordance phylogenetic species recognition (GCPSR)

Whole scaffold sequences for each sample were obtained using the FastaAlternateReferenceMaker tool in GATK. Sequences of a given scaffold from different samples altered at the SNPs positions were aligned using an in-house awk script. The alignments of all scaffolds were concatenated to generate a whole-genome alignment. The alignments were used to build maximum-likelihood phylogenetic trees with RaxML based on the model GTRCAT with 1000 bootstrap replicates. A tree for each scaffold with variants was inferred with the same parameters. All trees were rooted with O. kimflemingiae [START_REF] De Bekker | Ant-infecting Ophiocordyceps genomes reveal a high diversity of potential behavioral manipulation genes and a possible major role for enterotoxins[END_REF], by aligning this latest genome to that of O. polyrhachis-furcata using the nucmer command from MUMmer 3 [START_REF] Delcher | Fast algorithms for large-scale genome alignment and comparison[END_REF][START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]. The concordance between different genomic regions was also evaluated under the framework of Bayesian concordance analyses using BUCKy v. 1.4.4 (Larget et al., 2010), after Bayesian inferences via MrBayes 3.2.2 [START_REF] Ronquist | MrBayes 3.2: efficient bayesian phylogenetic inference and model choice across a large model space[END_REF]) over all scaffolds with SNPs bigger than 100 Kb chunked by 100 Kb separately. The MCMC chains were allowed to sample across the substitution model space. Chunks of 100 Kb were used due to the high linkage disequilibrium found (see the Section 3.7). The posterior tree distributions from the Bayesian inferences were fed into BUCKy to estimate the concordance of the nodes among the tree topologies, given as a concordance factor (CF) ranging from 0 to 1 (absence of concordance to complete concordance). The concordance factor thus represents the genome-wide genealogical support for the monophyly of clades.

Linkage disequilibrium and recombination footprints

The correlation coefficient r 2 between SNPs of less than 1 Mbp with minor allele frequency of at least 0.1 and per-sample inbreeding coefficients for diploid samples were calculated with plink [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF]; r 2 is a correlation coefficient representing the association and thus linkage between SNPs. In case of regular recombination, physically close SNPs are expected to be linked thus to have higher r 2 than farther SNPs. In case of low frequency of sexual reproduction, linkage disequilibrium, and thus r 2 , will be high across the whole genome. The r 2 coefficient was averaged over 100 bp intervals and plotted against the pairwise physical distance to obtain linkage disequilibrium (LD) decay plots. Pairwise Homoplasy Index (PHI) tests of recombination [START_REF] Bruen | A simple and robust statistical test for detecting the presence of recombination[END_REF] were performed for the whole dataset and for each species on each scaffold. The test is based on the detection of incompatibility between pairs of sites regarding whether there is a genealogical history that can be inferred parsimoniously that does not involve any convergent mutations; if any incompatibility exists, under an infinite-site model of evolution, it must come from recombination. PHI-test returns the probability of observing the data under the null hypothesis of no recombination. Bonferroni corrections were applied to the p-values. Population-scaled recombination rates (rho = 2Ne × r; Ne = effective population size, r = recombination per site per generation) were calculated for each species using the R package FastEPRR [START_REF] Gao | New Software for the Fast Estimation of Population Recombination Rates (FastEPRR) in the Genomic Era[END_REF]. The per-sample inbreeding coefficients were calculated for diploid individuals as method-of-moments estimate F coefficients based on observed and expected homozygosity for each diploid sample.

Neighbor-net networks were inferred using SplitsTree [START_REF] Huson | Application of phylogenetic networks in evolutionary studies[END_REF] for the whole dataset and for each species to gain insights into the level of recombination between and within the various O. unilateralis species.

Principal component analysis (PCA), molecular diversity and genetic differentiation

The nucleotide diversity (Pi) of all species/populations and inter-specific/inter-population F ST and d xy values were calculated using PopGenome package in R [START_REF] Pfeifer | PopGenome: an efficient swiss army knife for population genomic analyses in R[END_REF]. These calculations were performed based on SNPs called using the O. polyrhachis-furcata reference to use a common set of SNPs. The correlation between Pi and the latitude of the populations was tested with a Spearman's correlation test. We investigated the pattern of isolation by distance (IBD) by testing the existence of a correlation between F ST /(1 -F ST ) and the log of geographic distance [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF], using a Mantel test with 10,000 permutations.

Genetic clustering, inference of gene flow and introgression

A Bayesian clustering was conducted via FastStructure [START_REF] Raj | fastSTRUCTURE: variational inference of population structure in large SNP data sets[END_REF] with admixture model for the whole dataset to infer genetic structures between species. At the intra-specific level, we conducted the same analysis using the SNPs generated from each species reference. The most relevant number of clusters (K) were determined based on the plots of marginal likelihood and visual assessment of the barplots. The introgression footprints along the genome were investigated with a clustering analysis of genetic similarity on 500 SNPs sliding windows among the samples based on Gower's distance, using Mclust package in R [START_REF] Scrucca | mclust 5: clustering, classification and density estimation using Gaussian finite mixture models[END_REF]. Neighbor-joining trees based on the F84 distance [START_REF] Kishino | Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in hominoidea[END_REF] were also constructed to illustrate the relationships among the different clusters.

To infer events of gene flow between species, we used TreeMix [START_REF] Pickrell | Inference of population splits and mixtures from genome-wide allele frequency data[END_REF] which estimates the most likely topology representing the evolutionary history between populations and mixing based on genome-wide allele frequency data. TreeMix compares the likelihood of models with different numbers of migration events between populations.

Results

Detection of single nucleotide polymorphisms (SNPs)

By using O. polyrhachis-furcata BCC54312 as the reference for mapping all genomes, we retained 688,264 SNPs in the whole dataset. Only a small proportion of them was shared between the three host species (1144 SNP ~0.166%, Fig. S1). Based on the phylogenetic species recognition analysis (see the Section 3.2), we excluded a sample from the host C. leonardi with an unclear placement (NK122ss). There were 443,670 SNPs in the remaining dataset. The samples classified as O. camponoti-leonardi had more SNPs (219,017) than the other species (92,133 for O. camponoti-saundersi and 53,552 for O. polyrhachis-furcata). By using the O. polyrhachis-furcata reference genome to detect SNPs in all species, we could expect to find higher diversity in O. polyrhachis-furcata due to a bias favoring detection of SNPs from this species. However, we observed in contrast fewer SNPs in O. polyrhachis-furcata, which is likely due to genuine lower genetic diversity in this species (see Section 3.4). When using their respective references to call SNPs within each species, we found as expected more SNPs than when using the O. polyrhachis-furcata reference in the two other species [START_REF] Taylor | Clonal reproduction in fungi[END_REF]301,966 SNPs for O. camponoti-leonardi and 159,449 SNPs for O. camponoti-saundersi). Kobmoo et al. Molecular Phylogenetics and Evolution xxx (xxxx) xxx-xxx

U N C O R R E C T E D P R O O F N.

Species delimitation using the phylogenetic criterion of species recognition

The maximum likelihood (ML) tree based on whole genomes supported the existence of the three host-specific species. Indeed, we found three clades supported by 100% bootstrap support corresponding to the three hosts sampled (Fig. 2a). The only exceptions were NK604, isolated from a C. saundersi ant but placed within the O. camponoti-leonardi clade, and NK122ss, isolated from C. leonardi but placed as a sister group to the O. camponoti-saundersi clade with full bootstrap support (Fig. 2a). The trees inferred from the different scaffolds displayed congruent topologies (Fig. 2b) except for the scaffolds 6. 2,16,18,20,21,25,29,36,41,53, 105 (Fig. S2). Some of these scaffolds however suffered from a lack of polymorphism and/or yielded trees with low bootstrap supports (Table 2S).

The Bayesian concordance analysis comparing topologies between genealogies built from 100 kb fragments showed high concordance factors (CFs above 0.9; Fig. 3a), further indicating that most genomic regions supported the three host-specific clades and thus the existence of distinct, host-specialized species. The Bayesian concordance analysis also confirmed the position of NK122ss as the sister group of O. camponoti-saundersi (Fig. 3a).

Inter and intra-specific genetic differentiation

Based on the SNPs inferred using O. polyrhachis-furcata as the reference, the Bayesian clustering also yielded clusters grouping samples according to host species, except again for NK604 (Fig. 1c). The sample NK122ss, from the host C. leonardi, had admixed membership up to K = 4 but gained a total membership to a distinct cluster from K = 5, where two other clusters were distinguished in this host (Fig. S3a), a pattern corroborated by PCA (Fig. 1b). The sample NK122ss may belong to a species different from O. camponoti-leonardi, as also suggested by the phylogenies, and was thus excluded from further analyses. No additional cluster was found at higher Ks (Fig. S3a). The marginal likelihood curve in fact reached a plateau at K = 5 (Fig. S3b). The sample NK604, isolated from the host C. saundersi, clustered with those from the host C. leonardi in the PCA (Fig. 1b), in agreement with the results on species delimitation, and was thus considered to belong to O. camponoti-leonardi in the following analyses. It is likely a spill over, i.e. an occasional infection from a non-specific species without sustainable persistence. Furthermore, two genetic clusters could be clearly distinguished within O. camponoti-leonardi (OCL1 and OCL2 in Fig. S1b and1c). The genetic differentiation between the three main species was high (F ST > 0.5, Table 1).

We inferred genetic clusters within each of the three main species based on SNPs called using specific reference per species. For O. camponoti-leonardi, the results were slightly different from those when the O. polyrhachis-furcata reference was used (see above), as we retrieved three clusters (OCL1, OCL2, OCL3) for most K values (Fig. 4a and Fig. S4a). Only at K = 6 were four clusters observed but the fourth disappeared at higher K values and the marginal likelihood was lowest at K = 6 (Fig. S4b). We therefore considered three genetic clusters in O. camponoti-leonardi; two clusters were parapatric-one in the extreme North (OCL1) and another one extended to the center (OCL2). The samples NK140 and NK561, which appeared to belong to OCL2 using SNPs called with the O. polyrhachis-furcata reference (Fig. 1b andc), were assigned differently when using SNPs called with the O. camponoti-leonardi reference; the sample NK140 displayed an admixed genetic membership but still mostly belonged to OCL2 while the sample NK561 from the south of Thailand was assigned to another cluster (OCL3).

Regarding O. camponoti-saundersi, three sympatric clusters (OCS1, OCS2, OCS3) were found from K = 3 and above, with a plateau reached at K = 3 for the marginal likelihood (Fig. S5); the three clusters were localized mainly in the north-east of Thailand (Fig. 4b). For O. polyrhachis-furcata, the marginal likelihood oscillated between K = 2 and K = 10 (Fig. S6b); we identified three main genetic clusters, one in the north-east of Thailand (OPF1) and two in the South (OPF2 and OPF3), OPF3 being however represented by a single strain (NK157) (Figs. 4c andS6a). Another individual from the South of Thailand (NK151) had different cluster assignments across the various K values (Figs. 4c andS6a) and was thus classified into none of the genetic clusters. NJ trees and PCAs retrieved the same patterns of genetic subdivisions within each of the three species (Figs. S7-S9). The intra-specific differentiation level (between clusters of a given named species) was overall lower than inter-specific differentiation level (Table 2). Neighbor-net networks using SplitsTree also recovered different species (Fig. 5) and the same well differentiated genetic clusters within species (Fig. S10).

Nucleotide diversity and isolation by distance

We used the SNPs inferred using the best assembled reference genome (O. polyrhachis-furcata) to compute nucleotide diversity (Pi) to compare between species (Table 1) as well as between intra-specific clusters (Table 2). The nucleotide diversity was highest for O. camponoti-leonardi (0.026), followed by O. camponoti-saundersi (0.015) and O. polyrhachis-furcata (0.007). Based on the Bayesian clustering, we considered O. camponoti-leonardi to be composed of two genetic clusters, OCL1 and OCL2 (excluding the strain NK561 that appeared to belong to a distinct cluster; Figs. S4 and S7), O. camponoti-saundersi to be composed of three genetic clusters (OCS1, OCS2 and OCS3), and O. polyrhachis-furcata to be composed of three clusters (OPF1, OPF2 and OPF3), excluding the genetically ambiguous individual NK151. The genetic diversity was variable among genetic clusters within each species (Pi values in Table 2). OCS1 had much lower diversity compared to OCS2 and OCS3, as shown by the Pi values and the shorter branches on the NJ tree (Fig. S8a). The correlation between genetic diversity (Pi) and latitude for O. camponoti-leonardi was not significant (Spearman's correlation test: r = -0.3, p-value = 0.683; Table S3 and Fig. S11a). Isolation by distance was not significant in O. camponoti-leonardi (Mantel's test: r = 0.236, p-value = 0.25; Fig. S11b). Not enough sampling sites were available for the other species to test for such correlations.

Concordance among genealogies for within-species clusters

We then investigated whether the different clusters found within each of the identified cryptic species so far could be considered as distinct species based on the criterion of concordance among multiple genomic genealogies. The Bayesian concordance analysis based on 100-kb DNA fragments showed that the three clusters in O. camponoti-saundersi formed respectively clades with high CFs (above 0.85; Fig. 3b). For O. camponoti-leonardi, the clade OCL1 was strongly supported with a high CF (0.993) while the clade OCL2 as defined previously was poorly supported (CF = 0.477; OCL2 sensu lato in Fig. 3b). However, when excluding NK140 from OCL2, the clade became strongly supported (CF = 0.901), while NK140, like NK561 which was previously shown to be admixed, had unclear status. For O. polyrhachis-furcata, the clusters OPF1 and OPF2 formed each a monophyletic clade but with low CFs (0.433 and 0.485 respectively: Fig. 3b). Fig. S11

Introgression and gene flow between species

Previous analyses showed that a few individuals had admixed genetic membership. We therefore performed genetic clustering analyses along the genome to identify possible footprints of localized introgression. The general pattern of host-specific clusters was found throughout the genome (Fig. 6). However, introgression footprints was found along the genome for two individuals, NK151 and NK157, that showed admixed or distinct membership in comparison to the other O. polyrhachis-furcata samples in Bayesian clustering. The O. camponoti-leonardi individual NK561 appeared here highly admixed between OCL1 and OCL2. We did not recover in this analysis the different genetic clusters in O. camponoti-saundersi (Fig. 6). The Gower's genetic distance on which this clustering is based may not be able to detect finer genetic differences. SplitsTree analyses showed little footprint of gene flow between species (Fig. 5) and recovered the same well differentiated genetic clusters within species without footprints of extensive gene flow, except perhaps between OPF1 OPF2, and the two highly admixed individuals NK151 and NK157 (Fig. S11).

Gene flow events were inferred using Treemix based on the allelic frequency spectra of the different species and clusters retrieved from the SNPs called on the O. polyrhachis-furcata reference genome. We excluded the individual with unclear status, i.e., NK122ss. To account for linkage disequilibrium, blocks of 500 SNPs were used. Setting the number of migration events to one, we inferred one gene flow event between OCL1 and OCS (Fig. 7). However, the model with gene flow was not significantly more likely than the model without migration (log likelihood ratio test: p-value 0.235).

Linkage disequilibrium and recombination footprints

We investigated recombination footprints within species and genetic clusters by looking at linkage disequilibrium (LD) decay (Fig. 8). None of the species had their LD decaying to less than 0.2 at the distance of 1 Mb, and the LD of both O. camponoti-leonardi and O. camponoti-saundersi remained flat after an initial rapid decrease, resembling patterns in clonal species [START_REF] Sepúlveda | Genome sequences reveal cryptic speciation in the human pathogen Histoplasma capsulatum[END_REF]. Using half of the maximum LD as the criteria for linkage distance [START_REF] Taylor | Clonal reproduction in fungi[END_REF], O. unilateralis species seemed to maintain linkage over more than 1 Mb. Nevertheless, the PHI-tests revealed significant deviations from the expectation of no recombination for most scaffolds in O. camponoti-leonardi and O. saundersi, and their respective intra-specific genetic clusters. Evidence of recombination was found on fewer scaffolds for O. polyrhachis-furcata, both when considering the whole species and considering separately its intra-specific clusters (Table S4). For all species and their clusters, no footprints of recombination could be recovered from small scaffolds with little polymorphism (Table S4). The population-scaled recombination rates (rho) were all of similar order of magnitude (1.88e-5 for O. polyrhachis-furcata, 1.49e-5 for O. camponoti-leonardi and 5.95e-5 for O. camponoti-saundersi). SplitsTree analyses also suggested low levels of recombination within species (Fig. 5) and within genetic clusters (Fig. S11); reticulations were observed but seemed overall limited compared to the whole network except for O. polyrhachis-furcata, in which the reticulations appeared to be restricted to the two admixed samples (NK151 and NK157) (Fig. S11b). The per-sample inbreeding coefficients of the diploid samples were extremely high in all species, approaching 1; this means that most sites were homozygous (Table S5), suggesting a very high level of inbreeding. Nevertheless, one individual was found in an ant species different from the typical host of its species, as previously documented [START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF], and some others showed signals of genetic admixture in Bayesian clustering which were also reflected into ambiguous or poorly supported placements in phylogenetic trees. Using an allelic frequency spectrum-based inference, a single event of gene flow from O. camponoti-leonardi to O. camponoti-saundersi was detected; however, the model with gene flow was not statistically more likely than the model without gene flow. Occasional spill-over or hybridization thus did not seem to lead to long term introgression among species, potentially due to selection against hybrids.

Possible additional cryptic species and the one ant-one fungus paradigm

Our population genomics study further revealed the existence of differentiated genetic clusters within each host-specific species. The clusters within O. camponoti-leonardi and O. camponoti-saundersi were particularly well differentiated and supported through different analyses, with little footprints of gene flow. Multiple genealogies were congruent in separating the genetic clusters within camponoti-saundersi and O. camponoti-leonardi indicating they should be considered as distinct species according to the GCPSR criterion [START_REF] Taylor | Phylogenetic species recognition and species concepts in fungi[END_REF], and not just as different populations, and they would deserve further morphological examinations in future studies. The clusters within O. polyrhachis-furcata in contrast did not appear supported as distinct species.

This study thus revealed hidden diversity beyond well described species, which could not be detected previously [START_REF] Kobmoo | New species of Ophiocordyceps unilateralis, an ubiquitous pathogen of ants from Thailand[END_REF][START_REF] Kobmoo | Molecular phylogenies reveal host-specific divergence of Ophiocordyceps unilateralis sensu lato following its host ants[END_REF], probably due to the limited quantity of markers and ex amined samples [START_REF] Balasundaram | How many DNA markers are needed to reveal cryptic fungal species?[END_REF][START_REF] Taylor | Phylogenetic species recognition and species concepts in fungi[END_REF]. In addition, a sample (NK122ss) did not cluster with other individuals parasitizing the same host species and was instead placed as a sister taxon to O. camponoti-saundersi, without presenting admixture footprints, suggesting that it could belong to a different species or be admixed with another species not present in the sample. The finding that there could be more than one Ophiocordyceps species parasitizing a given host challenges the well-accepted "one ant-one fungus" paradigm [START_REF] Araújo | Zombie-ant fungi across continents: 15 new species and new combinations within Ophiocordyceps I. Myrmecophilous hirsutelloid species[END_REF]. It should however be noted that the ant taxonomy is itself still debated. The ants Colobopsis leonardi and C. saundersi were previously classified in the genus Camponotus before a recent resurrection of the genus Colobopsis [START_REF] Ward | A revised phylogenetic classification of the ant subfamily Formicinae (Hymenoptera: Formicidae), with resurrection of the genera Colobopsis and Dinomyrmex[END_REF] and both are cryptic species with unclear taxonomy [START_REF] Laciny | Colobopsis explodens sp. n., model species for studies on "exploding ants" (Hymenoptera, Formicidae), with biological notes and first illustrations of males of the Colobopsis cylindrica group[END_REF]. The genetic clusters found within described Ophiocordyceps species in our study could actually be associated to different ant species in agreement with the one ant-one fungus hypothesis. A taxonomic work with a strong population genomics aspect on both the hosts and the parasites will be necessary to elucidate species delimitation and cophylogeny analyses to assess whether co-divergence occurred [START_REF] De Vienne | Cospeciation vs host-shift speciation: methods for testing, evidence from natural associations and relation to coevolution[END_REF]Morand et al., 2015).

Reproductive mode and mating system in Ophicordyceps unilateralis

Recombination likely occurs in O. unilateralis sensu lato, but probably very rarely. The rho values were of the same order of magnitude (-5) for all three host-specific species and much lower than in fungi with regular outcrossing such as the heterothallic Zymoseptoria tritici and Z. ardabilae (rho values being respectively 0.0217 and 0.0045: [START_REF] Stukenbrock | Fine-scale recombination maps of fungal plant[END_REF], or other fungi with diploid-selfing obligate sexual reproduction like Microbotryum species (order of magnitude -3 to -2: [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF] and Saccharomyces paradoxus (0.001-0.003: [START_REF] Tsai | Population genomics of the wild yeast Saccharomyces paradoxus: quantifying the life cycle[END_REF]. The per-sample inbreeding coefficients for the diploid individuals were actually all very high in all species and clusters, indicating that most SNPs were homozygous. As all the strains used in this study were isolated from sexual spores resulting from meiosis, O. unilateralis regularly undergoes sexual reproduction. Altogether, these findings suggest for the whole species complex a high level of inbreeding or even frequent intra-haploid mating, i.e. mating between genetically identical haploid cells [START_REF] Billiard | Sex, outcrossing and mating types: Unsolved questions in fungi and beyond[END_REF][START_REF] Giraud | Mating system of the anther smut fungus Microbotryum violaceum: selfing under heterothallism[END_REF]. Intra-haploid mating is only possible in homothallic fungi as the existence of separated mating types in het erothallic fungi prevents such mating. Whether O. unilateralis is homothallic or heterothallic is still unclear; the genome of O. polyrhachis-furcata contains at least a mating type locus (MAT-1-2) and a putative mating type-switching gene (Wichadakul, unpublished data) which suggests that it may be homothallic but no experimental observation has been done.

Such high inbreeding is consistent with the life history traits of these fungi. Fungal development into the transmission stage inside ant nest was not observed and infection outside the nest is likely favored to avoid social immunity [START_REF] Loreto | Long-term disease dynamics for a specialized parasite of ant societies: a field study[END_REF]. Field observations and a theoretical model established that effective infections were rare; only a weak proportion of samples observed on field went through sexual reproduction [START_REF] Andersen | Disease dynamics in a specialized parasite of ant societies[END_REF][START_REF] Mongkolsamrit | Life cycle, host range and temporal variation of Ophiocordyceps unilateralis/Hirsutella formicarum on Formicine ants[END_REF]. The probability to encounter compatible spores from distinct diploid fungal individuals may be low, thus promoting selfing in the fungi. Low probability of outcrossing can be found in other pathogenic fungi such as the plant castrating anther-smut fungi Microbotryum in which mating occurs between cells resulting from the same meiosis issued from a dispersing diploid spore [START_REF] Giraud | Mating system of the anther smut fungus Microbotryum violaceum: selfing under heterothallism[END_REF][START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF]. The infective spores of O. unilateralis are haploid which could render the probability to find a compatible mate in an insect even lower, favoring intra-haploid mating.

The biogeography of Ophiocordyceps unilateralis in Thailand

While the genetic clusters detected within O. camponoti-saundersi were sympatric and all restricted to the Northeast of Thailand, O. camponoti-leonardi and O. polyrhachis-furcata had widely distributed genetic clusters, with different geographical ranges, mainly reflecting the separation between the north/center and the south of Thailand. Such genetic differentiation between the North and the South of Thailand has been documented in many organisms [START_REF] Chen | Large-scale phylogenetic analyses provide insights into unrecognized diversity and historical biogeography of Asian leaf-litter frogs, genus Leptolalax (Anura: Megophryidae)[END_REF][START_REF] Karanth | Molecular phylogeny and biogeography of langurs and leaf monkeys of South Asia (Primates: Colobinae)[END_REF][START_REF] Pramual | Phylogeography of the black fly Simulium tani (Diptera: Simuliidae) from Thailand as inferred from mtDNA sequences[END_REF][START_REF] Sittipraneed | Genetic differentiation of the honey bee (Apis cerana) in Thailand revealed by polymorphism of a large subunit of mitochondrial ribosomal DNA[END_REF][START_REF] Tu | Comparative phylogeography of bamboo bats of the genus Tylonycteris (Chiroptera, Vespertilionidae) in Southeast Asia[END_REF]Yu et al., 2019), which was suggested to reflect different biogeographic entities, one from the Indochina and another from the Malay peninsula (Woodruff, 2010). The existence of northern and southern distinct biogeographic compartments also explains that the genetic diversity does not increase toward south, contrary to what is found in the northern hemisphere for many organisms [START_REF] Hewitt | Genetic consequences of climatic oscillations in the Qua-ternaryPhilosophical Transactions of the Royal Society of London[END_REF]. Enlarging the sampling of both the hosts and the pathogens across Southeast Asia would give novel insights into the host-pathogen phylogeography of this region.

Conclusion

In this study, population genomics supported the host-specificity of previously described cryptic species in ant pathogenic fungi and further revealed genetic clusters sharing the same host but with contrasting geographical ranges as candidates to be erected as new species in the zombie ant fungal complex O. unilateralis sensu lato. Whole genome analyses also revealed individuals with admixed genetic membership but not leading to persistent gene flow among species or genetic clusters. The recognition of cryptic species has direct impacts on species diversity assessment and ecological research and population genomics can be powerful to reveal hidden diversity, as shown in our study. More generally, population genomics can reveal the existence of cryptic species or different populations of pathogens, contributing to our understanding of host-parasite/pathogen specificity and phylogeography. Furthermore, we studied for the first time recombination footprints in O. unilateralis sensu lato which suggested limited recombinations and a high level of inbreeding. Kobmoo et al. Molecular Phylogenetics and Evolution xxx (xxxx)

U N C O R R E C T E D P R O O F U N C O R R E C T E D P R O O F U N C O R R E C T E D P R O O F N.

Fig. 1 .

 1 Fig. 1. (a) Sampling sites of Ophiocordyceps unilateralis sensu lato across Thailand (1 = Ban Huai Baba, 2 = Doi Inthanon National Park, 3 = Nam Nao National Park, 4 = Phu Khiao Wildlife Sanctuary, 5 = Khao Yai National Park, 6 = Ko Chang Island, 7 = Khlong Nakha Wildlife Sanctuary, 8 = Khao Luang National Park, 9 = Hala Bala Wildlife Sanctuary). The pie charts represent the proportions of fungal samples collected on different ant species (CL red = Colobopsis leonardi, CS green = C. saundersi, PF blue = Polyrhachis furcata). (b) Principal component analysis (PCA) based on whole genome single nucleotide polymorphisms (SNPs); colors correspond to the host species as in Fig. 1a. (c) Bayesian clustering based on whole genome SNPs (K = 5). Vertical bars represent individuals and the colors their assignment probabilities to five genetic clusters. The arrows indicate the samples NK122ss and NK604. Numbers below the barplots represent the sample sites as in Fig. 1a. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig. 2. Maximum likelihood-based tree inferred using (a) the whole genome altered at the SNP positions, (b) separate scaffolds. The colors at the tips represent the identities of ant hosts (red = Colobopsis leonardi, green = C. saundersi, blue = Polyrhachis furcata). Thick branches correspond to those leading to nodes with 100% bootstrap support. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

 1 Genetic differentiation between species (F ST : above the diagonal, d xy : below the diagonal; OCL = Ophiocordyceps camponoti-leonardi, OCS = O. camponoti-saundersi, OPF = O. polyrhachis-furcata). The nucleotide diversity (Pi) within species is given in the diagonal in bold; d xy = absolute divergence between two species; F ST = fixation index measuring relative divergence.

Fig. 4 .

 4 Fig. 4. Barplots of Bayesian clustering using single nucleotide polymorphisms (SNPs) inferred using conspecific genome references, with the geographical location of samples assigned to the different clusters, for each of the three species at K = 3: (a) Ophiocordyceps camponoti-leonardi, (b) O. camponoti-saundersi and (c) O. polyrhachis-furcata.

Fig. 5 .

 5 Fig. 5. A neighbor-net phylogenetic network of Ophiocordyceps unilateralis (OCL = O. camponoti-leonardi, OCS = O. camponoti-saundersi, OPF = O. polyrhachis-furcata).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Clustering of Ophiocordyceps unilateralis sensu lato, based on Gower's distance, along the O. polyrhachis-furcata BCC54312 reference genome using 500 SNPs sliding windows (OCL = O. camponoti-leonardi, OCS = O. camponoti-saundersi, OPF = O. polyrhachis-furcata). The red asterix ( ) indicates the sample NK561 with a footprint of admixture between OCL1 and OCL2. The red arrows ( ) point to the samples NK151 and NK157 of OPF which have a clear signature of admixture. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Linkage disequilibrium (LD) decay of Ophiocordyceps camponoti-leonardi (red), O. camponoti-saundersi (green) and O. polyrhachis-furcata (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2

 2 Genetic differentiation between genetic clusters (F

	Fst d xy	OCL1	OCL2	OCL3	OCS1	OCS2	OCS3	OPF1	OPF2	OPF3
	OCL1 OCL2 OCL3 OCS1 OCS2 OCS3 OPF1 OPF2 OPF3	0.011 0.144 0.175 0.262 0.261 0.268 0.159 0.164 0.204	0.711 0.016 0.117 0.216 0.214 0.221 0.110 0.116 0.157	0.907 0.781 0 0.246 0.244 0.250 0.144 0.149 0.186	0.934 0.877 0.996 0.001 0.079 0.063 0.155 0.161 0.202	0.915 0.853 0.976 0.913 0.010 0.081 0.152 0.158 0.200	0.904 0.841 0.962 0.834 0.813 0.006 0.159 0.165 0.207	0.884 0.749 0.985 0.979 0.947 0.927 0.001 0.013 0.062	0.870 0.736 0.966 0.962 0.931 0.912 0.425 0.003 0.067	0.884 0.837 1.000 0.995 0.971 0.955 0.965 0.924 0

ST : above the diagonal, d xy : below the diagonal). The nucleotide diversity (Pi) within clusters are given in the diagonal in bold. OCL3 and OPF3 include each a single individual (respectively NK561and NK157).

  xxx-xxx case of the Ancylus fluviatilis species complex. Ecol. Evol. 8, 1063-1072. doi:10.1002/ ece3.3706. Woodruff, D.S., 2010. Biogeography and conservation in Southeast Asia: how 2.7 million years of repeated environmental fluctuations affect today's patterns and the future of the remaining refugial-phase biodiversity. Biodivers. Conserv. 19, 919-941. doi:10.1007/s10531-010-9783-3. Yu, H., Tian, E., Zheng, L., Deng, X., Cheng, Y., Chen, L., Wu, W., Tanming, W., Zhang, D., Compton, S.G., Kjellberg, F., 2019. Multiple parapatric pollinators have radiated across a continental fig tree displaying clinal genetic variation. Mol. Ecol. 28, 2391-2409. doi:10.1111/mec.15046.

Acknowledgements

We would like to thank Alodie Snirc for advice concerning DNA extraction, Antoine Branca for help with bioinformatic protocols; Kanoksri Tasanathai, Donnaya Thanakitpipattana, Wasana Noisripoom and Artit Khonsanit for help in samples collection; Navavit Ponganan for making the maps in this publication. We also thank the Department of National Parks, Wildlife and Plant Conservation of Thailand for their cooperation and support of our project.

Funding This work was supported by a Marie-Skłodowska Curie Action [grant number 655278] to NK, a Biotec Fellow's Research Grant [grant number P1950231] to JJL, and by the National Science and Technology Development Agency, Cluster and Program Management Office [grant number P1300854] to SM.

Appendix A. Supplementary material

The raw reads in this study were deposited at the NCBI Sequence Read Archive (SRA) under the accession PRJNA530395. The information regarding the fungal isolates used in this study were deposited at the NCBI BioSample under the accessions reported in the Table S1 of the Supplementary Materials. Supplementary data to this article can be found online at https://doi.org/10.1016/j.ympev.2019.106580.