Chromosomal dynamics in space and time: evolutionary history of Mycetophylax ants across past climatic changes in the Brazilian Atlantic coast

Fungus-farming ants of the genus Mycetophylax exhibit intra and interspecific chromosome variability, which makes them suitable for testing hypotheses about possible chromosomal rearrangements that endure lineage diversification. We combined cytogenetic and molecular data from Mycetophylax populations from coastal environments to trace the evolutionary history of the clade in light of chromosomal changes under a historical and geographic context. Our cytogenetic analyses revealed chromosomal differences within and among species. M. morschi exhibited three distinct karyotypes and considerable variability in the localization of 45S rDNA clusters. The molecular phylogeny was congruent with our cytogenetic findings. Biogeographical and divergence time dating analyses estimated that the most recent common ancestor of Mycetophylax would have originated at about 30 Ma in an area including the Amazon and Southern Grasslands, and several dispersion and vicariance events may have occurred before the colonization of the Brazilian Atlantic coast. Diversification of the psammophilous Mycetophylax first took place in the Middle Miocene (ca. 18–10 Ma) in the South Atlantic coast, while “M. morschi” lineages diversified during the Pliocene-Pleistocene transition (ca. 3–2 Ma) through founder-event dispersal for the Northern coastal regions. Psammophilous Mycetophylax diversification fits into the major global climatic events that have had a direct impact on the changes in sea level as well as deep ecological impact throughout South America. We assume therefore that putative chromosomal rearrangements correlated with increased ecological stress during the past climatic transitions could have intensified and/or accompanied the divergence of the psammophilous Mycetophylax. We further reiterate that “M. morschi” comprises a complex of at least three well-defined lineages, and we emphasize the role of this integrative approach for the identification and delimitation of evolutionary lineages.

main driving force in long-term chromosome evolution and are effective at minimizing genetic risks and increasing the potential for genetic divergence 14-16 . In contrast, in a population-based cytogenetic study with fungus-farming ants of the genus Mycetophylax, the role of centric fusions in the karyotypic diversification of those ant lineages was proposed 17 . Mycetophylax was formerly composed of three endemic species of coastal sand dune environments, known as Restinga 18,19 : M. conformis Mayr, 1884; M. morschi Emery, 1888; and M. simplex Emery, 1888 (see Fig. 1 for details of occurrence areas). Differences in both the number and structure of their chromosomes were found, including intraspecific variation within M. morschi, showing populations with two distinct karyotypes (2n = 26 and 2n = 30). It has therefore been suggested that M. morschi comprises a complex of cryptic species or even different lineages that have been subjected to consecutive CRs 17 . It is noteworthy that through phylogenetic analyses, members of the former Cyphomyrmex strigatus group, including Mycetosoritis asper and Mycetosoritis clorindae, were shown to form a well-defined clade with strong statistical support along with Mycetophylax sensu stricto. Thus, these species were transferred to the genus Mycetophylax, which currently comprises 21 species 20 .
The integration of cytogenetic data, phylogenetic trees, and molecular dating considerably improved the findings on evolutionary history within lineages as well as providing implications that perhaps could not be found independently. However, few studies have combined such approaches, but what has emerged from this integration are interesting and even unexpected results [21][22][23][24] . Divergence times and phylogenetic relationships can be estimated simultaneously under a Bayesian inference. A method for divergence time estimation and fossil calibration using a stochastic branching process and relaxed-clock model in a Bayesian framework was recently developed 25,26 . The Fossilized Birth-Death (FBD) process has the advantage of assuming that fossils and molecular sequences from extant species comprise the same diversification process [25][26][27] . On the other hand, under a molecular cytogenetic approach, such as Fluorescence in situ Hybridization (FISH) chromosome mapping, useful diagnostic characters for morphologically conservative species may be provided as variations in the number and/ or position of ribosomal DNA (rDNA) clusters or CR indications given by changes in chromosome number. Hence, the use of FISH data has become an important tool for describing and delimiting new taxa, especially  60 .
for insects [28][29][30] . In this way, we can accurately test the assumptions made previously on the likely chromosomal evolution within the Mycetophylax lineages.
Additionally, the endemism of the three former Mycetophylax species to the Atlantic coastal environments remains a peculiar and attractive feature since the correlation between chromosomal and genetic variability needs to be examined in the context of geographic distributions of the species 31,32 . The Brazilian Atlantic coastal area (often referred to as the Restinga ecosystem) can be physically and biologically defined as coastal plains of marine sedimentary origin 33 that were subjected to ecological and geological stress during the major climatic transitions associated eustatic change 34,35 . Global past climate changes have certainly had biogeographic impacts due to the effects of alterations in temperature and moisture on species 36,37 . Ancestral range estimation contributes to a better understanding of species distributions over evolutionary time scales, and it is possible to infer dispersion and/ or vicariance events from previous climatic and/or geological events in a given geographic area 38 . The Restinga ecosystem is a coastal environment within the Atlantic Forest biome that is floristically and geomorphologically heterogeneous. It is composed of a range of organisms that have intrinsically adapted to adverse conditions 39 . These habitats have been deeply influenced by the long-term effects of climate change that may have contributed to species diversification due to dynamism. Therefore, we hypothesized that Mycetophylax populations may have diversified through chromosomal changes triggered by the climatic events that have continuously affected their habitat.
In order to test the assumptions raised above, we correlated cytogenetic, molecular, and biogeographic data from the additional Mycetophylax species with the phylogenetic tree of the fungus-farming ants, and under a historical and geographical context, we aimed to describe the evolutionary relationships and chromosomal changes that raised the karyotype observed today. In addition, with this integrative approach, we aimed to decipher the taxonomic status of M. morschi, which is currently considered a single species, but which certainly represents a complex of independent evolutionary lineages.

