Systematics and evolution of the African butterfly genus Mylothris (Lepidoptera, Pieridae)

Abstract. We study the systematics and evolutionary history of the Afrotropical butterfly genus Mylothris (Lepidoptera: Pieridae) based on six gene regions (COI, EF1a, GAPDH, MDH, RpS5 and wingless). We find that the genus can be placed into five species groups, termed the jacksoni, elodina, rhodope, agathina and hilara groups. Within these species groups, we find that many species show very little genetic differentiation based on the markers we sequenced, suggesting they have undergone rapid and recent speciation. Based on secondary calibrations, we estimate the age of the crown group of Mylothris to be about 16 million years old, but that many of the species level divergences have happened in the Pleistocene. We infer that the clade has its origin in the forests of the Eastern part of Central Africa, and has spread out from there to other regions of Africa.


Introduction
Studies on the systematics of Afrotropical butterflies using molecular data have been increasing in number in the past decade (Kodandaramaiah and Wahlberg 2007;Aduse-Poku et al. 2009;Nazari et al. 2011;Price et al. 2011;van Velzen et al. 2013;Aduse-Poku et al. 2017), giving us a better understanding of their evolutionary histories in one of the diversity hotspots of the planet. Many genera are endemic to the continent, although there are clear connections to the Oriental region (Kodandaramaiah and Wahlberg 2007;Aduse-Poku et al. 2009;Aduse-Poku et al. 2015;Sahoo et al. 2018;Toussaint et al. 2019) and even to the Neotropical region (Silva Brandão et al. 2008). Genera in the family Pieridae have not been the focus of many studies to date. So far only two pierid genera have been studied using molecular approaches in any detail, Colotis Hübner, 1819 (Nazari et al. 2011) and Pseudopontia Plötz, 1870 (Mitter et al. 2011). Here we investigate the systematics and evolutionary history of the Afrotropical pierid genus Mylothris Hübner, 1819.
The range of Mylothris includes most of the African continent south of the Sahara, extending to nearby islands such as Madagascar and the Comoros in the east, and to Bioko and São Tomé e Príncipe in the west. A single species flies alongside other pierid species belonging to the Afrotropical fauna in south west Arabia. It is speciose: a revision of the genus by one of the authors (HWG) is in progress and will list 105 species, a number of them new, many of them in turn divided into several subspecies (Warren-Gash in press).
A distinguishing feature of the genus, as with the related Catasticta Butler, 1870 in the Neotropics (Braby and Nishida 2010) and Delias Hübner, 1819 in the Australasian region (Braby 2006), is that the hostplants, with one known exception, are members of parasitic mistletoes belonging to the family Loranthaceae. This in turn determines where and how plentifully Mylothris species are found. Their foodplants are arboreal epiphytes and so they are insects of the canopies -unlike most Afrotropical pierids, which breed on low-growing shrubs such as Capparis Linnaeus, 1758, and many of which occur in more open, drier biotopes. Thus Mylothris species diversity is richest in the main African forest belt, with some species spreading into the forest-savannah mosaic and anthropogenic biomes including suburban gardens.
Two other consequences follow. The first is that, since Loranthaceae occur most frequently near the crown of the tree, many of the tropical forest species fly high, and some of them are seldom seen, accentuating their rarity in institutional collections. Secondly, they are particularly vulnerable to primary forest degradation.
Typically, Mylothris eggs are laid in a group on the underside of a leaf of the hostplant, and the larvae when they hatch are gregarious. The length of the lifecycle depends on the species and the biome or region. In southern Africa they follow the seasons in line with their foodplant, as they do in colder (usually montane) climates further north. In less seasonal biotopes, adults can be on the wing at any time of year (Warren-Gash in press).
In appearance Mylothris are a remarkably homogeneous group. In adults of the M. jacksoni species group, the hindwing is yellow and the forewing usually mostly white. The other clades typically show some orange (sometimes yellow) at the base of the forewing on one or both surfaces becoming white distally, with a greater or lesser black distal margin. There are exceptions, but they are similar enough to have baffled taxonomists over the years, leading to groups of species being lumped together under a single name in some cases and, no less frequently, names erected to describe what are no more than individual variations. If that were not enough, cryptic rings of different species that in adult morphology are externally almost indistinguishable, can be found flying together, most frequently in the two main areas of diversity: the Albertine Rift and central Cameroon (Warren-Gash in press).
The forthcoming revision (Warren-Gash in press) will look at all these issues in more detail based on morphological characters. Here we aim to investigate whether molecular data can shed light on species delimitations and their relationships. As yet undescribed species are here named nsp1 to nsp4. They will be formally described in Warren-Gash (in press).

