Efficient Bayesian inference under the multispecies coalescent with migration.

Analyses of genome sequence data have revealed pervasive interspecific gene flow and enriched our understanding of the role of gene flow in speciation and adaptation. Inference of gene flow using genomic data requires powerful statistical methods. Yet current likelihood-based methods involve heavy computation and are feasible for small datasets only. Here, we implement the multispecies-coalescent-with-migration model in the Bayesian program bpp, which can be used to test for gene flow and estimate migration rates, as well as species divergence times and population sizes. We develop Markov chain Monte Carlo algorithms for efficient sampling from the posterior, enabling the analysis of genome-scale datasets with thousands of loci. Implementation of both introgression and migration models in the same program allows us to test whether gene flow occurred continuously over time or in pulses. Analyses of genomic data from Anopheles mosquitoes demonstrate rich information in typical genomic datasets about the mode and rate of gene flow.

Benefits and Limits of Phasing Alleles for Network Inference of Allopolyploid Complexes.

Accurately reconstructing the reticulate histories of polyploids remains a central challenge for understanding plant evolution. Although phylogenetic networks can provide insights into relationships among polyploid lineages, inferring networks may be hindered by the complexities of homology determination in polyploid taxa. We use simulations to show that phasing alleles from allopolyploid individuals can improve phylogenetic network inference under the multispecies coalescent by obtaining the true network with fewer loci compared to haplotype consensus sequences or sequences with heterozygous bases represented as ambiguity codes. Phased allelic data can also improve divergence time estimates for networks, which is helpful for evaluating allopolyploid speciation hypotheses and proposing mechanisms of speciation. To achieve these outcomes in empirical data, we present a novel pipeline that leverages a recently developed phasing algorithm to reliably phase alleles from polyploids. This pipeline is especially appropriate for target enrichment data, where depth of coverage is typically high enough to phase entire loci. We provide an empirical example in the North American Dryopteris fern complex that demonstrates insights from phased data as well as the challenges of network inference. We establish that our pipeline (PATÉ: Phased Alleles from Target Enrichment data) is capable of recovering a high proportion of phased loci from both diploids and polyploids. These data may improve network estimates compared to using haplotype consensus assemblies by accurately inferring the direction of gene flow, but statistical non-identifiability of phylogenetic networks poses a barrier to inferring the evolutionary history of reticulate complexes.

Likelihood-Based Tests of Species Tree Hypotheses.

Likelihood-based tests of phylogenetic trees are a foundation of modern systematics. Over the past decade, an enormous wealth and diversity of model-based approaches have been developed for phylogenetic inference of both gene trees and species trees. However, while many techniques exist for conducting formal likelihood-based tests of gene trees, such frameworks are comparatively underdeveloped and underutilized for testing species tree hypotheses. To date, widely used tests of tree topology are designed to assess the fit of classical models of molecular sequence data and individual gene trees and thus are not readily applicable to the problem of species tree inference. To address this issue, we derive several analogous likelihood-based approaches for testing topologies using modern species tree models and heuristic algorithms that use gene tree topologies as input for maximum likelihood estimation under the multispecies coalescent. For the purpose of comparing support for species trees, these tests leverage the statistical procedures of their original gene tree-based counterparts that have an extended history for testing phylogenetic hypotheses at a single locus. We discuss and demonstrate a number of applications, limitations, and important considerations of these tests using simulated and empirical phylogenomic data sets that include both bifurcating topologies and reticulate network models of species relationships. Finally, we introduce the open-source R package SpeciesTopoTestR (SpeciesTopology Tests in R) that includes a suite of functions for conducting formal likelihood-based tests of species topologies given a set of input gene tree topologies.

Inferring the Direction of Introgression Using Genomic Sequence Data.

Genomic data are informative about the history of species divergence and interspecific gene flow, including the direction, timing, and strength of gene flow. However, gene flow in opposite directions generates similar patterns in multilocus sequence data, such as reduced sequence divergence between the hybridizing species. As a result, inference of the direction of gene flow is challenging. Here, we investigate the information about the direction of gene flow present in genomic sequence data using likelihood-based methods under the multispecies-coalescent-with-introgression model. We analyze the case of two species, and use simulation to examine cases with three or four species. We find that it is easier to infer gene flow from a small population to a large one than in the opposite direction, and easier to infer inflow (gene flow from outgroup species to an ingroup species) than outflow (gene flow from an ingroup species to an outgroup species). It is also easier to infer gene flow if there is a longer time of separate evolution between the initial divergence and subsequent introgression. When introgression is assumed to occur in the wrong direction, the time of introgression tends to be correctly estimated and the Bayesian test of gene flow is often significant, while estimates of introgression probability can be even greater than the true probability. We analyze genomic sequences from Heliconius butterflies to demonstrate that typical genomic datasets are informative about the direction of interspecific gene flow, as well as its timing and strength.