Results
Karyotype determination and FISH mapping. The detailed cytogenetic analysis with Mycetophylax populations revealed the presence of intra-and inter-specific karyotype variations related to the number of diploid chromosomes (2n), the karyotypic formula (KF), and the fundamental number (FN). FISH signals with the 18S rDNA probe showed a stable number of rDNA clusters, denoting only one pair of chromosomes, but its localization and position along the chromosome differed between karyotypes (Fig. 2). The distribution pattern of the TTAGG (6) telomeric repeat was restricted to the terminals of both chromosome arms in all karyotypes and intercalated stronger and weaker signals at random (Fig. S1). No signals for interstitial telomeric sites (ITS) were detected. Lineage A. Colonies from the southern (Torres/RS) and southeast (Cabo Frio/RJ and Ilha Comprida/SP) Brazilian coast presented 2n = 26 chromosomes. This karyotype is composed of nine metacentric, three submetacentric, and one acrocentric chromosome pairs; the karyotype formula is KF = 9 M + 3SM + 1 A; and the fundamental number (FN) is 50. The rDNA clusters are located in the pericentromeric region of the long arms of submetacentric chromosome pair 11 (Fig. 2).
Lineage B. Colonies from the northeast Brazilian coast (Ilhéus/BA) presented 2n = 28 chromosomes. This unexpected karyotype is composed of ten metacentric and four submetacentric chromosome pairs with a karyotype formula of KF = 9 M + 5SM and an FN of 56. The rDNA clusters are located in the pericentromeric region of the short arms of submetacentric chromosome pair 7 (Fig. 2). Note that, unlike the others, this lineage does not present acrocentric chromosomes.
Lineage C. Colonies of both localities of the southern Brazilian coast (Torres/RS and Araranguá/SC) presented 2n = 30 chromosomes. This karyotype is composed of nine metacentric, five submetacentric, and one acrocentric chromosome pairs with a karyotype formula of KF = 9 M + 5SM + 1 A and an FN of 58. The rDNA clusters are located in the terminal region of the short arms of the acrocentric chromosome pair (Fig. 2). For the first time, the presence of two cytotypes of "M. morschi" was recorded at the same geographic area. It is therefore conclusive that they live sympatrically.

