Morphological identifications in this study were made based on illustrations of type material in Moore (1857), Melvill (1893), Leech (1892), de Niceville (1897), Seitz (1909, 1927), Tytler (1915), Oberthür (1920), Lang (2012) and Sondhi et al. (2016). Male genitalia were dissected, prepared and photographed using traditional techniques (Supplementary Material S1: Sbordoni V. personal communication).

DNA was extracted from one leg of each individual on Biomek FX liquid handling robot using a semi-automated DNA extraction protocol (Ivanova et al. 2006) on glass fiber plates (PALL Acroprep 96 with 3µm GF membrane over 0.2µm Bioinert membrane). DNA was eluted in 35-40µl of ddH₂0 pre-warmed to 56°C and stored at -80°C. Polymerase chain reactions (PCR) and DNA sequencing were carried out following standard procedures for Lepidoptera (Hajibabaei et al. 2005, deWaard et al. 2008). A 658-bp fragment of cytochrome c oxidase subunit I (COI) was targeted for amplification using the primers LepF (5’ ATTCAACCAATCATAAAGATATTGG3’) and LepR (5’TAAACTTCTGGATGTCCAAAAAATCA3’).

The entire nuclear RpS5 was initially amplified using the primers RpS5f/r (Wahlberg and Wheat 2008) and afterwards a fragment covering 483 bp was obtained by nested PCR using RpS5f and a novel specific primer (5’ACMGTDCGYCGTCARGCTGTRGAYGTBTCWCC3’), developed from conserved regions. The PCR conditions were carried out as in Todisco et al. (2010).

Sequences were obtained by using an ABI 3730xl sequencer following the manufacturer’s recommendations and they were edited and assembled using CodonCode Aligner 6.0.2. All sequences of Calinaga were submitted to GenBank (Accession Numbers in Table 1).

Two sequences from C. buddha (COI: EF683684, AY090208; RpS5: EU141406) and one from C. davidis (COI: HQ658143), as well as COI and RpS5 sequences for eight outgroups for species belonging to closely related lineages of Nymphalidae (Charaxinae and Satyrinae) from Wahlberg and Wheat (2008) were retrieved from GenBank and added to our dataset (see Table 1).

We used NETWORK 5 (Bandelt et al. 1999) to calculate a Median Joining (MJ) network representing the genealogical relationships among mtDNA haplotypes (Fig. 1B). Genetic distances were calculated in MEGA 6 under the Kimura 2-parameter model of base substitution (Tamura et al. 2013). The program JModelTest 2.1.7 (Guindon and Gascuel 2003, Darriba et al. 2012) was employed to evaluate models of nucleotide substitutions. Only the best-fit models were subsequently used for estimations of tree topology.

We used the BEAST package 1.8.3 (Drummond et al. 2012) for Bayesian phylogeny reconstruction and divergence time estimation for concatenated genes. The analysis was allowed to run for 50 million generations and the first 5000 trees were discarded as burnin. A consensus tree was calculated with posterior probability limit of 0.5, maximum clade credibility tree and median heights. Monophyly was enforced on a) Satyrinae, b) Charaxinae, and c) Satyrinae + Charaxinae (Fig. 2). The tree prior was set to Yule Process Speciation model with a random starting tree. The analysis was performed using an uncorrelated relaxed clock (lognormal relaxed distribution not selected) for both genes. Trees were visualized using FIGTREE 1.4 (Rambaut 2009).

Recent phylogenetic studies calibrated with fossil data (Peña and Wahlberg 2008; Wahlberg et al. 2009) suggest that Nymphalidae began to diversify in the late Cretaceous around 90 Mya, satyrines (Calinaginae + Satyrinae + Charaxinae) around 75 Mya and Satyrinae + Charaxinae about 70 Mya. Thus, in our BEAST analyses the tree root was calibrated with normal distribution at 75 Mya (stdev 3), Satyrinae + Charaxinae at 70 Mya (stdev 3.5) and Charaxes + Euxanthe at 22 Mya (stdev 1) (see Fig. 2).

