Phylogeny and divergence time estimation of Io moths and relatives (Lepidoptera: Saturniidae: Automeris)

Taxon sampling and extraction

We sampled and extracted 115 species and subspecies of Automeris (98), Leucanella (five), Pseudautomeris (six), and Molippa (six) available in the McGuire Center for Lepidoptera and Biodiversity (MGCL), Florida Museum of Natural History (FLMNH). Specimens were obtained from two collection types: dry, pinned specimens and those stored in ethanol at −80 °C, which have been wing-vouchered following Cho et al. (2016) and are specifically stored for use in molecular studies.

One specimen per taxon was selected from the collections. Pinned specimens were carefully selected based on their condition, focusing on those with intact wings and body parts. Priority was given to more recently collected specimens during the selection process to enhance the likelihood of successful DNA extraction. Two legs were pulled from each sample and placed in a 1.5 mL microcentrifuge tube. Corresponding identification labels for samples were recorded on the respective tubes, linking them to those collection data present in the MGCL specimen database. Similarly, molecular cold storage samples were taken from their respective boxes, and partial thoracic tissue was dissected and placed in a labeled 1.5 mL microcentrifuge tube, and all relevant information was recorded in the collective database. The extracted thorax tissue was then placed in a labeled 1.5 mL microcentrifuge tube, and all relevant information was recorded in a collective database. DNA extraction of legs from the MGCL pinned collection or thorax tissue from the molecular collection was chosen based on availability of species within the MGCL pinned and molecular collection. Due to the likelihood of higher DNA yield in thorax tissue, specimens from the molecular collection were given priority. If specimens were available in both molecular and pinned collections, specimens were chosen from the molecular collection.

Extractions were completed following the specified protocol of the Qiagen DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany). For quality control, extracted DNA was checked for concentration and fragmentation with a Qubit 2.0 fluorometer and electrophoresis gels (Fisherbrand Electrophoresis Power Supplies, FB200; Thermo Fisher Scientific, Waltham, MA, USA). After DNA extraction, if the desired concentration of DNA was not achieved, samples were re-extracted if tissue was still available. Otherwise, species with multiple samples were re-extracted with a different individual. Once the quantity of DNA was determined (8 ng/uL minimum and 100 ng/uL maximum), extracted samples were sent for Anchored Hybrid Enrichment (AHE) library preparation at RAPiD Genomics in Gainesville, Florida, USA. Higher concentration samples were prioritized. Out of the 115 extracted samples, 113 showed DNA concentration of ≥8 ng/uL. Library preparation, hybridization enrichment, and Illumina HiSeq 2500 sequencing (PE100) was carried out at RAPiD Genomics.

DNA sequencing and assembly

We used the “BOM1” Anchored Hybrid Enrichment probe set (Hamilton et al., 2019), which was developed to target 921 loci across Bombycoidea. We used this probe set because it has been proven effective in capturing sequence data from both dry-pinned and ethanol-preserved specimens (Hamilton et al., 2019; Dowdy et al., 2020; Li et al., 2022) This probe kit includes 58 loci across 24 “legacy” Sanger sequencing-based genes (Regier et al., 2008), eight bombycoid vision-related loci, and 855 loci designated as the Lepidoptera Agilent Custom SureSelect Target Enrichment “LEP1” probe kit from Breinholt et al. (2018). Anchored Hybrid Enrichment (AHE) is a sequencing technique specifically developed to target and capture a large number of orthologous loci from the genome. This method is well-suited for resolving evolutionary relationships, both at deep and shallow levels. The probes used in this technique bind to conserved regions flanked by variable regions that are distributed randomly throughout the genome (Hamilton et al., 2019). This method generates a varied and informative set of loci, containing exons, introns, intergenic, and conserved regions of the genome (Lemmon, Emme & Lemmon, 2012).

Raw AHE sequences were assembled following the methods of Breinholt et al. (2018), which implements an Iterative Baited Assembly (IBA) approach. This involves using the original sequencing probes to identify matching raw reads, which are then assembled into novel probes. The newly assembled probes are subsequently used as a query to match against the remaining raw reads, and the process is repeated iteratively until confident assemblies can no longer be obtained. The pipeline also checks for quality and cross contamination due to barcode leakage and removes paralogs. Resultant assemblies extend beyond the boundaries of the initial sequencing probes, thereby leading to the production of two distinct datasets—one comprising sequences solely from the probe region of the assembly (Probe dataset), and the other comprising sequences from the complete assembly, encompassing both the probe and outer flanking regions (PF dataset).

Because the IBA approach often resulted in multiple assembled sequences for each locus per specimen, sequences were aligned using MAFFT v7.245 (Katoh & Standley, 2013) and a 50% consensus generated using FASconCAT-G v1.02 (Kück & Longo, 2014) with the ‘-c -c -c’ command to result in one sequence per locus per specimen. To minimize the extent of missing data in the final concatenated dataset, loci that were only obtained from three or fewer species were excluded from the datasets. Loci were concatenated across all species into one supermatrix using FASconCAT-G with the ‘-s’ command.

