2.3.1. Single-Locus Tree-Based Methods

We predicted species limits utilizing the Generalized Mixed Yule Coalescent (GMYC) approach as described by Fujisawa and Barraclough [63] using the R package splits (default settings, including single threshold). We first generated an ultrametric tree in BEAST v1.10.4 [64–66], using the following settings according to Lücking et al. [32]: (a) strict

molecular clock; (b) GTR substitution model with base frequencies estimated and Gamma and invariant sites with six Gamma categories; (c) speciation through a Yule process with the "yule.birthRate" prior set to an exponential distribution with 4.0 as mean; (d) the "ucld.mean" prior set to an exponential distribution with 0.001 as mean and all other priors with default values; and (e) 100,000,000 MCMC generations.

We also inferred putative species boundaries applying the Poisson Tree Process (PTP), including its Bayesian implementation (bPTP) [67] on the bPTP Web Server (https://species. h-its.org/ptp/ accessed on 15 January 2022), with the following settings: 500,000 MCMC generations, 100 Thinning, and 0.25 Burn-in.

### 2.3.2. Distance-Based Methods and Barcoding Gap Analyses

For the barcoding gap analysis, we used a subset of the ITS sequence data limited to the genus *Cora*, retaining only terminals with less than 30% missing data, spanning most or all of the ITS1 and ITS2 region (716 terminals; File S4). This subset corresponded to 175 species-level taxa plus one subspecies delimited ad hoc (see above). We computed a distance matrix using the Kimura 2-parameter model implemented in BioEdit 7 through DNADIST 3.5 [68,69]. Within-species and between-species distances were then assessed by arranging the terminals according to the ITS-based ML phylogeny (see above) and comparing the distances between subsequent terminals, using our ad hoc species hypotheses as grouping variables. We computed a threshold distance value by comparing within-species and between-species distances, selecting the value that best discriminated between the two, retaining a maximum of within-species and a minimum of between-species pairs. The resulting threshold was then used to visualize distance patterns in the full two-dimensional distance matrix between all terminals, enabling the assessment of the predefined specieslevel clades into three categories: (1) species-level clades well-delimited by the threshold distance; (2) species-level clades with internal variation frequently greater than the threshold distance; and (3) species-complexes in which hypothesized species-level clades were not well-resolved by the threshold distance.

We further employed Automatic Barcode Gap Discovery (ABGD) [70] and Assemble Species by Automatic Partitioning (ASAP) [71] on the same dataset of 716 *Cora* terminals (File S4), using the corresponding web servers [https://bioinfo.mnhn.fr/abi/ public/abgd/abgdweb.html; https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html accessed on 15 January 2022]. For ABDG, we used the following settings: Kimura (K80) 2-parameter model, Pmin = 0.001, Pmax = 0.1, Steps = 10, X (relative gap width) = 1.5, 1.0. 0.5, 0.1 (stepwise), Nbins = 20. For ASAP, we used settings as follows: Kimura (K80) 2-parameter model, Split groups below this probability = 0.01, Keep 10 best scores, Use fixed seed value = −1, Highlights results between the genetic distances = 0.05 and 0.05.