Figure 2. 

Bayesian phylogeny of Calinaga estimated in BEAST using concatenated data. Purple squares are calibration points (root: 75 ± 3; Satyrinae + Charaxinae 70 ± 3.5, Charaxes + Euxanthe 22 ± 1). Monophyly was enforced on nodes marked with orange squares. The inset map shows the biogeographic regions used in DIVA analysis: A) Southwestern China ecozone, B) Himalaya-Tibetan plateau region, C) Northern Sino-Himalaya, D) Southern Sino-Himalaya, E) Indochina. Colored dots correspond to haplogroups on the tree.

Haplotype and nucleotide diversity were calculated using DnaSP 5.0 (Librado and Rozas 2009). Demographic equilibrium was tested in each mtDNA haplogroup by calculating F_S (Fu 1997) and R₂ (Ramos-Onsins and Rozas 2002) statistics, which have been shown to provide sensitive signals of population expansion (Ramos-Onsins and Rozas 2002). Arlequin 3.0 (Excoffier et al. 2005) and DnaSP 5.0 were used to compute F_S and R₂, respectively, and to test their statistical significance by simulating random samples (10,000 replicates) under the null hypothesis of selective neutrality and constant population size using coalescent algorithms (both modified from Hudson 1990). P-values for the two statistics were obtained as the proportion of simulated values smaller than or equal to the observed values (α= 0.05). Expected mismatch distributions and parameters of sudden expansion τ= 2 μ T were calculated using Arlequin 3.0 by a generalized least-squares approach (Schneider and Excoffier 1999), under models of pure demographic expansion and spatial expansion (Ray et al. 2003, Excoffier 2004). The probability of the data under the given model was assessed by the goodness-of-fit test implemented in Arlequin 3.0. Parameter confidence limits were calculated in Arlequin 3.0 through a parametric bootstrap (1000 simulated random samples).

Dispersal-vicariance optimization of ancestral areas implemented in DIVA 1.1 (Ronquist 1996) was used to infer the biogeographic history of the group. The distributions of populations were divided into five biogeographic units designated following Che et al. (2010) and Procheş and Ramdhani (2012), as follows: A) Southwestern China ecozone, B) Himalaya-Tibetan plateau region, C) Northern Sino-Himalaya, D) Southern Sino-Himalaya, E) Indochina (Fig. 2). The least number of dispersal events (12) was achieved when maxareas was set to 3.

The morphological study on the wing pattern of our specimens, based on illustrations in the most relevant literature, suggested the presence of four species: C. aborica, C. davidis, C. lhatso and C. sudassana (Table 1). Many of the samples procured as C. buddha proved to be misidentifications of either C. davidis or C. sudassana. In addition, the samples of C. buddha (COI: EF683684, AY090208; RpS5: EU141406) and C. davidis (COI: HQ658143) retrieved from GenBank were re-identified respectively as C. davidis and (possibly) C. lactoris (Table 1).

Male genitalia in Calinaga were found to be very uniform and thus genitalic characters were of little use in discrimination of species. Although Lang (2012) emphasized the shape and ‘pointiness’ of the dorsal apex of the valve as an important character, we did not find this trait useful or consistent within species. All the adult and genitalia images are available in Supplementary Material S1.

For COI (658 bp), 54 sequences (51 sequences analyzed in this study plus 3 sequences retrieved from GenBank) showed 23 different haplotypes (Fig. 1B; Table 1) with 71 variable sites. The 24 RpS5 sequences (483 bp) displayed 4 haplotypes and 6 variable sites.

