Differentiation patterns of emperor moths (Lepidoptera: Saturniidae: Saturniinae) of a continental island: divergent evolutionary history driven by Pleistocene glaciations

Material collection

Two sets of taxa for different analyses were assembled. First, to examine the monophyly of each Formosan genus and species belonging to the subfamily Saturniinae, one nuclear gene of 28S rDNA and two mitochondrial genes of COI and 16S rDNA were applied, and two members in Salassinae, one subfamily belonging to Saturniidae, were applied as outgroups. A total of 52 individuals of seven genera and 18 taxa from various parts of Taiwan, Vietnam and China were analyzed (Table S1). Second, the barcoding initiative provides a multitude of Southeast and East Asian saturniine COI sequences, COI sequences in the first dataset and those belonging to the Formosan saturniine genera from South and East Asia were downloaded from GenBank and the Barcode of Life Data System (BOLD) to elucidate the speciation timing and possible origin of Formosan saturniines. The specific names of the analyzed saturniine species follow the recent review in Kitching et al. (2018). Crucial information, collected in 519 sequences and the aligned sequences of 114 taxa is provided in Table S1 and Fig. S1. To prevent the misidentification of the downloaded sequences, information based on the monophyly and distribution area of each saturniine species was employed to accurately retrieve sequences from GenBank.

Voucher specimens used only in the present study were stored at −20 °C at the Department of Entomology, National Chung Hsing University. A document permitting the collection of saturniine insects was acquired from the Yushan, Shei-Pa, Taroko, and Kinmen national parks in Taiwan.

DNA extraction, amplification, and direct sequencing

One leg of each voucher specimen extracted for DNA was modified from Yeh et al. (2004) using the Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA). Crude DNA dissolved in 100 μl of ddH₂O was used as a template in subsequent polymerase chain reaction (PCR). Two mitochondrial genes, COI and 16S rDNA, and one nuclear gene, 28S rDNA, were amplified. K762UBombyx and Jerry D primers (Caterino & Sperling, 1999) (i.e., barcode primers) were used to amplify and sequence the COI fragment. The primer sets used for amplification and sequences of 16S rDNA and 28S rDNA were 16SR21-16S22 and 28Sa-28Sb, respectively (Yeh, Yang & Kang, 1997; Lin, Wang & Yeh, 2003). For COI, 35 cycles of amplification were performed in a final volume of 50 μl containing 10 mM Tris-Cl (pH 9.0), 50 mM KCl, 1.5 mM MgCl₂, 0.01% gelatin, 0.1% Triton-X100, 0.2 mM of each dNTP, 10 pmoles of each primer, 2 μl of the DNA template, and 2 units of SuperTaq polymerase (Protech Technology, Taipei, Taiwan). The PCR assays for 16S rDNA and 28S rDNA were performed in a volume of 25 μl mixture containing 1 μl of template DNA, 2.5 μl of 10× Taq buffer, 0.2 μl of 25 mM dNTP, 0.5 μl of Taq DNA polymerase, and 0.5 μl of primers. The PCR programming conditions were 95 °C for 2 min as the first denaturation followed by 35 cycles of denaturation at 95 °C for 1 min, annealing at 42–52 °C for 1 min, and extension at 72 °C for 50–60 s, with a final extension at 72 °C for 10 min. The amplified DNA fragment was excised and purified using the QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany) after resolution on agarose gel. The resultant DNA product was sequenced using a Taq dye terminator cycle sequencing kit (Applied Biosystems, Waltham, MA, USA) and ABI 377A sequencer. Sequences amplified in the present study have been deposited in GenBank under AB533547–AB533597 for the COI gene, AB841436–AB841480 for the 16S rDNA gene, and AB841488–AB841532 for the 28S rDNA gene (Table S1).

DNA analysis

Sequences of saturniines were compiled using the BioEdit program (Hall, 1999) and aligned with pertinent sequences retrieved from GenBank and the BOLD before being manually checked (Fig. S1). Pairwise distance was estimated using an uncorrected proportional divergence in MEGA7 software (Kumar, Stecher & Tamura, 2016).