Estimation of Cross-Species Introgression Rates Using Genomic Data Despite Model Unidentifiability.

Full-likelihood implementations of the multispecies coalescent with introgression (MSci) model treat genealogical fluctuations across the genome as a major source of information to infer the history of species divergence and gene flow using multilocus sequence data. However, MSci models are known to have unidentifiability issues, whereby different models or parameters make the same predictions about the data and cannot be distinguished by the data. Previous studies of unidentifiability have focused on heuristic methods based on gene trees and do not make an efficient use of the information in the data. Here we study the unidentifiability of MSci models under the full-likelihood methods. We characterize the unidentifiability of the bidirectional introgression (BDI) model, which assumes that gene flow occurs in both directions. We derive simple rules for arbitrary BDI models, which create unidentifiability of the label-switching type. In general, an MSci model with k BDI events has 2k unidentifiable modes or towers in the posterior, with each BDI event between sister species creating within-model parameter unidentifiability and each BDI event between nonsister species creating between-model unidentifiability. We develop novel algorithms for processing Markov chain Monte Carlo samples to remove label-switching problems and implement them in the bpp program. We analyze real and synthetic data to illustrate the utility of the BDI models and the new algorithms. We discuss the unidentifiability of heuristic methods and provide guidelines for the use of MSci models to infer gene flow using genomic data.

Bayesian Phylogenetic Inference using Relaxed-clocks and the Multispecies Coalescent.

The multispecies coalescent (MSC) model accommodates both species divergences and within-species coalescent and provides a natural framework for phylogenetic analysis of genomic data when the gene trees vary across the genome. The MSC model implemented in the program bpp assumes a molecular clock and the Jukes-Cantor model, and is suitable for analyzing genomic data from closely related species. Here we extend our implementation to more general substitution models and relaxed clocks to allow the rate to vary among species. The MSC-with-relaxed-clock model allows the estimation of species divergence times and ancestral population sizes using genomic sequences sampled from contemporary species when the strict clock assumption is violated, and provides a simulation framework for evaluating species tree estimation methods. We conducted simulations and analyzed two real datasets to evaluate the utility of the new models. We confirm that the clock-JC model is adequate for inference of shallow trees with closely related species, but it is important to account for clock violation for distant species. Our simulation suggests that there is valuable phylogenetic information in the gene-tree branch lengths even if the molecular clock assumption is seriously violated, and the relaxed-clock models implemented in bpp are able to extract such information. Our Markov chain Monte Carlo algorithms suffer from mixing problems when used for species tree estimation under the relaxed clock and we discuss possible improvements. We conclude that the new models are currently most effective for estimating population parameters such as species divergence times when the species tree is fixed.

Inference of Gene Flow between Species under Misspecified Models.

Genomic sequence data provide a rich source of information about the history of species divergence and interspecific hybridization or introgression. Despite recent advances in genomics and statistical methods, it remains challenging to infer gene flow, and as a result, one may have to estimate introgression rates and times under misspecified models. Here we use mathematical analysis and computer simulation to examine estimation bias and issues of interpretation when the model of gene flow is misspecified in analysis of genomic datasets, for example, if introgression is assigned to the wrong lineages. In the case of two species, we establish a correspondence between the migration rate in the continuous migration model and the introgression probability in the introgression model. When gene flow occurs continuously through time but in the analysis is assumed to occur at a fixed time point, common evolutionary parameters such as species divergence times are surprisingly well estimated. However, the time of introgression tends to be estimated towards the recent end of the period of continuous gene flow. When introgression events are assigned incorrectly to the parental or daughter lineages, introgression times tend to collapse onto species divergence times, with introgression probabilities underestimated. Overall, our analyses suggest that the simple introgression model is useful for extracting information concerning between-specific gene flow and divergence even when the model may be misspecified. However, for reliable inference of gene flow it is important to include multiple samples per species, in particular, from hybridizing species.

Molecular Clocks without Rocks: New Solutions for Old Problems.