Bayesian phylogeny analysis was used to reconstruct phylogenetic relationships (Fig. 2) for the concatenated genes. The resulting tree was based on an alignment of 1310 bp and it was rooted using outgroup sequences of members of two subfamilies: Charaxinae (Charaxes castor, Anaea troglodyte, Archaeoprepona demophon, Euxanthe wakefieldi) and Satyrinae (Caligo telamonius, Morpho peleides, Bicyclus anynana, Melanitis leda) (Table 1). The selected model of evolution for RpS5 nuclear gene was GTR+G+I (Rodriguez et al. 1990) with estimated base frequencies, 4 gamma categories, partitioned into positions (1+2),3 and for COI was HKY (Hasegawa et al. 1985) with estimated base frequencies, no heterogeneity model and not partitioned.

The Bayesian tree for concatenated genes resulted in a well-supported Calinaga clade (Fig. 2). Beside distinct lineages for C. lhatso (N= 1) and C. aborica (N= 2), two major clades were observed. The first clade, consisting of specimens identified as C. buddha (N= 20), C. davidis (N= 6) and C. buphonas (N= 2), was called Haplogroup A “davidis” because all samples were re-identified as C. davidis (Table 1). The latter also included two GenBank sequences for C. buddha (COI: EF683684, AY090208; RpS5: EU141406). Examination of the voucher images for these sequences showed that these also were indeed C. davidis (Min et al. 2008; Peña and Malm 2012). The second clade, consisted of specimens initially identified as C. sudassana (N= 14), C. buddha (N= 5) and C. lhatso (N= 1) and it was called Haplogroup B “sudassana” because all samples were subsequently re-identified as C. sudassana (Table 1).

The MJ network analysis (Fig. 1.B) of the mtDNA haplotypes was fully consistent with our Bayesian phylogenetic analysis. Forty-five mutational steps joined Haplogroup A “davidis” (including all C. davidis sequences; Fig. 1B) and Haplogroup B “sudassana” (including all C. sudassana sequences with C. lhatso sequences CLHA-CHSFT1). The genetic differentiation between the two haplogroups, estimated as the mean distance (Kimura 2 parameter uncorrected p-distance) between groups, was 10.3%. The analysis also highlighted a star-like configuration of Haplogroup A “davidis” with the H₁₅ haplotype shared among Yunnan populations localized in the headwater area for the three major rivers flowing from Tibet (Salween, Mekong and Yangtse; Fig. 1A, Fig. 2). In this geographical area Callinga displays its highest genetic variability, and Haplogroup A “davidis” and B “sudassana” coexist (Fig. 2). Thirty-four and thirty-three mutational steps joined C. aborica with Haplogroup A and B respectively. In addition, in Haplogroup A “davidis”, a GenBank sequence for C. davidis (HQ658143) from Tianma National Nature Reserve in Jinzhai County, Anhui, China (Xia and Hao 2011), possibly representing the taxon C. lactoris Fruhstorfer 1908, was divergent (Fig. 1B). Finally, in the Haplogroup B “sudassana”, the single specimen CLHA-CHSFT1 from Tsinling Mts, Foping nature reserve in Shaanxi, initially identified as C. lhatso based on morphological similarity, showed identical COI and RpS5 sequences to C. sudassana. Extraction and amplification of this sample was conducted twice in two different laboratories and generated the same results. Eleven mutational steps (Fig. 1B) separated the DNA barcode sequence of this individual from the other C. lhatso from Yunnan (CLHA-CHYZH1) despite their similar wing patterns.

According to Fs and R₂ statistics, calculated for Haplogroup A “davidis” (N= 20) and Haplogroup B “sudassana” (N= 30) mtDNA sequence sets (Table 2), the null hypothesis of constant population size could be rejected for these two phylogeographic units (bold in Table 2). Mismatch distribution of these groups was examined according to sudden expansion model (Table 2), and goodness of fit tests did not show significant deviations from expected distributions, so that parameter τ= 2μt could be used to estimate the time (t) elapsed from population expansion: estimated values of τ and their 5% and 95% confidence limits are shown in Table 2. According to Brower (1994) proposed mutation rates in Lepidoptera (µ= 0.01 substitutions/site/lineage/Ma), demographic expansion of the Haplogroup B “sudassana” could be traced back to t= 260 kya (126-403 kya). A similar estimate was obtained for Haplogroup A “davidis” t= 168 kya (70-362 kya).