Material and methods
We attempted to sample as many species as possible, and for selected species, as many populations as possible. We included 52 out of 105 species and a total of 235 individuals in our analyses (see Suppl. material 1: Table S1). Most specimens come from the private collection of HWG, with a few specimens taken from the collections of the African Butterfly Research Institute in Nairobi, Kenya. The specimens are part of a large taxonomic revision of Mylothris (Warren-Gash in press). In addition we included 23 outgroups from the subtribe Aporiina, taken from a previous publication (Wahlberg et al. 2014).
DNA was extracted from 2-3 legs of dried museum specimens using the Nucleospin Tissue Kit (Macherey-Nagel, Düren, Germany). From these specimens, we sequenced the mitochondrial COI gene and up to 5 nuclear genes (EF1a, GAPDH, MDH, RpS5 and wingless). Primers and laboratory protocols followed Wahlberg and Wheat (2008) in all cases. Successful PCR products were cleaned enzymatically with exonuclease and FastAP (Thermo Fisher Scientific, Hampton, NH, USA) and sent to Macrogen Europe (Amsterdam, the Netherlands) for Sanger sequencing. The resulting chromatograms were examined by eye using BioEdit (Hall 1999). All six genes are protein-coding, thus alignments were trivial. The aligned sequences were curated and managed in VoSeq (Peña and Malm 2012).
Phylogenetic analyses were carried out using IQ-TREE 1.6.10 (Nguyen et al. 2015) in a maximum likelihood framework. The full dataset was analysed. The data were partitioned by gene and analysed with the partition finding (Chernomor et al. 2016) and model finding (Kalyaanamoorthy et al. 2017) algorithms of IQ-TREE (using the command MFP+MERGE). Robustness of the results were assessed using UFBoot2 (Hoang et al. 2018) and a SH-like approximate likelihood ratio test (Guindon et al. 2010). Analyses were run on the CIPRES server (Miller et al. 2010).
Timing of divergence was estimated using a reduced dataset. One specimen per species with the most DNA sequence data available was chosen for this analysis. The dataset was analysed using BEAST v1.8.3 (Drummond et al. 2012), partitioned according to the results of the IQ-TREE analysis. Since BEAST does not offer the possibility to use the TIM2 model, the model for COI was set to the GTR+G model, the other two partitions were assigned models according to the IQ-TREE results (see below). The nucleotide models and clock models were unlinked between partitions. The tree prior was set to the birth-death model and a lognormal relaxed clock was implemented. Monophyly of Aporiina was enforced. Four secondary calibration points were taken from Chazot et al. (2019) and were modeled with a normal prior. The root was set to 42 million years ago (Mya) ± 5 million years (My), the crown age of Aporiina to 36 Mya ± 4 My, the most common recent ancestor of Mylothris and Aporia Hübner, 1819 to 28 Mya ± 3.5 My, and the most recent common ancestor of Catasticta and Melete Swainson, 1831 to 20 Mya ± 3 My. The analyses were run four times independently at 10 million generations per run sampling at 10,000 generation intervals. Convergence was examined using Tracer v1.7 (Rambaut et al. 2014). The four runs were combined after a burnin of 1 million generations and a summary tree was generated using the maximum a posteriori tree and mean estimated times of divergence.
The historical biogeography of the group was inferred using the R package BioGeoBEARS v1.1.1 (Matzke 2014). The chronogram resulting from the BEAST analysis (see results section) with the outgroups pruned, was used for this analysis. Owing to the open conceptual debate on model comparisons in BioGeoBEARS (Ree and Sanmartín 2018), our ancestral range analysis were performed under the dispersal -extinction -cladogenesis (DEC) model (Ree et al. 2005;Ree and Smith 2008). DEC is a likelihood-based procedure that models range evolution as a discrete-valued process, specifying instantaneous transition rates between geographical ranges along the branches of a phylogenetic tree. Consequently, DEC requires each sampled taxon be assigned with one or more geographic range(s), reflecting its present distributional range. The geographic distributions of Mylothris were extrapolated from the literature and the degree of endemism of local species, which is highest in West Africa and Madagascar. The following regions were used in the analysis: A, West Africa; B, West Central Africa; C, East Central Africa; D, East Africa; E, South Africa; F, Malagasy Region (see inserted map in Figure 6 of Results section). B and C are closer than the others, corresponding roughly to the central forest belt, but with two main centres of endemism: the Cameroon highlands and the montane zone just west of the Albertine Rift. In all dispersal rate scaler matrices, dispersal events between adjacent areas were not scaled (dispersal d = 1.0). A scaler of 0.5 was applied for dispersal events between areas separated by an ocean barrier. A scaler of 0.25 was applied for dispersal events between areas separated by a third region and a water barrier (see Suppl. material 2: Table S2). The maximum number of areas per ancestral state was fixed to four.