Molecular phylogeny.
Our phylogenetic analyses were based on Bayesian Inference (BI), which exhibited strong statistical support for most nodes (Figs. 2 and S2). We use the aligned dataset provided by Sosa-Calvo et al. 20 and inserted additional Mycetophylax sequences obtained in this work. The phylogenetic tree generated showed three well-supported clades of "M. morschi": the first one corresponded to the cytotype 2n = 26 (lineage A) (Bayesian posterior probability (PP) = 1.0), and the others diversifying from it formed a sister group of cytotypes 2n = 28 (lineage B) and 2n = 30 (lineage C) (PP = 0.96) (Fig. 2). These three lineages were grouped with unknown species of the former Cyphomyrmex strigatus complex, composing a clade with intriguing weak support. M. simplex and M. conformis fell in a sister group (PP = 0.75), which, together with another unknown species of the former C. strigatus complex, formed a major clade with "M. morschi" lineages and the unknown species from the C. strigatus complex (PP = 1.0) (Fig. 2). This major clade will henceforth be called "psammophilous Mycetophylax". FBD-based divergence dating. The FBD-based divergence dating analysis recovered the stem-group age (i.e., the earliest possible origin) of fungus-farming ants as 64.6 million years ago (Ma) and the crown-group age (i.e., the latest possible origin) as 61. Historical biogeography. The three independent estimates of the ancestral range showed similar results. In the BGB analysis, the BAYAREALIKE + j model yielded the best statistical fit for the data (AIC = 101.9; p = 0.06), although it had no significant difference to the model without the founder-event speciation parameter (Table S1). All models pointed out that the most recent common ancestor (MRCA) of Mycetophylax would have originated in an area including the Amazon and Southern Grasslands (BGB: p = 0.4, BayArea: p = 0.43, and BBM: p = 0.5) (Fig. 3). Most of the nodes estimated dispersal while vicariance was recovered in a few nodes, including the split of the "M. morschi" lineages.
The colonization of the coastal regions probably involved major dispersal events that would have directed the Mycetophylax ancestors to the Atlantic Forest from the Amazon and subsequently spread to the South Coast region at about 17 Ma (BGB: p = 0.37, BayArea: p = 0.59, and BBM: p = 0.62). All models strongly estimated the ancestral range of all psammophilous Mycetophylax, including the remaining ants from the "strigatus complex", which apparently lived on the Atlantic coast of Southern Brazil and later scattered farther north in both coastal and rainforest regions, including the Amazon (BGB: p = 0.98, BayArea: p = 0.89, and BBM: p = 0.89) (Fig. 3). No dispersion events were recovered, primarily indicating diversification within the area or sympatric range inheritance.   www.nature.com/scientificreports www.nature.com/scientificreports/ Mycetophylax showed remarkable chromosomal polymorphisms among populations 17 , which provoked us to seek additional knowledge about their evolutionary relationships. Our data corroborate the previous karyotypic findings 17 on the species M. conformis and M. simplex, both in terms of the number and structure of their chromosomes. Likewise, our data support the findings of dissimilar karyotypes within M. morschi populations with the unexpected discovery of a new karyotype (2n = 28) in a geographic area previously unknown to this species, the northeast Brazilian Atlantic coast. No universal rule seems to govern karyotype evolution, but a chromosomal variation between closely related lineages, such as psammophilous Mycetophylax, strongly suggests the role of CRs in the lineage diversification. In fungus-farming ants, few species have had their karyotypes described so far (45 out of 245), but available data show extensive variation among species (2n = 8 to 2n = 50) 4 . This aspect, which appears to be typical of Formicidae, also corresponds to the clade of the fungus-farming ants and may be the outcome of diverse CRs, such as Rb fusion/fission.
In an attempt to partially understand the speciation history of the psammophilous Mycetophylax, we integrated cytogenetic and phylogenetic data. Our cytogenetic information has provided some evidence that psammophilous Mycetophylax comprises distinct and new evolutionary lineages, since different chromosomal numbers eventually do not match during meiosis, acting as post-zygotic barriers. Besides, as rDNA is considered one of the most conservative fraction from the eukaryotic genome 40,41 , variations in the location and chromosome bearer of the rDNA clusters, as observed in the karyotypes analyzed here, may be related to the differentiation and divergence of natural populations 42 . In fact, rDNA clusters are even considered hotspots of chromosomal changes and lineage diversification in mice 43 . Here, we found that all of the Mycetophylax karyotypes and cytotypes showed only two signals of 45S rDNA but at distinct locations and homologue pairs suggesting that changes in ribossomal cluster appears to be related to the diversification of these ants. Some chromosomal mechanisms have been suggested to explain such mobility of rDNA clusters within sister lineages, including ectopic recombination 44 , transposition 45 , and centric fission 46 . The fact that there is no variation in the number of these clusters can be attributed to homologous recombination between repetitive sequences dispersed throughout genome (i.e., ectopic recombination) or more likely to transposition of these sequences into arbitrary position and chromosomes. However, we cannot state any further inferences, only that the intraspecific chromosomal diversity in M. morschi gives us further indications that these lineages may indeed be considered distinct species.
The application of molecular phylogenetics on the study of karyotype evolution has enabled us to deduce the direction of changes that give rise to chromosome variation in the most objective way 47,48 . For example, just because two species of a genus share the same chromosome number does not mean that they are phylogenetically closer to each other than two species with different chromosome numbers 5 . However, the base chromosome number (i.e., the haploid number present in the initial population of a monophyletic clade) may be directly related to the chromosomal variability inside that clade and its sister group 5 . Our phylogenetic reconstruction is congruent with our cytogenetic findings showing that the Mycetophylax cytotypes are phylogenetically close. The lineage A diverged prior to the other "M. morschi" lineages, which would imply that 2n = 26 (n = 13) was probably the ancestral karyotype. However, our results on chromosome evolution reconstructed its ancestral karyotype as 2n = 28 (n = 14), whereas it was 2n = 32 (n = 16) for the clade M. conformis + M. simplex. Previous estimates have shown slightly contrasting results, with an ancestral diploid number for "M. morschi" lineages (n = 15) and the clade M. conformis + M. simplex (n = 17), leading the authors to propose fusion as the main mechanism of chromosome evolution in Mycetophylax 17 .
Our FISH-based comparative approach indicates that is more likely that both Rb fissions and fusions contributed to the karyotypic diversification in psammophilous Mycetophylax and similarly in Formicidae rather than a larger role for Rb fissions, as has been proposed [14][15][16] . In fact, specifying which chromosomes are involved in such rearrangements is a difficult task, especially since large-scale CRs can modify karyotypes in a striking way. For example, a karyotypic reduction observed in Cucumis species (from 2n = 24 to 2n = 14) involved at least 59 CRs, including five fusions, four translocations, and 50 inversions 49 . Not least, the presence of only a single pair of chromosomes bearing the 45S rDNA clusters in fungus-farming ants appears to be uniform. Taking into account the phylogenetic relationships of fungus-farming ants from the basal species Mycocepurus goeldii 50 and cross referencing our data with the most derived species of leafcutter ants from Acromyrmex and Atta genera 51,52 , we noticed that there are only a pair of chromosomes bearing the 45S rDNA cluster, which would imply that this chromosomal feature may be plesiomorphic.
Fusions can be directly observed using telomeric probes. After in tandem fusions, the canonical telomeric sequences may remain in the interstitial chromosome regions, serving as clues of fusion and inversions involving the telomeres. Telomeric sequences can either be eliminated or not be eliminated by chromosome breakage, and when complete elimination does not occur, interstitial telomeric sites (ITS) can be observed 53 . Considering the fusion hypothesis proposed by Cardoso et al. 17 as the mechanism shaping Mycetophylax karyotypes, signals for ITS would be expected, but we did not find any evidence for ITS, suggesting that if tandem fusion did take place the telomeric regions would have been lost. Lastly, we added more Formicidae species containing the TTAGG as the ending motif on their chromosomes, and thus we reinforce the hypothesis that this sequence could have been the ancestral motif of the ant telomeres 54,55 . evolutionary history of psammophilous Mycetophylax under a cytogenetic perspective. Since the beginning of the Cenozoic (ca. 65 Ma), the Earth's biota has experienced both rapid and gradual climate changes 34,35 that have had profound impacts on species diversification and composition 56 . Indeed, the major climate transitions has triggered higher rates of species turnover, and the Paleocene-Eocene Thermal Maximum "PETM" (ca. 55 Ma), the Terminal Eocene Event "TEE" (ca. 34 Ma), Middle Miocene Climatic Optimum "MMCO" (ca. , and Pliocene-Pleistocene intervals (from ca. 3 Ma) have been reported as striking events for mammals 56 . Coincidentally, our time-calibrated phylogenetic reconstruction suggested that the basal divergence of the genus Mycetophylax took place in the Eocene-Oligocene boundary (ca. 38-22 Ma), while psammophilous Mycetophylax diverged into the Middle Miocene (ca. 18-10 Ma) and thereafter diversified into the existing lineages. Both inferred ages fit with one of these major events mentioned above, which makes them relevant to other taxa than mammals, including here invertebrates.
The fungus-farming ants originated about 60 Ma in South America 57 , likely during the post-extinction event recovery period and shortly before the PETM. After that, an abrupt drop in temperature began at around 52 Ma and continued to the end of the Eocene (i.e., TEE) 34 that had left hallmarks characterized by cooler and more seasonal climates, increase topographic heterogeneity, and persistent phylogenetic and ecologic diversification of lineages 58 . It was suggested therefore that this global cooling event spurred both fungus-farming ants' lineages diversification and ant-fungus coevolution 57 . Likewise, this assumption agrees with our estimated age for diversification of the genus Mycetophylax in the Miocene (ca. 23-5.3 Ma). The Miocene was a climatically dynamic interval from Cenozoic, when long-term global climatic cooling was punctuated by extreme climatic optima (i.e., MMCO) 59 . After this period of global warming -the warmest one in the last 24 Ma 34 -a climate transition associated with major Antarctic ice sheet expansion and global cooling began, the so-called Middle Miocene Glaciation "MMG" 59 , which drastically extended to the Pleistocene glacial intervals. As these global climate oscillations were directly associated with marine transgressions and regressions 35 , the evolutionary history of psammophilous Mycetophylax on sand dunes could plausibly be related to the heterogeneity of the Brazilian Atlantic coast as well to their topographic and climatic gradients.
The basal split between M. conformis and M. simplex took place in the Middle to Late Miocene, ca. 10 Ma (HPD = 14.3-5.6 Ma), corresponding to the initial cooling period (i.e., MMG) and a consequent drop in sea levels reaching near to the current levels 35 . These species are morphologically distinguishable 18 and are not found to inhabit the same coastline area, except in Cabo Frio, Rio de Janeiro 60 . M. conformis is further distributed in the north of the Atlantic coast with the southern limit of its occurrence at about 24°S, while M. simplex occurs from the southern coast of São Paulo (~25°S) to the border of southern Brazil (~34°S) 60 . Our ancestral range estimates place the species ancestor as having inhabited the South Atlantic Coast, and we identified a major dispersal event. www.nature.com/scientificreports www.nature.com/scientificreports/ The putative biogeographic scenario would be that M. simplex inherited the ancestral range, while M. conformis would have colonized the North Atlantic Coast via founder-event dispersal. Founder-event speciation refers to changes in the lineage distribution during split (cladogenesis), where one daughter lineage inherits the ancestral range and the other lineage "jumps" to a new area 61,62 . Considering the current contrasting geographic range and respective karyotypes of M. simplex and M. conformis, the fixation of different CRs either by genetic drift or geographic isolation could have been triggered by the recurrent climatic events that continuously changed the environment, leading CRs to play some role in lineage diversification. Moreover, relict populations of M. simplex living sympatrically with M. conformis could further indicate that species were sympatric in the past and that during successive Pleistocene climatic transitions they could have become allopatric majority. Phylogeographic analyses with M. simplex populations corroborate this, as it was identified that the species maintained a stable population size until about 75,000 years ago when it underwent a gradual demographic expansion 63 .
Continuous climate changes have been more abrupt since the Early Pliocene leading to extreme periods of glaciation (most temperate regions) and aridification (tropical regions) later in the Pleistocene 58 . Particularly in South America, these fluctuations were responsible for chorological changes that expanded and contracted the distribution areas of taxa, communities, and biomes 36 . Indeed, increasing climatic oscillations during the past 3 Ma have driven many range changes for a range of taxa, and higher speciation rates were identified during the Pliocene and Pleistocene 64,65 , which may well be associated with chromosomal changes promoting lineage diversification and potential reproductive isolation 32,66 . The chromosomal diversity observed in psammophilous Mycetophylax is in agreement with this scenario, especially considering "M. morschi" lineages. The divergence of these lineages was recovered at the Pliocene-Pleistocene boundary (ca. 3.6-1.8 Ma). Historical biogeography analyses have estimated the ancestral range as the South Atlantic Coast with no dispersion events, suggesting that at least the initial diversification of this group may have occurred sympatrically. Interestingly, it has been argued that sympatric sister taxa have a higher average karyotypic diversity than allopatric sister taxa 67 . Populations that have multiple accumulated CRs can persist more easily in secondary sympatry than populations whose karyotypes are insufficiently divergent, thus leading sympatric populations to reproductive divergence 7,8,67 . On the other hand, the split of "M. morschi" lineages B and C would have involved both dispersion and vicariance events, allowing a possible colonization route of diversification to be traced.
The Pliocene-Pleistocene intervals had large-scale effects on habitat change, including the Brazilian coast, which was remarkably influenced by the sea level 35,58 . Sea level fluctuations related to climate change and geological disturbances would have dramatically increased over these periods to −100 m from current sea level 35,58 , thus allowed the opening of new previously unexploited regions by psammophilous taxa, including Mycetophylax ants. A direct consequence of this severe drop in sea level during the glacial periods would be the emergence of the Brazilian continental shelf 68 . This displacement of the southeastern South American coastline would have reached hundreds of kilometers to the east, enabling suitable climatic conditions for forest-adapted species 68 and presumably extending the geographic range of psammophilous taxa to areas previously flooded by the sea. Populations of "M. morschi" subjected to putative CRs would have been confronted with adaptively adverse conditions as a result of climatic and ecological stress during colonization, favoring CRs that would increase adaptation to these conditions whenever such conditions arose 10,11 . This pattern of chromosomal changes and local adaptation has been reported for Anopheles gambiae 11 and supports the diversity and karyotypic evolution in the psammophilous Mycetophylax clade. Additionally, the perspective of CRs fixation in isolated and/or peripheral founder populations through genetic drift and inbreeding is raised by many of the chromosomal speciation models [9][10][11] . In this way, the northernmost dispersion of the Brazilian Atlantic coast by the psammophilous Mycetophylax where molded by several historical events coupled with chromosomal changes, likely Rb fissions and fusions that shaped the karyotypes and lineage differentiation.
In conclusion, our integrative approach combining cytogenetic and molecular data has once again proved to be effective in supporting systematics and taxonomic research. We emphasize that the lack of cytogenetic data is still a hindrance to better application of these approaches. Further, we shed light in the cryptic diversity within "M. morschi" and, as previously suggested 17 , should be treated as phylogenetically isolated lineages or even as distinct species as they form well-resolved clades. The entire genus Mycetophylax has undergone systematic and taxonomic changes over the years, reinforcing the need for taxonomic revision. Psammophilous Mycetophylax represents a cytogenetically diverse group, showing great variation in terms of the number and morphology of the chromosomes and the localization of highly conserved sequences, such as the 45S rDNA cluster. Fossil-calibrated molecular dating suggests that the genus Mycetophylax probably originated at the Eocene-Oligocene boundary in the Amazon and/or Southern Grasslands, while psammophilous Mycetophylax emerged in the Middle Miocene in the South Atlantic coast, and most diversification occurred during subsequent major climatic events in the sand dune areas of the Brazilian Atlantic Forest biome. Putative CRs in the context of global climatic events and subsequent ecological stress could have somehow triggered and accompanied the diversification of the lineages in this highly dynamic environment.