The most recent common ancestor of Calinaga has been estimated to have split from Satyrinae + Charaxinae sometime between 67-84 Mya (Wahlberg et al. 2009). Results of our BEAST analysis reveal that the first major split within Calinaga occurred in the Eocene, around 43 Mya (63-25 Mya) (Fig. 2). The time of the most recent ancestor (tMRCA) of C. aborica + “davidis group” was estimated to have occurred in the Oligocene ~30 Mya (47-17 Ma), whereas it was ~21 Mya (38-8 Ma) for C. lhatso + “sudassana group”. Finally, the tMRCA for both davidis and sudassana groups was estimated at ~12 Mya (24-4 Ma).

Results of the ancestral area reconstructions showed that the common ancestor of Calinaga probably occurred before ~43 Mya in area ADE (Fig. 2), corresponding to Southwestern China ecozone, Southern Sino-Himalaya and Indochina. Before ~30 Mya this also represented the ancestral area for clade C. aborica + Haplogroup A “davidis”. Subsequently in the Miocene before ~15 Mya, in a more restricted area ranging between the Southwestern China ecozone and Southern Sino-Himalaya (AD), the common ancestor of Haplogroup A “davidis” was present with subsequent dispersal (in the late Miocene and Pleistocene) into the Himalaya-Tibetan Plateau region. The common ancestor of clade C. lhatso + Haplogroup B “sudassana” colonized the Southern Sino-Himalaya and Indochina (AE) before ~20 Mya. Finally, before ~10 Mya Indochina (E) represented the ancestral area of the Haplogroup B “sudassana”.

Despite six species initially identified (Fig. 1, Table 1), our molecular results identified only four main groups (C. aborica, C. sudassana, C. davidis and C. lhatso). These were also consistent with subsequent careful morphological identifications (Fig. 2, Table 1). This is quite reflective of the persistent confusion about the taxonomic status and number of species within Calinaga. Historically, various authors have recognized different numbers and sets of valid Calinaga species, from one (Ehrlich 1958b) or three (Fruhstorfer 1914) to seven (D’Abrera 1985, 1993, Beccaloni et al. 2003), ten (Savela 2015) or eleven species (Sondhi et al. 2016). Many of the samples we procured as C. buddha proved to be misidentifications of either C. davidis or C. sudassana. A thorough examination of the type material of all the names under Calinaga is required to clarify whether some of these names should be synonymized. Within our examined specimens we found that C. lhatso and C. aborica were morphologically and genetically distinct entities. However, the identity of the specimen identified as C. lhatso (CLHA-CHSFT1), sharing its haplotype with C. sudassana, remains unclear. Many of the specimens initially identified to sub-species under C. davidis showed little or no genetic differentiation (e.g. C. davidis nubilosa, C. davidis genestieri, C. davidis cercyon, C. davidis buphonas etc.).

The question of where and when Calinaga originated and how its species are related have not been previously investigated. Although no fossils are known for the genus, previous phylogenetic studies on Nymphalidae that have used fossils for molecular calibration have dated the last common ancestor of Calinaga + (Satyrinae + Charaxinae) at 67-84 Mya (Peña and Wahlberg 2008; Wahlberg et al. 2009). Since these studies used only a single Calinaga species, the age of the crown group of Calinaga and divergences within the genus have remained unknown. Our results show that the first split in Calinaga occurred in Eocene at around 43 Mya (Fig. 2). The large gap between this date and the age of the most recent common ancestor of Calinaga and Satyrinae + Charaxinae (67-84 Mya) suggests diversification after the Cretaceous/Tertiary boundary with subsequent extinction of sister lineages (Ree and Smith 2008). The initial collision between India and Asia in the early Cenozoic (50-55 Ma) set in motion major orogenic associated environmental changes that continued through the Oligocene, Miocene and into modern times which were likely a major driving force in evolution of Calinaga, and numerous other plant (e.g. Liu et al. 2002, Mao et al. 2010) and animal lineages (e.g. Rüber et al. 2004, Zhang et al. 2006, Nazari et al. 2007, Che et al. 2010, Zhao et al. 2015) found in the Himalayan-Tibetan Plateau. In fact, the evolutionary history of organisms in this region helps elucidate parallel geological and climatic developments, and often aids the recognition of unrecorded geological events.

