Revisiting the multispecies coalescent model fit with an example from a complete molecular phylogeny of the Liolaemus wiegmannii species group (Squamata: Liolaemidae)

Departures from the Multispecies Coalescent (MSC) assumptions could cause artefactual topologies and node heights estimates, and therefore, trees inferred without MSC model fit testing could potentially misrepresent an accurate approximation of the evolutionary history of a group. The current implementation of the MSC model testing for non-genomic level molecular markers cannot process trees estimated from BEAST 2, limiting its application for large datasets of sequence-based markers. Here we recode functions of the R package P2C2M to assess model fit to the MSC and apply this new implementation, which we named P2C2M2, to test the MSC model in a 16-loci dataset of 42 lizard species focused on the Liolaemus wiegmannii group. We found strong evidence of model departures in several loci, possibly due to historical gene flow, which could also be causing an unexpected position of the L. wiegmannii group within the L. montanus section of Eulaemus, when hybridization is not accounted for. The L. anomalus group is inferred as the closest to the L. wiegmannii group when gene flow is incorporated via a Multispecies Network Coalescent model, and a reticulation, suggesting historical gene flow between the L. wiegmannii and L. montanus groups is inferred, which has not been previously reported. We argue that there are at least three sources of discrepancy between the literature and the node ages estimated in our study: the use of strict molecular clocks without statistical justification, misplaced fossil calibrations, and the estimation of coalescent times instead of species divergence times. We encouraged systematists to routinely test the fit of the MSC model when estimating species trees using sequence-based markers, and to follow a phylogenetic network approach when both this test is significant and when historical gene flow is considered one plausible source of the departure from the MSC model.

Methods
Details on data collection and analysis can be found in the Material and Methods section of the article: Revisiting the problem of Multispecies Coalescent Model fit with an example from a complete molecular phylogeny of the Liolaemus wiegmannii species group (Squamata: Liolaemidae)

In general, the primary data available here includes DNA sequence data for 16 loci (13 nuclear loci and 3 mitochondrial) generated by Sanger sequencing. Species sampled for this dataset were focused on members of the Liolaemus montanus section of the Eulaemus subgenus of Liolaemus (Squamata: Liolaemidae), including all the known species/lineages of the Liolaemus wiegmannii group. These alignments comprise both published data from GenBank and originally generated sequences for this study.

Several XML files for replicating BEAST 2 runs are also included. XML files for infer input trees to run P2C2M2, and XML for marginal likelihood estimates were manually edited. Instructions on these editions can be found at https://github.com/arleyc/P2C2M2 and https://www.beast2.org/path-sampling/, respectively.

Phylogenetic networks inference were made using a maximum pseudolikelihood approach implemented in PhyloNet 3.8.2. Networks were inferred based on Bayesian gene trees, allowing the number of reticulations (h) from 0 to 6. The nexus control file supplied in this dataset will replicate a 10 runs inference for h=2, which was the number of reticulations chosen based on a loglikelihood criteria. Alternative networks for h values different than 2 can be replicated by replacing this parameter value from the control file provided here.