Phylogenetic analyses and ancestral range estimations

Phylogenetic relationships among the Jassa species were reconstructed using only unique haplotypes with maximum likelihood (ML) and Bayesian inference (BI) of the CO1 gene fragment. The ML was inferred with PhyML (Guindon et al., 2010) using 1,000 bootstraps for the branch support. The best substitution model (GTR + I + G) was tested with the SMS routine in PhyML using both AIC and BIC as optimality criteria (Lefort, Longueville & Gascuel, 2017).

A time-calibrated BI phylogeny was reconstructed in BEAST 2.5.2 (Bouckaert et al., 2014) on XSEDE (Towns et al., 2014). The GTR model of evolution, with proportion of invariant (I) and gamma shape parameters (G), determined with bModelTest (Bouckaert & Drummond, 2017), and the Yule speciation model were set for priors. An uncorrelated relaxed clock with log-normal distribution was applied following the specifications proposed by Copilaş-Ciocianu, Sidorov & Gontcharov (2019); that is, rates ranged from 0.7% to 1.7% Ma⁻¹ with a starting value of 1.2%. Three runs each of 20,000,000 iterations of MCMC sampled each 1,000 iterations were performed. All runs were examined using Tracer v1.7.1 and all sampled parameters achieved sufficient sample sizes (ESS > 200). Tree files were combined using LogCombiner on XSEDE (1.8.4) with 15% of burn-in. The maximum clade credibility tree was generated using TreeAnnotator on XSEDE. All the XSEDE analyses were performed on the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2010).

The ancestral ranges of the haplotypes were estimated with the R package BioGeoBears (Matzke, 2013; R Development Core Team, 2018). This package is used to perform biogeographic inferences such as the estimation of ancestral ranges of species or Operational Taxonomic Units (OTUs) under different hypotheses (e.g., dispersion, founder-event speciation). Furthermore, BioGeoBears allows for the implementation of a third parameter (J) which permits a “jump speciation” in the daughter lineage, resulting in a possible different area from the direct ancestor. This feature accounts for the biology of Jassa species in this study by virtually including anthropogenic dispersal. The analyses were run using the DEC and DEC+J models (Matzke, 2014) and a maximum of five (J. marmorata) and four (J. slatteryi) areas of possible occurence as this was the maximum number of regions where a single haplotype was found for each species respectively. In order to avoid overinterpretation, only the nodes with at least 70% of bootstrap support and 90% of posterior probabilities were considered. The best-fitting model was selected based on Akaike’s information criterion (AIC).