Molecular data have been used to date species divergences ever since they were described as documents of evolutionary history in the 1960s. Yet, an inadequate fossil record and discordance between gene trees and species trees are persistently problematic. We examine how, by accommodating gene tree discordance and by scaling branch lengths to absolute time using mutation rate and generation time, multispecies coalescent (MSC) methods can potentially overcome these challenges. We find that time estimates can differ - in some cases, substantially - depending on whether MSC methods or traditional phylogenetic methods that apply concatenation are used, and whether the tree is calibrated with pedigree-based mutation rates or with fossils. We discuss the advantages and shortcomings of both approaches and provide practical guidance for data analysis when using these methods.

Gene flow assessment helps to distinguish strong genomic structure from speciation in an Iberian ant-eating spider.

Although genomic data is boosting our understanding of evolution, we still lack a solid framework to perform reliable genome-based species delineation. This problem is especially critical in the case of phylogeographically structured organisms, with allopatric populations showing similar divergence patterns as species. Here, we assess the species limits and phylogeography of Zodarion alacre, an ant-eating spider widely distributed across the Iberian Peninsula. We first performed species delimitation based on genome-wide data and then validated these results using additional evidence. A commonly employed species delimitation strategy detected four distinct lineages with almost no admixture, which present allopatric distributions. These lineages showed ecological differentiation but no clear morphological differentiation, and evidence of introgression in a mitochondrial barcode. Phylogenomic networks found evidence of substantial gene flow between lineages. Finally, phylogeographic methods highlighted remarkable isolation by distance and detected evidence of range expansion from south-central Portugal to central-north Spain. We conclude that despite their deep genomic differentiation, the lineages of Z. alacre do not show evidence of complete speciation. Our results likely shed light on why Zodarion is among the most diversified spider genera despite its limited distribution and support the use of gene flow evidence to inform species boundaries.

Comparing inference under the multispecies coalescent with and without recombination.

Accurate inference of population parameters plays a pivotal role in unravelling evolutionary histories. While recombination has been universally accepted as a fundamental process in the evolution of sexually reproducing organisms, it remains challenging to model it exactly. Thus, existing coalescent-based approaches make different assumptions or approximations to facilitate phylogenetic inference, which can potentially bring about biases in estimates of evolutionary parameters when recombination is present. In this article, we evaluate the performance of population parameter estimation using three methods-StarBEAST2, SNAPP, and diCal2-that represent three different types of inference. We performed whole-genome simulations in which recombination rates, mutation rates, and levels of incomplete lineage sorting were varied. We show that StarBEAST2 using short or medium-sized loci is robust to realistic rates of recombination, which is in agreement with previous studies. SNAPP, as expected, is generally unaffected by recombination events. Most surprisingly, diCal2, a method that is designed to explicitly account for recombination, performs considerably worse than other methods under comparison.

Interspecific gene flow obscures phylogenetic relationships in an important insect pest species complex.

As genomic data proliferates, the prevalence of post-speciation gene flow is making species boundaries and relationships increasingly ambiguous. Although current approaches inferring fully bifurcating phylogenies based on concatenated datasets provide simple and robust answers to many species relationships, they may be inaccurate because the models ignore inter-specific gene flow and incomplete lineage sorting. To examine the potential error resulting from ignoring gene flow, we generated both a RAD-seq and a 500 protein-coding loci highly multiplexed amplicon (HiMAP) dataset for a monophyletic group of 12 species defined as the Bactrocera dorsalis sensu lato clade. With some of the world's worst agricultural pests, the taxonomy of the B. dorsalis s.l. clade is important for trade and quarantines. However, taxonomic confusion confounds resolution due to intra- and interspecific phenotypic variation and convergence, mitochondrial introgression across half of the species, and viable hybrids. We compared the topological convergence of our datasets using concatenated phylogenetic and various multispecies coalescent approaches, some of which account for gene flow. All analyses agreed on species delimitation, but there was incongruence between species relationships. Under concatenation, both datasets suggest identical species relationships with mostly high statistical support. However, multispecies coalescent and multispecies network approaches suggest markedly different hypotheses and detected significant gene flow. We suggest that the network approaches are likely more accurate because gene flow violates the assumptions of the concatenated phylogenetic analyses, but the data-reductive requirements of network approaches resulted in reduced statistical support and could not unambiguously resolve gene flow directions. Our study highlights the importance of testing for gene flow, particularly with phylogenomic datasets, even when concatenated approaches receive high statistical support.

Complexity of the simplest species tree problem.