Our results suggest that the common ancestor of Calinaga probably occurred before ~43 Mya (63-25 Mya) in the region now corresponding to the Southwestern China ecozone, Southern Sino-Himalaya and Indochina (ADE: Fig. 2). Since the Indochinese peninsula did not exist before ~22 Mya (Leloup et al. 2001), the ancestral distribution of Calinaga was likely confined to the Southwestern China ecozone (D) and Southern Sino-Himalaya (A). Since its center of genetic diversity seems to lie on the south-eastern border of the Tibetan Plateau, in particular where the three major rivers sourced in Tibet (Salween, Mekong and Yangtse rivers; Fig. 1A, Fig. 2) are in proximity, Calinaga may well have originated in the Tibetan region. Geological records indicate that around the Oligocene/Miocene boundary (~23 Mya), the rapid extrusion of Indochina and the cumulative Cenozoic deformation of South-East Asia caused by underthrusting of the Indian Plate beneath the Eurasian Plate, resulted in intense orogeny in South-East Tibet (Zhang et al. 2013). In particular in the Early Miocene (~20 Mya), driven by rapid deformation of the eastern Himalayan syntaxis (Robinson et al. 2014), large modifications were induced in the river systems of this region (e.g. Salween, Mekong, Yangtze and Yarlung Tsangpo-Irrawaddy river: Zhang and Hu 2012). These geological upheavals have likely provided ample opportunities for allopatric isolation and subsequent speciation in Calinaga.

Our molecular dating suggests that C. lhatso and Haplogroup B“sudassana” split into distinct lineages in Early Miocene ~21 Mya (38-8 Mya), possibly isolating the last common ancestor of C. lhatso in the Southern Sino-Himalaya (A). Southern Sino-Himalaya and Indochina (AE) possibly represented the ancestral area for C. lhatso + Haplogroup B“sudassana”. Considering that before ~22 Mya the Indochinese peninsula did not exist, it is more likely that the ancestral area of this clade was mainly in Southern Sino-Himalaya (A). The common ancestor of Haplogroup B“sudassana” began diversifying around ~12 Mya (24-4 Mya) and it probably expanded its range from Southern Sino-Himalayan region (A) to Indochina peninsula (E). Conversely, the origin of C. aborica + Haplogroup A “davidis” must be traced back to Oligocene (~30 Mya), when the first uplift of the Himalayan–Tibetan Plateau occurred (29 – 65 Mya: Zhao et al., 2015). Afterwards, in Miocene ~12 Mya (21-4 Mya), Haplogroup A “davidis” diversified and expanded its range into Central-East China.

In addition, our analyses show that Haplogroup A “davidis” and Haplogroup B“sudassana” experienced demographic expansion in the Pleistocene, ~168 kya (70-363 kya) and ~260 ka (16-403 kya) respectively. Therefore, the genetic diversity and the present-day distributions of Calinaga should be attributed not only to the dramatic geological events in Oligocene/Miocene, but also to the impacts of the cycles of glacial advance/retreat that occurred in this region during the Pleistocene (Quan et al. 2014).

Finally, the taxon C. formosana, for which we examined several specimens but could not obtain any sequences (Supplementary Material S1), is morphologically allied to the Haplogroup A “davidis”. Since Taiwan formed ~4-5 Mya at a complex convergent boundary between the Philippine Sea Plate and the Eurasian Plate (“The Penglai Orogeny”: Teng 1990, Liu et al. 2000), the most plausible hypothesis is that C. formosana must have dispersed from mainland China to Taiwan recently in Pleistocene.