A total of 113 specimens were successfully sequenced with AHE. Four outgroup species were chosen to provide secondary calibrations for the divergence time analysis and to provide a root for the tree. All outgroups chosen were genera included in the Hamilton et al. (2019) phylogeny and all had available raw transcriptomes or genomes. We chose two outgroup taxa in Saturniidae, Attacus atlas and Therinia lactucina, one from Sphingidae, Manduca sexta, and one from Bombycidae, Bombyx mori (The International Silkworm Genome Consortium, 2008; Breinholt & Kawahara, 2013; Kawahara & Breinholt, 2014; Kanost et al., 2016). The two transcriptomes (Saturniidae) and two genomes (Bombyx and Manduca) that we used for the outgroups were assembled to the BOM1 AHE probe regions using the methods described above. This enabled us to combine the transcriptome and genome datasets with our newly sequenced data that were also assembled to the BOM1 AHE probe regions.

Alistat v1.6 (Wong et al., 2020) was used to evaluate the completeness of the concatenated alignment for Probe and PF datasets. We created Probe only and PF datasets because PF datasets can sometimes yield greater robustness to phylogenetic analyses of AHE data (e.g., Kawahara et al., 2018; Hamilton et al., 2019; Homziak et al., 2019; St Laurent et al., 2021). Completeness scores (C-scores) were computed for each taxon (Cr), and those taxa with a C-score <0.15 were removed to avoid specimens with poor capture quality. We also removed loci that were captured across ≤ 3 taxa to improve dataset quality.

Phylogenetic analysis

Phylogenetic inference was conducted on concatenated supermatrices using a maximum likelihood analysis with IQTREE v. 1.5 (Nguyen et al., 2015). We constructed two datasets, one with just the probe region, and another with both probe and PF regions. Both datasets were analyzed as nucleotides, and we determined the best substitution model and partitioning scheme using the command ‘-m MFP+MERGE’ using ModelFinder (Lanfear et al., 2012) as implemented in IQTREE. The command ‘-B 1000 -bnni’ was used to perform 1000 ultrafast bootstrap (UFBS) replicates, while optimizing each bootstrap tree using a hill-climbing nearest neighbor interchange (NNI) search to reduce the risk of overestimating branch supports. All trees were rooted to Bombyx mori. We refer to high node support as those with UFBS ≥95.

Divergence time estimation

We performed divergence time estimation in BEAST v2.6.7 (Bouckaert et al., 2019) using the topology generated by IQTREE from the PF dataset as the starting tree. This specific topology was used since it yielded higher bootstrap values than those obtained using the Probe dataset. However, when running BEAST using the PF dataset for dating, the resulting 95% confidence intervals were disproportionately small, indicating the high likelihood that the large dataset caused the underlying Bayesian analysis to become easily stuck at a local minima (Kawahara et al., 2019; Rougerie et al., 2022). Due to this, we ultimately chose to use the Probe dataset for dating.

Loci were partitioned based on the best partitioning scheme as previously identified by ModelFinder in IQTREE and their corresponding site models (herein all set to the HKY substitution model) unlinked. For this analysis we opted to link the clock model (a relaxed clock with a log normal distribution) across all partitions because analyses running the dataset with unlinked clock models failed to converge particularly with regard to the resulting estimated divergence times. We applied a mean clock rate of 0.41 substitutions per site per 100 million years with “Mean In Real Space” checked based on the mutation rate that was recently estimated for Bombyx mori by Han et al. (2023) and set the ‘S’ parameter of our clock model prior to 0.1 so that the log normal distribution would closely match their reported 95% confidence interval of 0.33 × 10−8–0.49 × 10−8 per site per generation (we treat a generation as one year).

Among insects, there are disproportionately few Lepidoptera fossils (Labandeira & Sepkoski, 1993). Therefore, our tree was calibrated using ranges of dates obtained from Kawahara et al. (2019). We used four secondary calibration points with normal distributions to constrain the most recent common ancestors (MRCAs) of Bombycidae + (Sphingidae + Saturniidae) (the root node), Sphingidae + Saturniidae (node 7), Oxyteninae + remaining Saturniidae (node 6), and the MRCA of Saturniinae + Hemileucinae (node 5) (see Table 1). We used the Calibrated Yule model of speciation (Yule, 1924) as the tree prior, and ran two independent MCMC chains for 40 million generations each, sampling every 1,000 generations. TRACER v1.6.0 (Rambaut & Drummond, 2013s; Rambaut & Drummond, 2013b) was used to assess stationarity of and convergence between runs. Trees were then combined across runs with LogCombiner v2.6.3 after discarding a conservative 50% burn-in and a maximum clade credibility tree with median heights was recovered with TreeAnnotator v2.6.3 (Rambaut & Drummond, 2013s; Rambaut & Drummond, 2013b) from the posterior sampling of trees. All pipeline steps and phylogenomic analyses were conducted on the University of Florida HiPerGator high-performance computing cluster (http://www.hpc.ufl.edu/).