Materials and Methods
Biological material. Populations of Mycetophylax were previously collected throughout their area of geographic distribution along the Brazilian Atlantic coast (Fig. 1) www.nature.com/scientificreports www.nature.com/scientificreports/ Chromosome preparation and FISH. Mitotic metaphases were obtained by the method described by Imai et al. 14 using brain ganglia of pre-pupae larvae through dissection in colchicine-hypotonic solution, following the adjusts described by Cardoso et al. 70 . Conventional giemsa staining was used to determine the chromosome number and morphology. We classified the chromosomes following the standard nomenclature for centromere position proposed by Levan et al. 71 with the small modification on the term acrocentric referring to all the chromosomes with an arm ratio greater than 7.0. The arm ratio (r) is the ratio of the length of the long arm of a chromosome to that of the short arm. Thus, metacentric chromosomes (M) are those with r = 1.0-1.7, submetacentrics (SM) have r = 1.7-3.0, subtelocentrics (ST) have r = 3.0-7.0, and acrocentrics have values from 7.0 to above.
DNA extraction, amplification, and sequencing. Genomic DNA extraction was performed from one ant worker per colony according to the standard CTAB/chloroform techniques 75 with adaptations for ant extraction. We sequenced fragments of five nuclear protein-coding genes: elongation factor 1-alpha-F1 (EF1α-F1), elongation factor 1-alpha-F2 (EF1α-F2), wingless (Wg), long-wavelength rhodopsin (LW Rh), and topoisomerase 1 (Top1) (GenBank accession numbers: MN745212-MN745286). The primers used to generate the sequence data are listed in Table S3. PCR was performed in a final volume of 25 μL which included 12.5 μL of GoTaq G2 Hot Start Colorless MasterMix (Taq DNA Polymerase, dNTPs, MgCl 2 , and buffer), 9.5 μL of Nuclease-Free Water, 1 μL of each pair of primers, and 1 μL of DNA, according to the manufacturer's instructions (PROMEGA, MADISON, WI, USA). Amplifications reaction conditions and sequencing were based on the methodology outlined in previous studies 20,76,77 . The chromatogram quality was verified by a Phred quality score higher than 30, and then the sequences were manually trimmed in GENEIOUS R7 78 .
Bayesian phylogenetic analysis. The concatenated dataset consisting of ~2.2 kbp (279 bp for EF1α-F1, 333 bp for Wg, 387 bp for LW Rh, 459 bp for EF1α-F2, and 669 bp for Top1) was added to the aligned dataset of Sosa-Calvo et al. 20 and aligned by using the MEGA7 79 . Data were partitioned and modeled using PARTITIONFINDER 2 80 under the Bayesian Information Criterion (BIC), greedy search scheme and with 15 input data blocks consisting of the first, second, and third codon positions of each of the five gene fragments. Thereafter, nine partitions were employed in independent Bayesian analyses using MRBAYES v3.2 81 . The standard nucleotide model was used for all partitions, and the nucleotide substitution models were variable for each partition (Table S4). Two simultaneous independent Monte Carlo Markov chain (MCMC) runs were performed with seven heated chains and one cold; the temperature was set at 0.05; and 10 million generations were performed with sampling every 1000 generations. The convergence among runs was verified by the average standard deviation of split frequencies that had to reach <0.01, and the appropriated burn-in was determined to be 10%. The tree-generated and posterior probabilities were visualized in FIGTREE v1.4 82 .