The multispecies coalescent model provides a natural framework for species tree estimation accounting for gene-tree conflicts. Although a number of species tree methods under the multispecies coalescent have been suggested and evaluated using simulation, their statistical properties remain poorly understood. Here, we use mathematical analysis aided by computer simulation to examine the identifiability, consistency, and efficiency of different species tree methods in the case of three species and three sequences under the molecular clock. We consider four major species-tree methods including concatenation, two-step, independent-sites maximum likelihood, and maximum likelihood. We develop approximations that predict that the probit transform of the species tree estimation error decreases linearly with the square root of the number of loci. Even in this simplest case, major differences exist among the methods. Full-likelihood methods are considerably more efficient than summary methods such as concatenation and two-step. They also provide estimates of important parameters such as species divergence times and ancestral population sizes,whereas these parameters are not identifiable by summary methods. Our results highlight the need to improve the statistical efficiency of summary methods and the computational efficiency of full likelihood methods of species tree estimation.

A simulation study to examine the impact of recombination on phylogenomic inferences under the multispecies coalescent model.

Phylogenomic analyses under the multispecies coalescent model assume no recombination within locus and free recombination among loci. Yet, in real data sets intralocus recombination causes different sites of the same locus to have different genealogical histories so that the model is misspecified. The impact of recombination on various coalescent-based phylogenomic analyses has not been systematically examined. Here, we conduct a computer simulation to examine the impact of recombination on several Bayesian analyses of multilocus sequence data, including species tree estimation, species delimitation (by Bayesian selection of delimitation models) and estimation of evolutionary parameters such as species divergence and introgression times, population sizes for modern and extinct species, and cross-species introgression probabilities. We found that recombination, at rates comparable to estimates from the human being, has little impact on coalescent-based species tree estimation, species delimitation and estimation of population parameters. At rates 10 times higher than the human rate, recombination may affect parameter estimation, causing positive biases in introgression times and ancestral population sizes, although species divergence times and cross-species introgression probabilities are estimated with little bias. Overall, the simulation suggests that phylogenomic inferences under the multispecies coalescent model are robust to realistic amounts of intralocus recombination.

Phylogenetics in space: How continuous spatial structure impacts tree inference.

The tendency to discretize biology permeates taxonomy and systematics, leading to models that simplify the often continuous nature of populations. Even when the assumption of panmixia is relaxed, most models still assume some degree of discrete structure. The multispecies coalescent has emerged as a powerful model in phylogenetics, but in its common implementation is entirely space-independent - what we call the "missing z-axis". In this article, we review the many lines of evidence for how continuous spatial structure can impact phylogenetic inference. We illustrate and expand on these by using complex continuous-space demographic models that include distinct modes of speciation. We find that the impact of spatial structure permeates all aspects of phylogenetic inference, including gene tree stoichiometry, topological and branch-length variance, network estimation, and species delimitation. We conclude by utilizing our results to suggest how researchers can identify spatial structure in phylogenetic datasets.

Gauging ages of tiger swallowtail butterflies using alternate SNP analyses.

Divergence times underpin diverse evolutionary hypotheses, but conflicting age estimates across studies diminish the validity of such hypotheses. These conflicts have continued to grow as large genomics datasets become commonplace and analytical approaches proliferate. To provide more stable temporal intervals, age estimations should be interpreted in the context of both the type of data and analysis being used. Here, we use multispecies coalescent (MSC), concatenation-based, and categorical data transformation approaches on genome-wide SNP data to infer divergence ages within the Papilio glaucus group of tiger swallowtail butterflies in North America. While the SNP data supported previously recognized relationships within the group (P. multicaudata, ((P. eurymedon, P. rutulus), (P. appalachiensis, P. canadensis, P. glaucus))), estimated ages of divergence between the major lineages varied substantially among analyses. MSC produced wide credibility intervals particularly for deeper nodes, reflecting uncertainty in the coalescence times as a possible result of conflicting signal across gene trees. Concatenation, in contrast, gave narrower and more well-defined posterior distributions for the node ages; however, the higher precision of these time estimates is a likely artefact due to more simplistic underlying assumptions of this approach that do not account for conflict among gene trees. Transformed categorical data analysis gave the least precise and the most variable results, with its simple substitution model coupled with a relaxed clock tending to produce spurious results from large genome-wide datasets. While median node ages differed considerably between analyses (∼2 Mya between MSC and concatenation-based results), their corresponding credibility intervals nonetheless highlight common temporal patterns for deeper divergences in the group as well as finer-scale phylogeography. Age distributions across analyses support an origin of the group during the warm period of the early to mid-Pliocene. Late Pliocene climate aridification and cooling drove divergence between eastern and western groups that further diversified during the period of repeated Pleistocene glaciations. Our results provide a structured comparative assessment of divergence time estimates and evolutionary relationships in a well-studied group of butterflies, and support better understanding of analytical biases in divergence time estimation.