Results
The dataset consisted of 5168 aligned base pairs. The IQ-TREE partition finding analysis combined EF1a, GAPDH, MDH, and RpS5 into one partition and assigned the GTR+I+G model to it. COI was kept as its own partition with TIM2+I+G model assigned to it, as was wingless with the K2P+G model.
We find Mylothris to be a well supported monophyletic group (UFB = 100, SH-like = 100), that is sister to a clade of Neotropical and Palaearctic/Oriental taxa ( Figure 1). We identify five major clades that are strongly supported (the exception being the hilara clade with UFB = 87, SHlike = 90): the jacksoni, elodina, rhodope, agathina and hilara clades. We find the jacksoni clade to be sister to the rest of Mylothris, with elodina being the next lineage to diverge, followed by the rhodope clade, and finally the agathina and hilara clades being sister to each other.
Relationships Only M. elodina Talbot, 1944 is found in the elodina clade ( Figure 2). The eight individuals we sequenced show very little genetic differentiation, but a relatively long branch leads to them. The position of the lineage is not highly supported (UFB = 79, SH-like = 85.3), and there is a possibility that further data might place it as the sister to the rest of Mylothris.

Western Africa (A)
West Central Africa (B)
The hilara clade comprises 15 sampled species, and again similar patterns of little genetic differentiation are found among the species as in the other clades ( Figure 5). The species M. sjostedti Aurivillius, 1895, M. similis Lathy, 1906, M. poppea Cramer, 1777and M. erlangeri Pagenstecher, 1902 appear to be distinct from each other and the rest, while the rest of the species form a complex set of relationships with shared haplotypes and minor monophyletic groups (e.g. M. nsp1), which are very closely related genetically, even if morphological characters are clearly distinct.
Our timing of divergence analysis suggests that the crown clade of Mylothris began diversifying around 16 Mya ( Figure 6) after diverging from their sister group some 10 My earlier. The five major clades diverged from each other between 16 and 8 Mya, during the mid to late Miocene, but almost all species level divergences have happened in the Pleistocene. Our ancestral range analysis recovers an origin of Mylothris in the East Central Africa. All the five major clades of the group are also inferred to have been derived from ancestors distributed in the East Central Africa. The ancestors of the extant species in Madagascar and Comoros islands are estimated to have range expanded from mainland Africa only recently (~2 Mya).