Divergence time estimates.
For divergence time analyses, we used BEAST v2.5 83 under the Fossilized Birth-Death (FBD) model 26 using an uncorrelated, log normal relaxed clock model to describe the branch-specific substitution rates 84 . The FBD model has become the most appropriate way to calibrate divergence time estimates when the calibration dates represented fossil occurrence times 26 . Accordingly, we integrated the occurrence times of 13 sets of Myrmicinae ant species fossils 85 into the tree prior to impose a time structure on the tree and to calibrate the analysis to absolute time. We calibrated internal nodes with minimum-age constraints using Myrmicinae amber fossils (Table S5). The root age was set at 98.6 Ma based on the estimated age for the subfamily Myrmicinae 77 . We set the nucleotide substitution model as GTR + I + G for all genes. Independent MCMC analyses were run for 60 million generations, sampling every 1000 generations. Runs were evaluated using TRACER v1.7 86 with effective sample size (ESS) values for all parameters over 200. The first 20% of the sampled tree topologies from the analyses were discarded as burn-in, and the remaining trees were summarized in TREEANNOTATOR v2.5. Before this, all fossils were removed from the tree using the FullToExtantTreeConverter tool (implemented in BEAUTI v2.5). The generated trees and credible intervals were visualized in both FIGTREE v1.4 82  www.nature.com/scientificreports www.nature.com/scientificreports/ Mycetarotes, Mycetosoritis, Kalathomyrmex and Cyatta species. It should be stressed that some specimens of the "strigatus group" were pruned from the tree as they are unknown lineages with unknown geographic ranges as well. We ran the phylogenetic tree using BEAST v2.5 83 with 25 million generations, sampling every 1000 generations, and a final burn-in of 10%. For this, we used a Yule speciation process, an uncorrelated, log normal relaxed clock model, and GTR + I + G for all genes. We characterized seven geographic ranges by the Neotropics that corresponded to the species distributions used: South Atlantic Coast (A); North Atlantic Coast (B); Brazilian Atlantic Forest (C); Northern Grasslands, including the Caatinga and Cerrado biomes (D); Southern Grasslands, including Chaco and Pampas biomes (E); Amazon Rainforest (F); and Middle America, covering all of Central America and Southern USA (G). The species distributions were taken from www.antmaps.org 88 .
We used the RASP v4.1 89 software package to estimate the likely ancestral ranges of Mycetophylax. RASP implements the major model-based approaches currently used, which enables comparisons to be made. Three approaches were employed for such estimates. The first one was carried out with the BioGeoBEARS (BGB) 90 to compare alternative biogeographic models. BGB can implement six specific models: the dispersal-extinction-cladogenesis model (DEC) 91 , a likelihood version of the dispersal-vicariance analysis (DIVALIKE), and the Bayesian inference of historical biogeography for discrete areas (BAYAREALIKE); and all three models can also be run with the founder-event speciation parameter (+j) 38,61 . The model also implements the rate of "dispersal" (range expansion) and "extinction" (range contraction) along the internal branches of the phylogeny 38,61 . In order to analyze the connectivity between the areas, we set up dispersion multipliers so that adjacent areas received a dispersion multiplier of 1 and non-adjacent areas received a value of 0.5. AIC results were used to assess the statistical significance of data-to-model fit.
The second approach used was BayArea 92 as it re-calculates the ancestral state probabilities of each unit area with alternative burn-in values 89 . BayArea also incorporates geographic coordinate data, making it an attractive feature for biogeographic investigations. However, BayArea may need to be repeated several times to get a stable result 89 , so we carried out independent runs with 10 million generations each, setting different burn-in values. Lastly, the third approach used was the Bayesian Binary MCMC (BBM) 93 . BBM analysis was performed using an estimated model, F81 + G, with a null root distribution, and the MCMC was run for 10 million generations with sampling every 1000 generations.
Ancestral chromosome number reconstruction. The chromosome evolution model and the reconstruction of ancestral chromosome numbers was recovered with CHROMEVOL v2.0 94,95 using both the maximum likelihood (ML) and Bayesian approaches. The chromosome evolution models evaluated were implemented as specified by Cardoso et al. 17 . To avoid over-representing the OTUs, and as several karyotype data on fungus-farming ants are unavailable, we used a restricted phylogeny, comprising data with cytogenetic information and only two species per genus per clade, since some genera were paraphyletic. The phylogeny was generated under the same specifications cited above for the MRBAYES analyses. The chromosome models over the phylogeny estimated in Bayesian inference and its null hypotheses were analyzed with 10,000 simulations, and the model that best fit the data set was selected under the Akaike Information Criterion (AIC).