Phylogenetic inferences of Formosan saturniines were made using the maximum likelihood (ML) approach and Bayesian inference (BI). The ML tree was executed using the online IQ-TREE 1.6.11 (http://iqtree.cibiv.univie.ac.at/) (Nguyen et al., 2015). The Edge-unlinked partition model was selected for the partition type because each amplicon has its own branch length. The best substitution model was selected using the “Auto” option and the FreeRate heterogeneity (+R). One hundred bootstrap replicates were used to test the statistical support of each node. MrBayes v.3.2 (Ronquist et al., 2012) was used to conduct BI, and the substitution model was estimated using jModelTest 0.1 (Posada, 2008). The evolutionary hypotheses for the best-fit substitution models of COI, 16S rDNA, and 28S rDNA were TIM1+I+G, TPM2uf+G, and TPM2+I+G, respectively, according to the Bayesian information criterion. The BIs of the three partitioned genes were analyzed with three heat chains and one cold chain, and Markov chain Monte Carlo searches were conducted for 10⁷ generations with sampling every 100 generations. The analysis was based on an average standard deviation of split frequencies lower than 0.01. The initial 25% of trees were discarded as burn-in, and the remaining trees were used to generate a consensus tree. Moreover, all the relevant Southeast and East Asian saturniine COI downloaded from GenBank were used to construct a ML tree using IQ-TREE 1.6.11 (http://iqtree.cibiv.univie.ac.at/) (Nguyen et al., 2015). The best substitution model was selected using the “Auto” option and the FreeRate heterogeneity (+R). One thousand Ultrafast bootstrap analyses were used to assess the statistical support of each node. Sibling species and those closely related to Formosan saturniine species were analyzed based on the reconstructed COI tree.

The haplotype network of each Formosan saturniine COI sequence was analyzed based on sibling or closely related saturniine species by using the TCS program with a 95% connection limitation (Clement, Posada & Crandall, 2000).

Calibration dating for the coalescent time of the separation between Formosan saturniines and their mainland or sibling saturniine lineages was estimated using a strict molecular clock in BEAST v.1.6.1 (Drummond & Rambaut, 2007). Different substitution rates for COI applied in insects were employed to estimate the relative divergent time since no fossil record has been recorded for Formosan saturniines and their sibling lineages (Sohn et al., 2012). Two substitution rates for COI applied in Lepidoptera were employed. A substitution rate of 1.15% per lineage in 1 million year for COI was used because it has been calibrated by the geological event and optimally applied for closely related taxa (Brower, 1994); While the substitution rate of 1.55% proposed by Nazari & Sperling (2007) was applied as well. Markov chain Monte Carlo sampling was conducted on 10⁸ generations with sampling every 10⁴ generations, depending on the effective sample size of the parameter estimated using Tracer software v.1.5 for each saturniine species (Rambaut & Drummond, 2009). The trees were recorded every 10,000 generations, and 10% of the recorded trees were burn-in.

Sequence composition of COI, 16S rDNA, and 28S rDNA in Saturniinae species

COI, 16S rDNA, and 28S rDNA were successfully amplified for 52 individuals of 18 taxa in 7 genera, with fragment sizes of 643, 486, and 654 bp, respectively. The downloaded sequences of COI and 16S rDNA from GenBank were 467 and 43, respectively. The average base compositions, including the downloaded sequences for G, A, T, and C, are 14.9%, 28.8%, 39.5%, and 16.8% in the COI gene; 13.9%, 39.8%, 38.7%, and 7.6% in 16S rDNA; and 30.6%, 25.8%, 20.8%, and 22.8% in 28S rDNA; respectively.

Sequence variations in Saturniinae species

The uncorrected nucleotide divergence and distribution frequencies in Saturniinae for each COI, 16S rDNA, and 28S rDNA in three categories were as follows: (1) 0–8%, 0–10%, and 0–1%, respectively, among individuals within species; (2) generally 6–11%, 0–11%, and 0–3%, respectively, among species of a given genus; and (3) generally 10–15%, 8–13%, and 2–3%, respectively, among genera (Fig. 1). Most saturniine species have highly similar COI and 16S rDNA sequences within species, with intraspecific divergence generally being <2%, whereas the interspecific and intergeneric genetic divergence was 6–12% (Figs. 1A and 1B). The divergence of the conserved 28S rDNA was generally low (<2%), and for many congeneric species, no genetic variations were detected (Fig. 1C).

Approximately 90% of uncorrected nucleotide divergences among COI within species were <2% (Fig. 1A); all divergences exceeding 2% originate from comparisons between allopatric populations or subspecies. For example, a considerable divergence of 4.8% was discovered in the population of Att. atlas in Taiwan from other allopatric populations of Att. atlas in Indonesia and Vietnam. By contrast, subtle or few genetic divergences were noted in allopatric species or subspecies (data provided below). The Taiwanese subspecies Saturnia japonica arisana (Shiraki) has an identical sequence to its allopatric nominate subspecies Saturnia japonica japonica (Moore) from China (data provided below).

Genetic divergences >0.5% among individuals have been observed in several Taiwanese lineages, such as Loepa mirandula Yen, Nässig, Naumann & Brechlin, Actias ninpoana ningtaiwana Brechlin, Att. atlas, Saturnia pyretorum Westwood, and Samia wangi Naumann & Peigler (data provided below).

Phylogenetic inferences of Saturniinae species

Phylogenetic inferences based on COI + 16S rDNA + 28S rDNA using BI and ML methods revealed that the sampled Saturniinae formed a monophyletic group comprising monophyletic genera and species (Fig. 2); this provides essential information to document the evolutionary differentiation events between sibling or closely related saturniines. Moreover, 519 COI sequences in 114 saturniine species related to Taiwanese Actias, Antheraea, Attacus, Loepa, Rhodinia, Samia, and Saturnia were downloaded from GenBank to analyze the differentiation history of each Taiwanese saturniine (Fig. S1, Table S1). The phylogenetic topology (Fig. S2) revealed that these Asian saturniine species are primarily monophyletic, although some of them exhibited paraphyletic deduced allopatric closely related species. For example, the endemic Loepa formosensis Mell and Rhodinia verecunda Inoue in Taiwan form terminal lineages within the widespread Southeast Asian Loepa diffunorientalis Brechlin and East Asian Rhodinia fugax Butler, respectively.

Phylogenetic relationships between insular Taiwanese saturniines and their Asian allies

On the basis of COI phylogenetic relationships (Fig. S2), the potential sister saturniine species, subspecies, and populations related to Taiwanese saturniines were used to estimate the differentiation and divergence events of each Taiwanese saturniine. However, no species closely related to Actias neidhoeferi Ong & Yu was discovered (Fig. S2).

Differentiation in the species category of Taiwanese saturniines

All six saturniine species endemic to Taiwan were found to be monophyetic. Among them, Loe. mirandula exhibited the greatest divergence from its sibling taxa and Rho. verecunda exhibited the least divergence. In Loe. mirandula, the network analysis revealed >30 substitution steps to its possibly sibling Loepa miranda Moore (Yen et al., 2000) and the closely related Loepa vanschaychi Brechlin, Loepa sinjaevi Brechlin, Loepa orientomiranda Brechlin & Kitching, Loepa obscuromarginata Naumann, Loepa roseomarginata Brechlin, and Loepa meyi Naumann (Fig. 3A). Three haplotypes with a maximum of five substitution steps were discovered in Loe. mirandula. For Loe. mirandula, the intraspecific sequence divergence was 0.5% on average, and the divergence from its allied Asian species was 5.3–6.5% (Fig. 3B). The differentiation event of Loe. mirandula from its sibling saturniine species was estimated to have occurred approximately 2.7 Mya (Fig. 3C).

For Taiwanese Rho. verecunda, two haplotypes with four substitution steps were found, differentiating Rho. fugax into two haplogroups (Fig. S3). Average sequence divergence within Rho. verecunda was 0.3%, and that from Rho. fugax and Rhodinia newara (Moore) was 0.7% and 3.7%, respectively (Fig. S3B). Network and phylogenetic analyses revealed that the Taiwanese endemic Rho. verecunda produces Rho. fugax paraphyly. The splitting events of Rho. verecunda and one of the Rho. fugax lineage occurred approximately 0.23 Mya (Fig. S3).

For Taiwanese Loe. formosensis, a network analysis revealed four haplotypes with a maximum of five substitution steps and >17 substitution steps to the closely related Loe. diffunorientalis distributed in Vietnam (Fig. S4A). On average, intraspecific sequence divergence was 0.4% within Loe. formosensis and 3.5–5.3% from its related saturniines Loe. diffunorientalis, Loepa diffunoccidentalis Brechlin, and Loepa nepalensis Brechlin (Fig. S4B). A speciation event of Loe. formosensis and its sibling Loe. diffunorientalis was estimated to have occurred approximately 1.5 Mya (Fig. S4C).

A network analysis revealed no variation in the Taiwanese species of Antheraea superba Inoue but substantial sequence divergence of 3.8% from its sister species Antheraea yamamai Guérin-Méneville from Japan and Korea and that from China (Fig. S5). The differentiation event of the two species can be dated back to approximately 1.3 Mya (Fig. S5C).

The network analysis also indicated that the Taiwanese Antheraea formosana Sonan groups with Indian Antheraea assamensis (Helfer) were related to Antheraea paniki Nässig & Treadaway from the Philippines with >16 substitution steps (Fig. S6A). Antheraea formosana was noted to exhibit an average sequence divergence of 1.1% and 3.4% from Ant. assamensis and Ant. paniki, respectively (Fig. S6B). The speciation events for these three saturniine species seems to have occurred <1 Mya, whereas those for Ant. formosana and Indian Ant. assamensis occurred <0.5 Mya (Fig. S6C).

Network analysis for Taiwanese endemic Saturnia fukudai Sonan revealed 23 substitution steps to its sibling lineage of Saturnia boisduvalii Eversmann from Korea (Fig. S7A), reflecting a sequence divergence of 4% between the two (Fig. S7B). Lineage calibration revealed that the estimated divergence was approximately 1.3 Mya between them (Fig. S7C).

Differentiation in the subspecies category of Taiwanese saturniines

Three of four Taiwanese endemic subspecific saturniines were monophyletic. Moreover, a genetic admixture caused Sat. jap. japonica to be paraphyletic. Actias ningpoana (C. & R. Felder), commonly known as the Chinese luna moth, distributes throughout East China and the islands of Taiwan and Hainan. The Taiwanese endemic subspecies Act. nin. ningtaiwana comprises two genetic lineages, one of which (lineage II) has a close affinity to individuals of the subspecies Actias ningpoana ningpoana (C. & R. Felder) from China (Fig. 4). The genetic admixture in the two subspecies led to high intrasubspecific sequence divergence within the subspecies Actias nin. ningtaiwana of up to 0.7% (Fig. 4B). The calibration time revealed that the divergence between Taiwanese lineages I and II occurred approximately 0.45 Mya (Fig. 4C). Taiwanese lineage I can be found only in Taiwan, whereas lineage II has an admixture relationship with the Chinese population, which split approximately 0.1 Mya.

Two Saturnia species from Taiwan were recognized as endemic subspecies. Network analysis for Taiwanese endemic subspecies of Saturnia thibeta okurai (Okano) revealed approximately 23 substitution steps to its neighboring subspecies of Saturnia thibeta pahangensis Brechlin and Saturnia thibeta extensa Butler from Malaysia and India, respectively (Fig. S8A). Sequence divergences between the Taiwanese and other subspecies were approximately 4% (Fig. S8B), and lineage calibration revealed a divergence of approximately 1.7 Mya between them (Fig. S8C).

For the Taiwanese endemic subspecies of Sat. jap. arisana, an identical haplotype has been found in several mainland Sat. jap. japonica individuals, indicating a recent migration event of this Taiwanese subspecies (Fig. S9A). Moreover, the lineage Sat. jap. arisana was found to have a paraphyletic relationship with Sat. jap. japonica (Fig. S9C). In addition, the complete mitochondrial genome of Saturnia boisduvalii (Eversmann) (Accession No. EF622227) might be a misidentification of Sat. jap. japonica (Fig. S9C).

In the Taiwanese endemic Act. sin. subaurea, three haplotypes with a maximum of three steps were discovered; however, only one haplotype existed in its sibling Actias sinensis sinensis (Walker) from China and Vietnam (Fig. S10A). The haplotype network revealed that both Act. sin. subaurea and Act. sin. sinensis have more than 28 substitution steps related to their potential sister species Actias parasinensis Brechlin from Bhutan (Fig. S10A). A distinct genetic divergence of 2.1% was noted between Act. sin. subaurea and Act. sin. sinensis (Fig. S10B), with a split event at approximately 0.9 Mya according to lineage calibration (Fig. S10C).

Differentiation in the Taiwanese allopatric population

Population differentiation analysis revealed that all Taiwanese saturniine populations belonged to monotypic species. A recent migration event for Taiwanese populations would thus be expected. Additionally, substantial divergence was found in the Taiwanese Att. atlas population. The attractive giant monotypic Att. atlas distributed widely throughout tropical regions, from India eastward to the Philippines and Taiwan. The network indicated that two divergent haplotypes exist in the Taiwanese Att. atlas, with sequence divergences of 2% and 4.8% from the China–Vietnam–India region and Indonesia, respectively (Fig. 5). The allopatric Att. atlas populations from the three regions formed distinct lineages (Fig. 5C). Lineage calibration estimated that the speciation event of Att. atlas occurred approximately 2.5 Mya. The Taiwanese lineage diverged from the China–Vietnam–India lineage approximately 1.2 Mya.

The monotypic Antheraea pernyi (Guérin-Méneville), which has seven haplotypes in Taiwan, China, and Korea, exhibits average sequence divergences of 1.7%, 2.4%, and 2.6% from the allied species Antheraea vietnamensis Brechlin & Paukstadt, Antheraea jawabaratensis Brechlin & Paukstadt, and Antheraea roylei Moore, respectively (Fig. S11). The differentiation event of the Taiwanese Ant. pernyi population could be dated back to approximately 0.4 Mya (Fig. S11C). Calibration dating indicated high affinity among Ant. pernyi, Ant. jawabaratensis, Ant. vietnamensis, and Ant. roylei and with differentiation <1.4 Mya. However, the molecular identification of Ant. pernyi and its allies, such as Ant. roylei and Ant. vietnamensis, have been complicated by the domesticated hybrids of these species that were used for silk production (Nagaraju & Jolly, 1986; Peigler & Naumann, 2003).

Saturnia pyretorum Westwood distributed in Taiwan and China and has two distinct lineages, each of which has divergent individuals (Fig. S12). Sequence divergence between the two lineages was 1.4–2.2%. Calibration dating revealed two split events between Taiwanese and mainland lineages, approximately 0.25 and 0.1 Mya, respectively, indicating multiple migrations of Sat. pyretorum to Taiwan.

Samia wangi Naumann & Peigler is a monotypic species distributed in southern China and Taiwan. The shared haplotype between the Taiwan and mainland China populations indicated that a dispersal event may have occurred recently (Fig. S13). Calibration dating also revealed a recent migration event. The sequence divergences of Sam. wangi to the two Samia species of Samia cynthia (Drury) and Samia ricini (Jones) were approximately 1.9% and 3.8%, respectively (Fig. S13).

Sequence divergence for insular saturniines to their allopatric populations and closely related alliances

The present molecular study indicates that the taxonomy of several saturniine species endemic to Taiwan should be reassessed. Species recognition would be indefinite if the analyzed species were under speciation processes. Moreover, recognizing one differentiation population as a subspecies within a species or as a distinct full species is also challenging. In the past two decades, DNA barcodes have provided much data, revealing that a 2% COI divergence can be usefully applied to identify species. In the current study, phylogenetic topology revealed that the analyzed saturniine species were primarily monophyletic, with some being paraphyletic (Fig. S2). Several Formosan saturniine lineages exhibited COI divergences >2% from their sibling lineages, whereas others had divergences <2%. To document the taxonomic debates pertaining to Taiwanese saturniines, we followed the recent taxonomic review of Saturniidae by Kitching et al. (2018). The frequencies of sequence divergences of each Taiwanese saturniine species, subspecies, and population from its closely related alliances are presented together in Fig. 6.

Sequence divergences within Formosan endemic species were <0.8% and >0.8% with their closely related alliances, except for Rho. verecunda (Fig. 6A). Anteraea superba has been considered a subspecies of Ant. yamamai, commonly known as the Japanese oak silkmoth, is a temperate saturniine distributed primarily throughout northern China, Korea, and Japan (Peigler & Wang, 1996; Zhu & Wang, 1996). These two species have highly similar morphological characteristics, although the wing color pattern varies from bright yellow to reddish brown among Ant. yamamai populations (Peigler & Wang, 1996). A substantial sequence divergence (>3.4%) from the sister species of Ant. yamamai indicates that Ant. superba in Taiwan is a separate evolutionary unit, either recognized in species or subspecies categories (Fig. 6A). Formosan Sat. fukudai has been considered an endemic subspecies of Sat. jonasii (Peigler & Wang, 1996; Zhu & Wang, 1996), and Kitching et al. (2018) treated it as a distinct species. Sequence divergences of approximately 4% in Sat. fukudai from its sibling species, Sat. jonasii and Sat. boisduvalii, suggest that Sat. fukudai in Taiwan is a separate evolutionary unit (Fig. 6A). Rhoinia verecunda, endemic to Taiwan, produced a paraphyletic Rho. fugax, involving several subspecies in China (Fig. S3). Intraspecific sequence divergences within Rho. verecunda were 0–0.8% and between Rho. verecunda and Rho. fugax were 0.2–1.4% (Fig. 6A); however, the color patterns and wing bandings between them could be distinguished, thereby indicating inconsistencies between morphological characteristics and molecular divergence.

Sequence divergences were <2% within Formosan subspecies, >2.2% from their closely related species, and 0–4% from their allied subspecies (Fig. 6B). Actias nin. ningtaiwana, which is endemic to Taiwan, is known subspecies of Actias selene ningpoana C. & R. Felder of Actias selene (Hübner), more commonly known as the Indian luna moth, which is distributed throughout South Asia (Peigler & Wang, 1996). Brechlin (2012) highlighted that several Act. ningpoana subspecies were distinct from Act. selene. The present study identified sequence divergences larger than 2.2% (Fig. 6B) and a splitting time between Act. ningpoana and Act. laovieta Naumann, Nässig & Löffler of approximately 1.0 Mya (Fig. 4). Genetic variations with an overlapping of sequence divergences between intersubspecies and intrasubspecies were observed: 0–2% within Act. nin. ningtaiwana, 0–2% between Act. nin. ningtaiwana and Act. nin. ningpoana, and 0.7–2.1% between Act. nin. ningpoana and Act. nin. ningtaibaishana; therefore, their subspecific nomenclature was open to debate (Fig. 6B). For Act. sin. subaurea, Zhu & Wang (1996) suggested that Act. sin. subaurea should be a species distinct from Act. sinensis and be called Actias heterogyna (Mell). The present study revealed that a distinct genetic divergence of 2–2.4% existed between Act. sin. subaurea and Act. sin. sinensis, with a split event approximately 0.9 Mya (Fig. S10). Therefore, Act. sin. subaurea in Taiwan is a separate evolutionary unit that is either recognizable as a species or subspecies category (Fig. S10). Wing color and pattern are similar among Sat. thibeta subspecies, although the wing bandings are more pointed in Formosan Sat. thi. okurai (Roepke, 1953; Pinratana & Lampe, 1990; Peigler & Wang, 1996; Zhu & Wang, 1996). Sequence divergences of approximately 4% in Sat. thi. okurai to its related subspecies of Sat. thi. pathangensis may raise the issue of species splitting (Fig. 6B). A haplotype of a Formosan subspecies of Sat. jap. arisana identical to some mainland Sat. jap. japonica individuals shows no variation (Fig. 6B), although wing banding is clear and thicker in Sat. jap. arisana (Roepke, 1953; Pinratana & Lampe, 1990; Peigler & Wang, 1996; Zhu & Wang, 1996).

At the population level, overlapping sequence divergences were generally found within and between populations, which might be caused by recent dispersal or multiple migration events such as those for Sat. pyretorum and Sam. wangi (Fig. 6C); no taxonomic debate has been proposed except for one concerning in Att. atlas. The attractive giant Att. atlas is distributed widely throughout tropical regions, from India eastward to the Philippines and Taiwan, and several subspecies of varying size, coloration, and wing pattern have been proposed (Roepke, 1953; Pinratana & Lampe, 1990; Peigler & Wang, 1996; Zhu & Wang, 1996). However, Kitching et al. (2018) suggested that the proposed Att. atlas subspecies should be considered synonymous with Att. atlas. Sequence divergences with ranges of 1.2–5.2% among several allopatric Att. atlas populations suggest distinct evolutionary units among them (Fig. 6C). For the Formosan Sat. pyretorum, two distinct genetic lineages were likely to have multiple dispersal events that produces the genetic admixture between insular and mainland populations (Fig. 6C, Fig. S12); however, Sat. pyretorum was artificially introduced to Taiwan in the first half of the 20th century (Chu, 2005). Accordingly, the divergence between the natural Taiwanese population and the mainland population should be 0.8–1.4%. No sequence variation was found in some individuals of Sam. wangi from Taiwan and the mainland, implying a recent migration event (Fig. 6C).

Calibration dating and origin of island saturniines

No fossil records have been reported for Taiwanese saturniines and their close affinities to infer their absolute divergent time. Molecular calibration based on COI sequences for East and Southeast Asian saturniines can help in determine the spatiotemporal patterns of insular saturniines. Speculation suggests that organisms in Taiwan primarily originated from neighboring mainland China, Southeast Asia, and Northeast Asia, either through currents or the land bridges formed during glaciation (Voris, 2000; Ray & Adams, 2001; Chang & Chen, 2005; Jang-Liaw, Lee & Chou, 2008; Huang & Lin, 2010; Tsai, Wan & Yeh, 2014; Yu, Lin & Weng, 2014; Tsai & Yeh, 2016; Weng, Yang & Yeh, 2016; Kurita et al., 2017). However, phylogenetic analyses indicated that most Formosan saturniines were distinct and monophyletic. Therefore, migration through currents across the Taiwan Strait might not be the most common route of saturniines’ arrival to Taiwan.

Figure 7 present the split times accompanied by the geographical distribution of each Formosan saturniine species, subspecies, and population from its neighboring closely related alliances with Pleistocene epochs. In six Taiwan endemic saturniine species, calibration dating yielded sporadic splitting time within the range of 0.2–2.7 Mya (Fig. 7A). Loepa mirandula, Loe. formosensis, Ant. superba, and Sat. fukudai exhibited split events of >1.3 Mya (Early Pleistocene) and <0.5 Mya (Middle Pleistocene) in Ant. formosana and Rho. verecunda. Furthermoer, two subspecific endemic saturniines in Taiwan had splitting events 0.9–1.7 Mya (Early Pleistocene) and another two, Sat. ningpoana ningtaiwana and Sat. japonica arisana, might have had splitting events in the Middle and Late Pleistocene (Fig. 7B). However, the estimated divergence time of the subspecific Formosan Sat. thi. okurai might be <1.71 Mya when including mainland populations (Fig. S8); this is based on the fact that Sat. thibeta also exists in northern Vietnam and northern China (Peigler & Naumann, 2003). At the population level, the splitting time of Formosan populations from their mainland populations was estimated to be <0.4 Mya in the Middle and Late Pleistocene, except for Att. atlas, which had a splitting time of approximately 1.2 Mya (Early Pleistocene) (Fig. 7C). Kawahara et al. (2019) used phylogenomics to determine the evolutionary timing and pattern of butterflies and moths including Actias, Saturnia, and Attacus genera of Saturniinae. Phylogenomic analysis also revealed the sister relationships between Actias and Saturnia (Fig. 2). Figure 7A also indicates that the splitting time, 0.2–2.7 Mya, in species differentiation events was far less than the generic divergent time, 22–30 Mya, as proposed by Kawahara et al. (2019).

Merckx et al. (2015) indicated that the montane adaptive biota in Mount Kinabalu, at an elevation of 4,095 m, on the tropical island Sabah mostly originated from long-range dispersal from cool localities elsewhere. The present study indicates that some insular saturniines might have originated from temperate regions in Northeast Asia, whereas others may have originated from southern China and Southeast Asia (Fig. 7). Taiwanese saturniines, including Ant. superba, Sat. japonica arisana, Sat. fukudai, and Ant. pernyi, predominantly inhabited montane or subalpine regions in Taiwan at altitudes >1,500 m and have related lineages or sibling species distributed mostly in temperate areas such as northern China, Japan, and Korea (d’Abrera, 1998). Taiwanese saturniines distributed throughout plains and hills might not be related to temperate lineages. All species, subspecies, and populations of the Taiwanese Act. sinensis, Ant. formosana, Att atlas, Loe. formosensis, Loe. mirandula, and Sat. thibeta have closely related lineages distributed throughout India, Thailand, Vietnam, Malaysia, Indonesia, and southern China (Peigler & Wang, 1996). Some Taiwanese saturniines exhibited strong genetic relationships with saturniine lineages in southern China, including Rho. verecunda, Act. nin. ningtaiwana, Sat. pyretorum, and Sam. wangi (Fig. 7). Moreover, periodical Pleistocene glaciations may have led to the repeated colonization of Act. nin. ningtaiwana in Taiwan, as indicated by the coexistence of two divergent lineages in Taiwan (Figs. 4 and 7). However, artificial genetic admixtures occurring in Sat. pyretorum may explain the phenomenon of mixing genetic lineages. Because of a history of domestication for silk production, a genetic admixture between Sam. ricini and Sam. canningi (Hutton) as well as between Ant. pernyi and Ant. roylei (Peigler & Naumann, 2003) (Figs. S11 and S13) was evident; additionally, the domesticated hybrid Ant. proylei was reported to likely be derived from Ant. roylei and Ant. pernyi (Nagaraju & Jolly, 1986).

To address the speciation events of Lepidopteran Heliconius erato in the Amazonian region, Brower (1994) applied COI divergences and estimated the separation time to be 1.5–2 Mya. In the model insect of Drosophila, Coyne & Orr (1997) indicated that the time required for full reproductive isolation was generally 1.5–3.5 Mya. The application of COI sequences in molecular calibration has demonstrated that average speciation times were 1–2 Mya for the temperate biota of Europe and North America and approximately 2.6 Mya for the tropical Amazonian biota (Avise, 2000; Hewitt, 2000; Knowles, 2000; Rull, 2008; Hoorn et al., 2010; Rull, 2011). Studies including COI sequences have reported that the split time for most insects in Taiwanese speciation events was 1.0–2.8 Mya in the Early Pleistocene (Yeh et al., 2004; Huang & Lin, 2010; Tsai, Wan & Yeh, 2014; Tsai & Yeh, 2016; Weng, Yang & Yeh, 2016), although a few cases can be dated back to approximately 4 Mya (Weng, Yang & Yeh, 2016) or <0.3 Mya (Weng, Yeh & Yang, 2016) in the Middle Pleistocene.

The present study revealed that the split time for most endemic Taiwanese saturniine species and subspecies was 0.5–1.7 Mya in the Pleistocene. The inconsistency between the morphological characteristics and molecular divergence of the insular Sat. jap. arisana and Rho. verecunda may be attributable to genetic drift or local adaption when geographic discontinuity occurred in the absence of a land bridge (Yeh et al., 2004; Tsao & Yeh, 2008; Wu et al., 2010; Tsai, Wan & Yeh, 2014; Tsai & Yeh, 2016; Weng, Yang & Yeh, 2016). Moreover, this study was based on small datasets, which might possibly present only a partial differentiation history of Taiwanese saturniines. The differentiation history of each Taiwanese satruniid will be clearer, if genomic information is included in future analyses. However, studies over the past three decades have included small datasets and have contributed critical findings in systematics.