Genome-wide ultraconserved elements resolve phylogenetic relationships and biogeographic history among Neotropical leaf-nosed bats in the genus Anoura (Phyllostomidae).

AnouraGray, 1838 are Neotropical nectarivorous bats and the most speciose genus within the phyllostomid subfamily Glossophaginae. However, Anoura species limits remain debated, and phylogenetic relationships remain poorly known, because previous studies used limited Anoura taxon sampling or focused primarily on higher-level relationships. Here, we conduct the first phylogenomic study of Anoura by analyzing 2039 genome-wide ultraconserved elements (UCEs) sequenced for 42 individuals from 8 Anoura species/lineages plus two outgroups. Overall, our results based on UCEs resolved relationships in the genus and supported (1) the monophyly of small-bodied Anoura species (previously genus Lonchoglossa); (2) monotypic status of A. caudifer; and (3) nested positions of "A. carishina", A. caudifer aequatoris, and A. geoffroyi peruana specimens within A. latidens, A. caudifer and A. geoffroyi, respectively (suggesting that these taxa are not distinct species). Additionally, (4) phylogenetic networks allowing reticulate edges did not explain gene tree discordance better than the species tree (without introgression), indicating that a coalescent model accounting for discordance solely through incomplete lineage sorting fit our data well. Sensitivity analyses indicated that our species tree results were not adversely affected by varying taxon sampling across loci. Tree calibration and Bayesian coalescent analyses dated the onset of diversification within Anoura to around âˆ¼ 6-9 million years ago in the Miocene, with extant species diverging mainly within the past âˆ¼ 4 million years. We inferred a historical biogeographical scenario for Anoura of parapatric speciation fragmenting the range of a wide-ranging ancestral lineage centered in the Central to Northern Andes, along with Pliocene-Pleistocene dispersal or founder event speciation in Amazonia and the Brazilian Atlantic forest during the last âˆ¼ 2.5 million years.

Using natural history to guide supervised machine learning for cryptic species delimitation with genetic data.

The diversity of biological and ecological characteristics of organisms, and the underlying genetic patterns and processes of speciation, makes the development of universally applicable genetic species delimitation methods challenging. Many approaches, like those incorporating the multispecies coalescent, sometimes delimit populations and overestimate species numbers. This issue is exacerbated in taxa with inherently high population structure due to low dispersal ability, and in cryptic species resulting from nonecological speciation. These taxa present a conundrum when delimiting species: analyses rely heavily, if not entirely, on genetic data which over split species, while other lines of evidence lump. We showcase this conundrum in the harvester Theromaster brunneus, a low dispersal taxon with a wide geographic distribution and high potential for cryptic species. Integrating morphology, mitochondrial, and sub-genomic (double-digest RADSeq and ultraconserved elements) data, we find high discordance across analyses and data types in the number of inferred species, with further evidence that multispecies coalescent approaches over split. We demonstrate the power of a supervised machine learning approach in effectively delimiting cryptic species by creating a "custom" training data set derived from a well-studied lineage with similar biological characteristics as Theromaster. This novel approach uses known taxa with particular biological characteristics to inform unknown taxa with similar characteristics, using modern computational tools ideally suited for species delimitation. The approach also considers the natural history of organisms to make more biologically informed species delimitation decisions, and in principle is broadly applicable for taxa across the tree of life.

Site pattern probabilities under the multispecies coalescent and a relaxed molecular clock: Theory and applications.

The first step in statistical inference of the evolutionary histories of species is developing a probability model that describes the mutation process as accurately and realistically as possible. A major complication of this inference is that different loci on the genome can have histories that diverge from the common species history and each other. The multispecies coalescent process is commonly used to model one source of this divergence, incomplete lineage sorting, or ILS. In Chifman and Kubatko (2015), the authors computed the site pattern probabilities for four taxa under a full probability model based on the Jukes-Cantor substitution model when the molecular clock holds. This paper generalizes that work to a relaxed clock model, allowing for mutation rates to differ among species. This will enable better phylogentic inference in cases where the molecular clock does not hold.

A monograph of Aspergillus section Candidi.