Discussion
Mylothris are a distinct and well-supported lineage of African pierids, not closely related to any other pierids on that continent. There are five clear species groups within the genus. Of the five groups, the rhodope, agathina and hilara species groups are more closely related to each other than to the other two groups. Based on our sampling of genes, it appears that the jacksoni group is sister to the rest, but with poor support. Indeed, in some preliminary analyses of our data, the elodina group came out as sister to the rest, and this may be the case with further gene and taxon sampling. The very different morphology of elodina, reflected also in the male genitalia, is a possible indicator that it is sister to the rest of Mylothris.
It is highly probable that we have uncovered all major species groups. The unsampled species are all, as far as we can tell, clearly morphologically related to already sampled species. The sampling process was conducted in the context of a revision of the genus (Warren-Gash in press), and large quantities of material were examined and dissected, including many specimens too old to be readily amenable to molecular sequencing. The unsampled species showed no differences which were important enough in terms of morphology or genitalia to suggest the likely existence of other species groups.
Regardless of the relationships of the major species groups, the biogeographic history of the deeper clades appears to be clear. We infer from our analyses that the group originated in East Central Africa, and from there spread to other parts of Africa. The fragmentation of the forests in the Miocene and Pleistocene may have contributed to the diversification of the group, perhaps by causing populations to be fragmented and isolated from each other. This appears to be a common pattern in forest dwelling butterflies of Africa (Aduse-Poku et al. 2009;van Velzen et al. 2013;Dhungel and Wahlberg 2018), but the mechanisms leading to their diversification are still poorly understood.
Within the species groups of Mylothris we find genetic variation between some of the species to be minimal, and indeed some species share the same haplotypes for the markers sequenced in this study. The infra-specific variation within Mylothris species is, however, considerable, and that has been given full weight both in the preparation of this paper and the forthcoming revision. A number of morphs previously considered as species have been synonymised in consequence (Warren-Gash in press). Nonetheless, such close genetic similarity needs explanation if the specific separation of these groups is to be justified.
There is no single answer to why some species are genetically identical for the markers we have sequenced. In one case -Mylothris agathina and M. chloris -they are sympatric with no evidence of interbreeding, and in addition the female genitalia of the two are clearly and consistently different. In the M. jacksoni group, the wing patterns vary, as does the biome in which they fly. Many of the scarcer species, both sampled and unsampled, are montane vicariants; others fly sympatrically without interbreeding and in several cases where the early stages have been observed, the larvae are more distinct than the adults and use different larval foodplants (e.g. M. jacksoni and M. sagala). They not only look different but in their foodplant preferences behave differently.
In other groups, there are examples of species which fly together in at least a part of their range, which are morphologically similar in one sex (usually the male) but distinctively different in the other (e.g. M. sulphurea basalis Aurivillius, 1906 andM. aurantiaca Rebel, 1914, with overlapping ranges in the north-east of the DRC).
A group of more than a dozen species in the hilara group, which have in common a red basal patch at the base of the forewing and broadly similar (but not identical) male genitalia, is found across the width of the main forest belt. A few are very local. In the rugged terrain of the mountains west of the Albertine Rift, they are separated more by terrain than geographical distance, but in those few instances where they overlap, they do not interbreed. At lower altitudes it is quite usual for several species in the group to fly together in varying combinations, without any evidence of cross-fertilisation. In the Albertine Rift widely varying female forms are found which occur nowhere else. Morphology shows up what appear to be widely separated vicariants of similar species in east and west Africa, more akin to each other than those flying in between. They almost certainly use different larval foodplants. Some are tied to strictly forest environments (e.g. M. continua); others prefer more open country (e.g. M. rueppellii Koch, 1865) The case of the jaopura complex, with nine named and mostly allopatric species, offers a different challenge. Those sequenced are indistinguishable on that basis, and the genitalia, while distinctive within the complex, are very similar. However, morphology in terms of size, wingshape and markings is distinctive in each case, as very often is the biome in which they occur. Where they overlap (e.g. M. kiwuensis and M. nsp2), they do not interbreed. Again, and given the very different environments in which they occur, it is highly probable that their larval foodplants also differ. What they have in common is that in practically every case they mimic a Mylothris species in another group, sometimes so well that they can only be separated with certainty by dissection or sequencing.
Given the factors discussed above, the markers we have used are probably not sufficiently variable to differentiate the closely related species. It is likely that a genomic approach would shed further light on the species boundaries, as has been found in recent studies on nymphalid butterflies (Dinca et al. 2019;Campbell et al. 2020). Our results do however suggest that radiation within these groups is likely to be recent. Similar minimal genetic differentiation is found in the pierid genus Colias Fabricius, 1807 (Wheat and Watt 2008; Laiho and Ståhls 2013), but interestingly not in the closely related genus Delias (Braby et al. 2007;Müller et al. 2013). The reasons behind the close genetic relationships of several groups of species is not entirely clear, but may have to do with the repeated contraction and expansion of the forests in Africa, especially during the Miocene and Pleistocene. Similar patterns appear to be the case in other forest specialist genera of the family Nymphalidae: Charaxes Ochsenheimer, 1816 (Aduse-Poku et al. 2009), Cymothoe Hübner, 1819(van Velzen et al. 2013 and Euphaedra Hübner, 1819 (Dhungel and Wahlberg 2018). Mylothris and the nymphalid genera cited appear to be going through rapid speciation in African forests and they require further study to understand the driving mechanisms