Aspergillus section Candidi encompasses white- or yellow-sporulating species mostly isolated from indoor and cave environments, food, feed, clinical material, soil and dung. Their identification is non-trivial due to largely uniform morphology. This study aims to re-evaluate the species boundaries in the section Candidi and present an overview of all existing species along with information on their ecology. For the analyses, we assembled a set of 113 strains with diverse origin. For the molecular analyses, we used DNA sequences of three house-keeping genes (benA, CaM and RPB2) and employed species delimitation methods based on a multispecies coalescent model. Classical phylogenetic methods and genealogical concordance phylogenetic species recognition (GCPSR) approaches were used for comparison. Phenotypic studies involved comparisons of macromorphology on four cultivation media, seven micromorphological characters and growth at temperatures ranging from 10 to 45 °C. Based on the integrative approach comprising four criteria (phylogenetic and phenotypic), all currently accepted species gained support, while two new species are proposed (A. magnus and A. tenebricus). In addition, we proposed the new name A. neotritici to replace an invalidly described A. tritici. The revised section Candidi now encompasses nine species, some of which manifest a high level of intraspecific genetic and/or phenotypic variability (e.g., A. subalbidus and A. campestris) while others are more uniform (e.g., A. candidus or A. pragensis). The growth rates on different media and at different temperatures, colony colours, production of soluble pigments, stipe dimensions and vesicle diameters contributed the most to the phenotypic species differentiation. Taxonomic novelties: New species: Aspergillus magnus Glässnerová & Hubka; Aspergillus neotritici Glässnerová & Hubka; Aspergillus tenebricus Houbraken, Glässnerová & Hubka. Citation: Glässnerová K, Sklenár F, Jurjevic Z, Houbraken J, Yaguchi T, Visagie CM, Gené J, Siqueira JPZ, Kubátová A, Kolarík M, Hubka V (2022). A monograph of Aspergillus section Candidi. Studies in Mycology 102: 1-51. doi: 10.3114/sim.2022.102.01.

Taxonomy of Aspergillus series Versicolores: species reduction and lessons learned about intraspecific variability.

Aspergillus series Versicolores members occur in a wide range of environments and substrates such as indoor environments, food, clinical materials, soil, caves, marine or hypersaline ecosystems. The taxonomy of the series has undergone numerous re-arrangements including a drastic reduction in the number of species and subsequent recovery to 17 species in the last decade. The identification to species level is however problematic or impossible in some isolates even using DNA sequencing or MALDI-TOF mass spectrometry indicating a problem in the definition of species boundaries. To revise the species limits, we assembled a large dataset of 518 strains. From these, a total of 213 strains were selected for the final analysis according to their calmodulin (CaM) genotype, substrate and geography. This set was used for phylogenetic analysis based on five loci (benA, CaM, RPB2, Mcm7, Tsr1). Apart from the classical phylogenetic methods, we used multispecies coalescence (MSC) model-based methods, including one multilocus method (STACEY) and five single-locus methods (GMYC, bGMYC, PTP, bPTP, ABGD). Almost all species delimitation methods suggested a broad species concept with only four species consistently supported. We also demonstrated that the currently applied concept of species is not sustainable as there are incongruences between single-gene phylogenies resulting in different species identifications when using different gene regions. Morphological and physiological data showed overall lack of good, taxonomically informative characters, which could be used for identification of such a large number of existing species. The characters expressed either low variability across species or significant intraspecific variability exceeding interspecific variability. Based on the above-mentioned results, we reduce series Versicolores to four species, namely A. versicolor, A. creber, A. sydowii and A. subversicolor, and the remaining species are synonymized with either A. versicolor or A. creber. The revised descriptions of the four accepted species are provided. They can all be identified by any of the five genes used in this study. Despite the large reduction in species number, identification based on phenotypic characters remains challenging, because the variation in phenotypic characters is high and overlapping among species, especially between A. versicolor and A. creber. Similar to the 17 narrowly defined species, the four broadly defined species do not have a specific ecology and are distributed worldwide. We expect that the application of comparable methodology with extensive sampling could lead to a similar reduction in the number of cryptic species in other extensively studied Aspergillus species complexes and other fungal genera. Citation: Sklenár F, Glässnerová K, Jurjevic Z, Houbraken J, Samson RA, Visagie CM, Yilmaz N, Gené J, Cano J, Chen AJ, Nováková A, Yaguchi T, Kolarík M, Hubka V (2022). Taxonomy of Aspergillus series Versicolores: species reduction and lessons learned about intraspecific variability. Studies in Mycology 102 : 53-93. doi: 10.3114/sim.2022.102.02.