WIntaerern, Kaotiholin, aMl eJonduornzaa,l .o.. f& O Sduohnliangtology Evidence for widespread gene flow and migration in Pantala flavescens 2022, Vol. 25, pp. 43–55 doi:10.48156/1388.2022.1917166 Evidence for widespread gene flow and migration in the Globe Skimmer dragonfly Pantala flavescens Jessica Ware 1*, Manpreet Kaur Kohli 1, Ciara Mae Mendoza2, Daniel Troast3, Hiroshi Jinguji 4, Keith A. Hobson 5, Göran Sahlén 6, R. Charles Anderson7 & Frank Suhling 8 1 Division of Invertebrate Zoology, American Museum of Natural History, USA 2 The Hansjörg Wyss Department of Plastic Surgery, New York University School of Medicine, New York, USA 3 Stem Cell Initiative Flow Cytometry Core, Columbia University, USA 4 Faculty of Food and Agricultural Sciences, Fukushima University, Japan 5 Department of Biology, University of Western Ontario, Canada 6 Rydberg Laboratory for Applied Sciences, Halmstad University, Halmstad, Sweden 7 Manta Marine, PO Box 2074, Malé, Maldives 8 Institute of Geoecology, Working Group Landscape Ecology & Environmental Systems Analysis, Technische Universität Braunschweig, Germany *Corresponding author. Email: [email protected] Abstract. The global population structure and dispersal patterns of Pantala flavescens (Fa- bricius, 1798) are evaluated using a geographically extensive mitochondrial DNA dataset, a Research Article more limited samples of nuclear markers, wing isotopic (δ²H) data and a literature review. No spatial or temporal haplotype structure was recovered between the samples. Isotope OPEN ACCESS data suggest that most samples were immigrants at the collection locations. A literature This article is distributed review of migration events for the species confirms regular inter-and intra-continental mi- under the terms of the grations occur (the majority reported from Asia, Africa and Australasia), with individuals Creative Commons and swarms dispersing thousands of kilometers over land and oceans. Migrations coincide Attribution License, with prevailing winds and seasonal rains, which points to a mechanism we name the “pan- which permits unrestricted use, tropical Pantala conveyor belt”, suggesting widespread gene flow is possible for an aquatic distribution, and reproduction in insect with excellent flying ability linked to rapid larval development. any medium, provided the original author and source are credited. Key words. Odonata, deuterium, haplotype, isoscape, F , migration, ΦPT ST Published: 11 March 2022 Received: 29 July 2021 Accepted: 3 March 2022 Introduction Citation: The Globe Skimmer or Wandering Glider dragonfly, Pantala flavescens (Fabricius, Ware, Kohli, Mendoza, Troast, 1798) (Odonata: Anisoptera: Libellulidae), is a well-known long-distance migrant. Jinguji, Hobson, Sahlén, Anderson It has a circum-global distribution, and is the most widely distributed of all dragon- & Suhling (2022): fly species (Fraser, 1936; Russell et al., 1998). Pantala flavescens occurs most often Evidence for widespread gene flow in warmer climates, especially within the tropics, but also migrates seasonally into and migration in the Globe Skimmer dragonfly Pantala flaves- temperate zones (Borisov & Malikova, 2019; May, 2013). Its ability to disperse over cens. International Journal of great distances is demonstrated by its occurrence on offshore islands, for example Odonatology, 25, on many Pacific Ocean islands (Rowe, 2004; Schmidt, 1938) including New Zealand 43–55 (Corbet, 1979; Lieftinck, 1975; Rowe, 1980). Pantala flavescens has even reached doi:10.48156/1388.2022.1917166 the remote Sub-Antarctic Amsterdam Island (Devaud & Lebouvier, 2019) and is established on Rapa Nui (formerly Easter Island) (Dumont & Verschuren, 1991; Data Availability Statement: Moore, 1993; Samways & Osborn, 1998). Its longest, regular, well-documented All relevant data are within the paper and its migratory route is across the western Indian Ocean (Anderson, 2009; Hobson et Supporting Information files. al., 2012a), where it has been suggested to follow seasonal (monsoonal) rains, International Journal of Odonatology │ Volume 25 │ pp. 43–55 43 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens which allow females to take advantage of temporary leading to panmixia (Als et al., 2011; Palm et al., 2009; pools of water for oviposition (Anderson, 2009; Corbet, but see Wirth et al., 2001 for counter argument). Wide- 1999; Suhling et al., 2015). spread genetic mixing has also been recorded in several The ability of P. flavescens to migrate such long dis- other groups of organisms, including insects, birds, bats tances is facilitated by a unique combination of morpho- and plants (Als et al., 2001; Naro-Maciel, 2011; Neeth- logical characteristics. These include distinct wing mor- ling et al., 2008; Oomen et al., 2011; Ridgway et al., phology (Alvial et al., 2019; Li et al., 2014; Moore, 1993; 2001; Peel et al., 2013; White et al., 2011). Neverthe- Sacchi & Hardersen, 2013; Suárez-Tovar & Sarmiento, less, most reports of panmixia might best be described 2016; Outomuro & Johansson, 2019; Zhao, 2012), of as ‘local panmixia’, where random mating occurs across which the most obvious adaptation is their enlarged a population in a particular region rather than across hind wing bases (a feature associated with gliding, the entire species (e.g. Peel et al., 2013; White et al., Corbet, 1999), as well as unique thoracic musculature 2011). Studies on the population structure of long- (Bäumler et al., 2018). Physiological features that allow range migrant insects are still rare, with the exception individuals to migrate long distances likely involve their of locusts (Locusta migratoria, (Linnaeus, 1758) and ability to make efficient use of lipids as fuel rather than Schistocerca gregaria Forsskål, 1775), monarch butter- carbohydrates (e.g., Kallapur & George, 1973); it has flies (Danaus plexippus (Linnaeus, 1758)), painted lady been suggested that dragonfly fat reserves can sustain butterflies (Vanessa cardui (Linnaeus, 1758)), and the eight hours of flapping flight. In addition, P. flavescens green darner dragonfly, Anax junius (Drury, 1773) (e.g., may be able to feed on insects while migrating thereby Hallworth et al., 2018). replenishing their fat stores en route (Anderson, 2009). Here, we evaluate the global population structure of Given its wide distribution and extraordinary migra- P. flavescens, using genetic data to test the hypothesis of tory abilities, there is interest in the population struc- mixing across its global range. We have compiled genetic ture of P. flavescens. Studies in single countries have sequences available through the NCBI GenBank online revealed high haplotype diversity and an expected lack database, supplementing this with a limited number of of population structure (Cao et al., 2015; Hayashi & new samples, to form the largest such dataset yet as- Arai, 2004; Low, 2017). Over a much wider geographi- sembled, including samples from 29 countries across six cal range, albeit with a small sample size, Troast et al. continents. We estimated various population structure (2016) found that P. flavescens specimens from South measures, including testing for levels of gene flow. Fur- America, North America, and Asia shared mitochondrial ther, we measured stable hydrogen isotope (δ²H) values DNA haplotypes, suggesting high rates of gene flow and in a small sample of P. flavescens wings to establish what a single global panmictic population. Similarly, Alvial et proportion of these specimens originated in regions dif- al. (2017) found a lack of population structure between ferent from those where they were collected (Hayashi samples from across South and Central America using & Arai, 2004; Hobson et al., 2012a, 2012b; 2021). This both mitochondrial and nuclear DNA markers. In con- allowed us to test for the degree of migratory behaviour trast, a subsequent study by Alvial et al. (2019), which in P. flavescens, on the assumption that the presence of compared samples not only from Central and South a large number of immigrant individuals across a global America but also from islands in the Pacific and Indian range would confirm widespread migration. Finally, we Oceans, did find some suggestion of population struc- compiled published observations of migratory behav- ture between ‘continental’ and ‘insular’ samples using iour which served to identify possible connections and microsatellite loci, which also suggested that individuals pathways between and within the continents. from Rapa Nui formed a genetically distinct population. In addition, Pfeiler & Markow (2017) reanalysed pub- lished data (Low, 2017; Troast et al. 2016) and suggest- Materials and methods ed that ‘concluding global panmixia in P. flavescens may Study organism be premature.’ Thus, the results from previous studies are ambiguous, with some suggesting panmixia while Pantala flavescens (Fabricius, 1798) is a member of the others suggest some degree of population structure. dragonfly family Libellulidae (Odonata: Anisoptera). Global panmixia, or random mating and uninterrupt- The vernacular English names “globe skimmer” and ed gene flow across an entire species, appears to be un- “wandering glider” reflect this species’ long-distance common in widespread species. Barriers to genetic ex- flying ability and wide distribution. Adults are red, yel- change, which may be geographical, temporal or behav- low, or brown in colour, commonly observed “hawking” ioural in nature, tend to lead to the evolution of geneti- over open areas. Egg-laying often occurs in tandem, cally distinct regional populations. However, in some and females produce their eggs continuously (Koch et species dispersal or migration can lead to widespread al., 2011; Sharma, 2017; Ware et al., 2012). Eggs are gene flow despite the presence of obstacles (Neeth- distributed widely, even being oviposited among differ- ling et al., 2008; Oomen et al. 2011). A commonly cited ent ponds during a single oviposition event, which is in- example is that of the European Eel (Anguilla anguilla terpreted as risk-spreading to avoid loss of all offspring (Linnaeus, 1758)), in which individuals from separate should temporary habitats dry out before development river basins travel to the Sargasso Sea for reproduction, is completed (Schenk et al., 2004). Larvae mostly inhab- International Journal of Odonatology │ Volume 25 │ pp. 43–55 44 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens it temporary or ephemeral lentic waters including rain- We additionally reconstructed haplotype networks pools forming after heavy rainfalls. Offspring develop for P. flavescens sequences available from NCBI Gen- very rapidly, taking 4–6 weeks from egg hatch to adult Bank for mitochondrial 16S and Cytochrome B, as well emergence (e.g., Hawking & Ingram, 1994; Kumar, as nuclear 18S and Histone 3. Apart from 16S, these 1984; Suhling et al., 2004). sequence datasets were from single countries, and so although they do provide insights into genetic diversity they shed no direct light on global population struc- Sampling ture. We collected 23 specimens of P. flavescens from six countries during 2004 to 2015 (see Supplementary ma- DNA Extraction/PCR Amplification terial, Table S1). In addition we downloaded genetic se- quences of P. flavescens from NCBI GenBank (N = > 600 DNA was extracted from tissue of P. flavescens samples for COI; complete list in Supplementary material, Table which had been stored in ethanol or acetone-dried. S1). We included COI sequences from North America Extraction was completed using Qiagen DNeasy Blood (N = 16; Canada, USA), Central and South America (N = and Tissue Kits (Maryland, USA). We sequenced the 25; Brazil, Chile, Costa Rica, Guyana, Peru), Asia (N = mitochondrial Cytochrome Oxidase I (COI; ~ 658 base 190; Cambodia, China, India, Indonesia, Japan, Korea, pairs), a gene which was chosen because of a large Malaysia, Maldives, Pakistan, Philippines, Singapore, number of existing sequences available on NCBI, allow- Vietnam), Africa (N = 4; Guinea Bissau, Liberia, Sen- ing for this large, global dataset. Polymerase Chain Re- egal; unfortunately, samples from South Africa and action (PCR) was done using the primer set of Simon Namibia failed during the sequencing step), Australia et al. (1994): 1709f (TAATTGGAGGATTTGGAAATTG, T = m and the Pacific (N = 13; Rapa Nui [a territory of Chile], 56.2°C) and 2191r (CCYGGTARAATTARAATRTARACTTC, Fiji, French Polynesia, Hawaii, Mariana Islands, Tonga T = 56.7°C). PCR reaction total volume was 25 µL, com- m Islands). Additionally, we downloaded and ran analy- prising 12.5 µL Taq 2X Master Mix, 5–7 µL DNA, 3.5– ses on P. flavesc ens nuclear genetic data available from 5.5 µL RNA-ase free water, and 1–2 µL of the forward NCBI (Table S1). and reverse primers. PCR was done using an Eppendorf Thermocycler us- ing the following parameters: initial heat step at 94.0°C Genetic dataset partitions held for one minute, 20 cycles of denaturation (30 s, 94.0°C), annealing (45 s, 48.0°C), and extension (45 s, We ran our analyses with the geographic dataset de- 72.0°C), an additional 20 cycles of denaturation (45 s, scribed above for the COI. Additionally, we ran the COI 94.0°C), annealing (45 s, 50.0°C), and extension (45 s, haplotype analyses with either (1) fewer Japanese sam- 72.0°C), followed by a final extension at 72.0°C held for ples or (2) excluding samples from Alvial et al. (2017, five minutes. PCR products were visualized via gel elec- 2019). trophoresis. We used commercially available Macrogen 1. Inference of population structure can be sensitive to services (New York City, New York, USA) for sequencing sampling bias, (see Puechmaille, 2016). More than and purification. half of our samples came from Japan, while other re- gions had < 10 individuals. Thus, for additional hap- lotype analyses we limited the number of Japanese Sequences samples to 77 (“reduced” data set). We chose these samples to sample P. flavescens evenly from across Sequences were edited using Sequencher® (2016). the Japanese islands for which we had data. This These sequences, along with those from NCBI Gen- gave us 248 P. flaves cens sequences in total. Bank, were aligned using Clustal-X Version 2.1 (Larkin, 2. For NCBI samples from Alvial et al.’s (2017, 2019 2007) and Mesquite Version 2.75 (Maddison & Maddi- studies (N = 40)) it was impossible to distinguish son, 2017). For the COI, twenty-three de novo sequenc- which samples were which. Hence, these were es from Africa, Australia, the Maldives, and the USA coded as “Alvial” (i.e. grouping together all of her were sequenced and resulting sequences were depos- samples from Costa Rica, Chile, Peru, Tonga, Rapa ited in NCBI GenBank (Accession numbers OM945847– Nui, and the Maldives). Alvial et al.’s samples were OM945869). In total we had 25 samples from Central included in testing our Null hypothesis that there is and South America, 16 from North America, 13 from no global geographic structure. These samples had Australasia, 190 from Asia including the Maldives, and to be excluded from the rest of our geographic hy- four from Africa for the COI. For the 16S the 16 sequenc- pothesis testing scenarios as we could not infer prov- es were from Africa, India, Japan, Korea, Malaysia, Paki- enance (i.e., if the samples were from the Northern stan, and the United States of America, for the 18S the or Southern, Eastern, or Western Hemisphere, see four sequences were from Japan, for the Cytochrome B Inference of Population Structure section for more the 77 sequences were from China, and for the Histone details on scenarios). 3 the four sequences were from Japan. International Journal of Odonatology │ Volume 25 │ pp. 43–55 45 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens Haplotype analyses tion rates (4N M) and Ө = 4N μ under scenarios (c) and e e (d). We let the migration matrix stay at default settings. To assess variation among individuals of P. flavescens, We ran 10 replicate maximum likelihood Migrate-N we reconstructed minimum spanning networks for runs of 20,000 steps per chain each with a Brownian each gene fragment using PopArt version 1.7 (Leigh & motion model, and one long chain. We sampled every Bryant, 2015) with the minimum spanning networks al- 100 steps, with a burn-in of 10,000 steps and a bound- gorithm (Bandelt et al., 1999). ed-adaptive heating scheme with four temperatures (3 cold, 1 hot) of 1.00, 1.50, 3.00, and 10,000. To test the effect of temperature on our results, we ran an Inference of population structure analysis with two different “hottest temperatures” and the results did not differ. We evaluated the best model We also assessed population structure summary sta- using Log Bayes Factor (LBF) and highest model prob- tistics using the COI genetic data. All statistics were ability values in Migrate-N. Effective Sample Sizes (ESS) estimated from the “reduced dataset”, and required were evaluated in TRACER version 1.7.1 (Rambaut et assigning individuals to a fixed number of putative geo- al., 2018). graphically based hypotheses (populations), which we did in several ways. First, we estimated the number of genetic clusters/populations using K-means clustering Stable Isotope Analysis analysis in adegenet R package 2.1.3 (Jombart, 2008; Jombart & Ahmed 2011). We selected the optimal num- We excised wings from 29 individuals representing ber of clusters using the Bayesian Information Criterion each of the continents of North America, Central and (BIC); choosing the solution with the best BIC value. Dis- South America, Africa, Australia, and Asia plus the criminant Analysis of Principal Components (DAPC), as Pacific; samples were chosen from both islands and implemented in adegenet R package, was then used to mainland areas (Table S2). All samples had been dried assess genetic variation between identified population in acetone, which has been shown to have no impact clusters. We also determined membership probability on isotopic results (Hobson et al., 2012a; Hobson & of each sample for each of the predicted populations; Wassenaar, 2019). We submitted wings to the Cornell to do so we calculated the degree of genetic differentia- Isotope Laboratory (COIL) facility, where surface oils tion (F ) among predicted clusters using the hierfstat R ST were removed using a 2:1 chloroform:methanol solu- package version 0.5-7 (Goudet, 2005), and we ran an tion and air dried in a fume hood. Subsamples of wings AMOVA in the poppr package version 2.8.6 (Kamvar et were weighed into silver capsules, crushed, and loaded al., 2014, 2015) to test whether the observed variation into a zero-blank carousel attached to TC/EA (Thermo, was higher between or within each cluster. Bremen, Germany) where samples were combusted at We estimated population genetics statistics assuming 1350°C pyrolytically. Resultant H gas was analyzed un- different scenarios considering possible biogeographi- 2 der continuous-flow isotope-ratio mass spectrometry cally separated populations: (a) assuming the popula- (CFIRMS) on a Thermo Delta-V mass spectrometer us- tions were defined based on the results from K-means ing a helium carrier gas. Samples were calibrated us- clustering, (b) assuming each of the continents is a sepa- ing within-run measurements of keratin standards CBS rate population, (c) assuming Northern/Southern Hemi- (-197‰) and KHS (-54.1‰) according to the compara- spheres are separate populations, (d) assuming Eastern/ tive equilibration method of (Wassenaar & Hobson, Western Hemispheres are separate populations, and (e) 2003). An in-house (Spectrum) commercial standard a null hypothesis of one panmictic population. Insects was used to correct for instrument drift. Stable-hydro- generally do not follow geo-political boundaries, but we gen isotope values are reported here in standard delta also ran an additional scenario (f) assuming countries notation in parts per thousand (‰) deviation from the are separate populations. Population genetics statistics Vienna Mean Ocean Water (VSMOW) standard. Based under the various scenarios were estimated in the soft- on within-run replicate measurements of standards, we ware package R as described above. We calculated the estimated measurement precision to be ~ 2‰. genetic structure index, F and other statistics (haplo- ST, Assigning P. flavescens to geographical origins using type diversity, Hd; nucleotide differences, K; number of the stable-hydrogen isotope approach is problematic polymorphic sites, S; Tajima’s D, Fu and Li’s statistics, because this species is known to use ephemeral pools and nucleotide diversity, π) under each of these popula- that may be created from highly seasonal (often mon- tion scenarios as described above. For the single pan- soonal) rainfall events. Hobson et al. (2012a) modeled mictic population global scenario we could not estimate the origins of individuals captured on the Maldives en F as there was only one population. ST route to Africa using an amount-weighted mean annual precipitation model that provided the first isotopic evi- Effective migration rates dence for long-distance migration in this species and pointed to origins in northern India or possibly even For the COI, we used the program Migrate-N version further north and east. That approach was based on a 3.7.2 (Beerli et al., 2019) to estimate effective migra- calibration algorithm derived from known-origin drag- International Journal of Odonatology │ Volume 25 │ pp. 43–55 46 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens onflies in North America (Hobson et al., 2012b: Wing ern/Southern Hemisphere), we recovered the same δ²H = 0.91(MAD) -42.5, r² = 0.75, where MAD is the haplotype structure. For the other gene fragments stud- amount-weighted mean annual precipitation for a giv- ied, one or two main haplotypes were recovered (Fig- en site). Based on the findings of Lopez-Calderon et al. ure S3), with the exception of Cytochrome B, for which (2019), we used this relationship to evaluate whether several haplotypes were suggested among the samples a given dragonfly likely originated from sites of capture (all of which were from China). In the case of the mi- using the normal probability density function and with tochondrial marker 16S, we located 16 samples from an assumed standard deviation of the residuals of the three continents on NCBI GenBank, all of which shared calibration algorithm to be a highly conservative 19‰ a single haplotype (Figure S3). In the case of the nuclear (Hobson et al., 2012b). We considered three odds ratios marker H3, only four sequences were available on NCBI (2:1, 3:1, and 4:1) for each individual assignment as be- Genbank, three of which shared the same haplotype; ing a “local” or an “immigrant”. the fourth sample revealed a very different haplotype, and likely represents a different species, mislabeled as P. flavescens, or contamination. Results Haplotypes Population Genetics Our COI dataset consisted of an alignment with 686 characters. Haplotype minimum spanning analyses on The clustering method in the adegenet package esti- all datasets suggest that most P. flavescens samples mated K = 1 to three populations of P. flavescens as most share a single common haplotype for COI (Figure 1, net- likely for COI. We used the “elbow” method to choose work from reduced dataset); we found one main hap- K = 2 and K = 3 as the population cluster values (sce- lotype with 24–43 additional haplotypes varying by ≤ 3 nario (a); Jombart & Collins, http://adegenet.r-forge.r- nucleotides (43 = full dataset and 24 = reduced dataset). project.org/files/tutorial-genomics.pdf). However, as Regardless of what level of population assignment was inferred clusters higher than K ≥ 2 lacked any geograph- used (continent, Western/Eastern Hemisphere, North- ic structure, we consider these values as to be difficult Figure 1. PopART minimum spanning network; Dataset: Alvial’s samples not included, reduced number of Japanese samples. Haplotype networks coloured based on sampling country, continent, hemisphere and New/Old world. International Journal of Odonatology │ Volume 25 │ pp. 43–55 47 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens to interpret: they do not seem to share any clear pat- from West Africa, this result should be treated with tern that unites the individuals in these clusters (i.e., caution. For scenario (c), F was 0.007; for scenario (d) ST these clusters do not reflect timing of sampling, geo- F was 0.004. In short, for the geographical scenarios ST graphic region of sampling, etc.). In addition, the BIC (b)–(d) our dataset suggests high gene flow and no ge- scores did not vary greatly among the K = 1, K = 2, and netic structure between putative populations. K = 3 (BIC score range 300–~450) and DAPC scatterplots showed large overlap across clusters (Figure 2). We fur- ther tested whether differences between K = 2 and K = Effective migration rates 3 could be explained by temporal bias in sampling, but as each putative K population contained samples from Migrate-N found widespread gene flow in COI. Mi- across sampling periods (i.e., early 2000s and early/mid grate-N analytical results suggested that individuals 2010s), we reject this explanation for inferred patterns. migrating (in either direction) from Western Hemi- Lastly, under the K = 2 and K = 3 scenarios we found low sphere to Eastern Hemisphere was the best fit with genetic structure and F values between the popula- ln(Prob(D|Model)) = -4158.48 compared to the model ST tions, consistent with high rates of gene flow (most F for individuals migrating (in either direction) from the ST values were < 0.5). Northern/Southern Hemisphere (Northern/Southern F values between different biogeographical sce- Hemisphere ln(Prob(D|Model)) = -6480.57). Within ST narios varied, ranging from below 0.05 in scenario (b) the Eastern/Western Hemisphere analysis, movement (between Asia and Africa, Asia and Australasia, North in both directions had normal distributions of migra- America and Asia) indicating little population differen- tion rates, M, although relative variance around the tiation, to 0.55 between Africa and Australasia (Table 1) means varied. Mean migration rate for the Western suggesting greater population differentiation. However, Hemisphere-Eastern Hemisphere direction was 478.2, with only four samples from Africa, all of which were and the mean migration rate in the Eastern Hemisphere Figure 2. Density and cluster plots for defined populations; (A) density plot assuming continents are populations; (B) BIC val- ues for the putative number of K clusters, (C) cluster plot assuming countries are populations; In each cluster plot and density plot putative populations are colour-coded. International Journal of Odonatology │ Volume 25 │ pp. 43–55 48 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens Table 1. Pairwise F values under various population scenarios: (a) the K-means clustering analysis populations estimated in ST adegenet R package, (b) assuming continents are separate populations. Scenario (a) K=2 K=3 Population F 2 1 3 2 1 ST 3 – – 0.0 0.217 0.158 2 0 0.168 0.217 0 0.277 1 0.168 0 0.158 0.277 0 Within population variation 58.983 56.462 Among population variation 41.017 43.538 Scenario (b) Population F (b) Africa Asia Australasia N America S America ST Africa (N=4) 0 0.029 0.551 0.197 0.223 Asia (N=190) 0.029 0 0.002 0.031 0.062 Australasia (N=13) 0.551 0.002 0 0.158 0.111 North America (N=16) 0.197 0.031 0.158 0 0.131 South America (N=25) 0.223 0.062 0.111 0.131 0 Within population variation 95.889 Among population variation 4.111 Table 2. Population genetics statistics based on different population scenarios, with K: nucleotide differences, π: nucleotide diversity, S: # of polymorphic sites. p ≤ 0.05 in bold. Population Tajima’s D K π S Hd variance Fu and Li’s D Fu and Li’s F Fst Western Hemisphere -1.37923 0.300 0.001 3 0.236 0.009 -1.50777 -1.70431 0.004 (N=41) p>0.10 p>0.10 p>0.10 Eastern Hemisphere -2.18914 0.943 0.006 24 0.573 0.002 -3.24440 -3.38651 (N=207) p<0.01 p<0.05 p<0.02 Northern Hemisphere -2.16518 0.494 0.003 20 0.581 0.002 -3.69480 -3.69542 0.007 (N=228) p<0.01 p<0.02 p<0.02 Southern Hemisphere -1.79423 0.251 0.002 3 0.328 0.009 -1.64715 -1.97431 (N=20) 0.10>p>0.05 p>0.10 0.10>p>0.05 All data -2.35627 0.452 0.003 22 0.348 0.002 -4.45510 -4.34133 p<0.01 p<0.02 p<0.02 to Western Hemisphere direction was 40.0; the mean individuals had lower wing δ²H values than predicted migration rate in the North-South direction was 311.1, considering mean amount-weighted local precipitation and the mean migration rate in the South-North direc- δ²H (Figure S1), implying more temperate regions of or- tion was 19.5. In short, the results of this analysis sug- igin. In only three cases was the individual wing δ²H val- gest more migration going from west to east and north ue higher than expected from the local collection site. to south than in the opposite directions. Using an odds ratio of 2:1 or 3:1 resulted in assignment of at least 20 out of 29 individuals as immigrants vs. lo- cally produced (Table 3). This changed only slightly with Stable Isotope Analyses a more conservative 4:1 odds ratio (19 out of 29). Note that these estimates represent minimum numbers of Our isotope data indicated a broad span of isotopic ori- immigrants; specimens with δ²H values consistent with gins for our samples, ranging from xeric to mesic habi- local origin could have originated from a different area tats and different geographic origins (Figure S1). Most with a similar local precipitation δ²H. International Journal of Odonatology │ Volume 25 │ pp. 43–55 49 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens Table 3. Stable isotope (δ²H) results of wing chitin of individual Pantala flavescens taken from subsamples of collections across at the sites as indicated. Also indicated are the assignment to immigrant (im) or local (loc) based on given odds ratios as de- scribed in the Methods. ID Location Latitude Longitude Elevation Wing Precip. Predicted Odds Ratio (m) δ²H (‰) δ²H wing δ²H 2:1 3:1 4:1 (‰) (‰) 1 India 13.5662 78.499 0 -96.0 -58.7 -58.01 im im im 2 China 30.4919 119.615 1506 -121.7 -87.0 -108.06 loc loc loc 3 China 30.4919 119.615 1506 -94.4 -57.0 -108.06 loc loc loc 4 Japan 35.8896 139.652 12 -145.5 -113.1 -93.50 im im im 5 Japan 38.5587 140.85 208 -138.4 -105.3 -97.14 im im im 6 Australia -18.9 145.77 0 -141.9 -109.2 -55.28 im im im 7 Australia -18.9 145.77 0 -120.1 -85.2 -55.28 im im im 8 Australia -18.9 145.77 0 -70.6 -30.8 -55.28 loc loc loc 9 Australia -18.9 145.701 0 -83.2 -44.7 -55.28 im im im 10 Senegal 14.6433 -12.3363 23 -81.0 -42.3 -67.11 loc loc loc 11 Senegal 14.6433 -12.3363 23 -100.0 -63.1 -67.11 im im im 12 Senegal 13.9726 -15.0079 55 -107.5 -71.4 -68.93 im im im 13 Guinea-Bissau 11.683 -14.7977 14 -92.7 -55.1 -69.84 im im loc 14 Cruise ship 18.0294 64.1041 0 -89.7 -51.8 -58.01 im im im 18 Guyana 6.4833 -58.2167 0 -126.6 -92.4 -49.82 im im im 19 Guyana 6.4833 -58.2167 0 -108.5 -72.5 -49.82 im im im 20 Hawaii 21.6403 -158.063 5 -76.7 -37.5 -48.91 im im im 21 Hawaii 21.6403 -158.063 5 -62.8 -22.3 -48.91 loc loc loc 22 Hawaii 21.6403 -158.063 5 -65.5 -25.2 -48.91 loc loc loc 23 Hawaii 21.6403 -158.063 5 -77.1 -38 -48.91 im im im 24 Peru -7.449 -73.9403 295 -88.6 -50.6 -76.21 loc loc loc 25 Peru -9.8933 -76.32 1981 -96.7 -59.5 -108.06 loc loc loc 26 Peru -9.8933 -76.32 1981 -89.9 -52.0 -108.06 loc loc loc 27 Brazil -12.196 -38.9414 250 -137.3 -104.1 -48.00 im im im 28 Costa Rica 8.6143 -83.5411 304 -115.5 -80.2 -85.31 im im im 29 USA 35.96 -83.924 260 -130.9 -97.1 -73.48 im im im 30 USA 30.4642 -84.5953 43 -105.3 -69.0 -58.92 im im im 31 USA 30.3356 -97.9112 175 -89.2 -51.3 -61.65 im im im 32 USA 39.8683 -75.1877 4 -116.3 -81.1 -78.94 im im im Discussion across five continents for COI and all individuals collect- ed across three continents for 16S. Secondly, F values ST Our genetic data suggest a lack of population structure (i.e. estimates of genetic structure within and between across the global range of P. flavescens, which can be populations) for COI suggest little population differen- attributed to widespread migration. This interpretation tiation and high levels of gene flow between Northern/ is supported by our isotope data and by our literature Southern Hemispheres, Eastern/Western Hemispheres, review. These findings suggest the existence of a mech- and most continents (Table 2); country dataset results anism we call the pantropical Pantala conveyor belt. were not biologically meaningful, and are likely unre- Here, we have compiled the largest P. flavescens hap- liable given small sample sizes for some countries (cf lotype dataset to date, including for the first time indi- Excoffier, 2007). Thirdly, from our stable isotope analy- viduals from each of the five continents where this spe- sis it is clear that a majority of the individuals sampled cies is usually found. Our results demonstrate genetic were immigrants to their locations of capture, further homogeneity across all geographic regions, suggesting supporting widespread mixing of P. flavescens. either a recent shared common ancestry or widespread Nevertheless, incomplete lineage sorting, or recent dispersal. With this globally sampled dataset we find shared common ancestry, may also be associated with not only a lack of genetic structuring among regions, such low F values. In support of this possibility, Tajima’s ST but also high levels of gene flow, and this suggests D values (which estimate deviation from neutral expec- some level of panmixia. First and foremost, one main tations) were negative for all areas for COI, suggesting haplotype is shared among most individuals collected population expansion. These values were statistically International Journal of Odonatology │ Volume 25 │ pp. 43–55 50 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens significant in the Eastern and Northern hemispheres, We also note that among the individuals studied suggesting less polymorphism than expected under by Alvial et al. (2019) there was high relatedness only neutral mutation and genetic drift, which in turn implies among their Rapa Nui samples. This could indicate that the population is not in equilibrium. This provides that those samples shared parents, and as a result alternative hypotheses to high gene flow, namely that any interpretation of population structure should be the population may have experienced a genetic bottle- treated with caution as these are unreliable. Alvial neck or may be experiencing purifying selection, either et al.’s (2019) samples from Rapa Nui also suggested of which would also explain the apparent lack of any ge- signatures of a recent genetic bottleneck; we do not netic structure in P. flavescens across its range. have samples from Rapa Nui, but when their sam- K-means clustering (a means of partitioning observa- ples were included in analyses with our samples we tions into clusters, which in this case partitions individ- found no evidence of genetic structure in COI. If all uals into putative populations) suggests that there may the samples Alvial et al. (2019) collected during one be several populations, although none of the “popula- trip from Rapa Nui were of closely related individuals, tions” identified can be easily interpreted in terms of a that could explain the perhaps artificially high level biological pattern. For example, the clusters contained of structure between Rapa Nui and the Americas. Fu- individuals from several countries, continents, and ture work should include collection of a time series hemispheres, across different sampling years, with no of Pantala from Rapa Nui to determine whether Alvial apparent pattern. Similarly, we did not find any statisti- et al.’s (2019) results can be confirmed with samples cal support for putative populations defined on the ba- with lower relatedness. Indeed, future work should sis of various geographic scenarios. F values were low- examine more nuclear data to evaluate what biogeo- ST est when grouping individuals as Northern/Southern graphical patterns are reflected in the nuclear genome Hemisphere or Eastern/Western Hemisphere popula- of P. flavescens; there can be genetic discordance be- tions. While there is greater apparent genetic structure tween gene fragments due to phenomena such as in- among some continents, this may be biased by unequal complete lineage sorting and introgression, while pat- sampling across regions. terns of events such as hybridization and isolation can Migrate-N (a program which models migration us- lead to biogeographical discordance (Toews & Brels- ing genetic datasets) suggests that of two tested sce- ford, 2012). In addition to augmenting knowledge of narios, east-west migration is more probable compared Rapa Nui, future efforts should also aim to increase to north-south migration; this can be interpreted as genetic sampling across the continent of Africa. More more common, suggesting the east-west direction is a generally, we encourage further study with more and more common migration route than north-south migra- larger samples, and a range of genetic markers. We tion which makes sense given Pantala flavescens’ use note in particular that for COI (for which at least 71 of intertropical convergence zone winds. F values for haplotypes have been identified) very much larger ST COI recovered between P. flavescens groups were low sample sizes than used here would be desirable, and (North-South hemisphere scenario; East-West hemi- also that the study of faster-evolving genetic markers sphere scenario), which also suggests high gene flow than COI may be needed to fully elucidate the popula- among hemispheres (see also Table 2 for additional tion structure of P. flavescens. population genetics statistics). The F values were simi- Our stable isotope data support the general percep- ST lar to those of the highly migratory Anax junius (micro- tion that P. flavescens is a highly vagile species. The satellite F = 0.02-0.08, Matthews, 2007), and lower majority of individuals investigated here were classed ST than those of somewhat migratory Ischnura hastata as immigrants, which means we collected individuals (Say, 1839) (ND1 F = 0.105, Lorenzo-Carballa et al., from areas other than their natal regions regardless ST 2010), or Cordulegaster sarracenia Abbott & Hibbitts, of where or when we sampled. Stable isotope analy- 2011 (microsatellite F  = 0.423, Abbott et al., 2018). It sis has previously shown that P. flavescens collected in ST is possible that by using a more rapidly evolving marker, Maldives were all immigrants (Hobson et al., 2012a), such as SNPs or microsatellites, additional structuring while those collected in Japan were mostly immigrants between areas might be revealed. Alvial et al. (2019), (Hobson et al., 2021). The fact that we observed only however, showed very small differences in Hd, π, and small numbers of potentially local individuals in our number of haplotypes whether COI or microsatellites samples confirms that most adults disperse soon after were used for P. flavescens, although with microsatel- emergence (e.g. Corbet, 1964; Kumar, 1972). Increased lites some structure was found between what they char- isotope sampling in areas like the Andes, where we had acterised as “continental” and “insular” regions. Their only potentially local individuals, and particularly on continental samples were limited to Central and South isolated islands (coupled with paternity studies) could America. Microsatellites may evolve more rapidly than assist with establishing whether there are exceptions to COI (e.g., McDonald & Potts, 1997), and thus if P. flave- the patterns suggested by our results here, (i.e. wheth- scens colonized the Americas, most likely from Africa, in er there are any resident non-migratory populations the relatively recent past, it is possible that such popula- in other parts of the range of P. flavescens). However, tion structure might be apparent in micros atellite data even in Hawaii, two out of four sampled individuals (2019) but not in COI data (as demonstrated here). were classed as immigrants. International Journal of Odonatology │ Volume 25 │ pp. 43–55 51 Ware, Kohli, Mendoza, ... & Suhling Evidence for widespread gene flow and migration in Pantala flavescens Finally, a literature review confirmed large-scale Physiological and morphological features that may movements of P. flavescens between locations where it allow individuals to efficiently migrate long distances is commonly found (text S1, Table S3; Figure S2). There likely involve their ability to efficiently use lipids as fuel are several possible pathways of population exchange, rather than carbohydrates (e.g., Kallapur & George, often as large swarms, between continents and even 1973; it has been suggested that dragonfly fat reserves across large stretches of open ocean (text S1). Most of can sustain eight hours of flapping flight), as well as dis- these observations were made in association with par- tinct wing morphology (e.g., Alvial et al., 2019; Moore, ticular wind patterns such as the Inter-Tropical Conver- 1993; Sacchi & Hardersen, 2013; Suárez-Tovar & Sar- gence Zone (ITCZ) rain fronts, but also associated with mien to, 2016; Outomuro & Johansson, 2019; Zhao et typhoons and cyclones. Noted dispersal events imply al., 2012), and unique thoracic musculature (Bäumler that individuals are able to cover distances of several et al., 2018). Further, a major precondition (or trait) for thousand kilometers even over oceans and seas, e.g., being an obligate dragonfly migrant may be the abil- Anderson, 2009; Drake et al., 1995; Sparrow et al., ity to develop successfully in ephemeral freshwater 2020). This supports the idea that migration and dis- habitats, which are usually formed after seasonal rains persal may be wind-driven or wind-supported. We as- (Corbet, 1999; Suhling et al., 2019). This seasonality in sume that seasonal appearances at off-shore islands can breeding habitat is directly linked to the propensity of provide clues for reconstructing cross-ocean migration P. flavescens to follow the ITCZ. routes (cf. Anderson, 2009). One such migration route has been described by Anderson (2009) crossing the In- dian Ocean between India and Africa via the Maldives Conclusions and Seychelles Islands, and there is additional evidence for cross-Indian Ocean travels (text S1); Hobson et al. Using COI data from samples taken across most of P. fla- (2021) also describe oceanic crossing by individuals col- vescens’ range, we found high gene flow and a lack of lected on Japanese islands. As many Pacific islands are genetic structure, with one main haplotype shared also visited by P. flavescens, there is support for the globally. This result is supported by isotope data which idea of intercontinental pathways via the Pacific Ocean. suggest frequent long-distance migration. We also com- Whether P. flavescens can regularly cross the eastern Pa- piled extensive observational evidence from the litera- cific barrier remains unclear. The most likely route con- ture which supports these findings. We conclude that necting the Eastern and Western Hemispheres seems there is full genetic exchange among all continents, i.e., to be via the central Atlantic Ocean from West Africa, a globally panmictic population of P. flavescens seems where individuals would be supported by trade winds possible, even if very isolated islands may be excluded blowing towards northern South America and the Carib- from regular migration. The mechanism supporting pan- bean. This dispersal path has been demonstrated for two mixia is regular long-range migration within and between other migratory Odonata of African origin, Anax ephip- continents, even across oceans. Trans-oceanic migration piger (Burmeister, 1839) and Tramea basilaris (Palis ot sustained by trade winds and additional storm events de Beauvois, 1807) (Lambret et al., 2013; Hedlund et explain the ability of individuals to reach distant islands al., 2020; Meurgey & Poiron, 2012) as well as the desert and to travel between continents, while active flight and locust Schistocerca gregaria (Forsskål, 1775) (Lovejoy et movement of seasonal rains in relation to the ITCZ can al., 2006). Seasonal winds may also assist the movement explain patterns observed over continents. All evidence of swarms of P. flavescens as far away from the tropics as points at a mechanism we call the pantropical Pantala about 52°N in East Asia (Borisov, 2012; Cao et al., 2015), conveyor belt, a mechanism which makes it possible for North America (May, 2013) and, recently, also in parts of an insect with excellent flying ability and rapid develop- Europe (Buczyński et al., 2014); in these high latitude re- ment to achieve some level of global panmixia. gions they probably do not survive the winter (Borisov, 2012; Borisov & Malikova, 2019; May, 2013). In Asia and Africa, movements seem to be mainly associated with Acknowledgements the seasonal shifts of the ITCZ and the associated rain- We thank the authors of prior studies on Pantala for making their bearing winds, while reasons for massive overshoot of sequences available on NCBI. Andreas Martens and Oleg Kosterin the ITCZ in North America are not yet fully understood provided relevant literature and made valuable comments. Ware (May, 2013). In the Southern Hemisphere most larger would like to acknowledge funding from NSF DBI #1564386. Ciara landmasses except for Antarctica and southern Patago- Mendoza and Dan Troast began work on this project while students nia regularly receive migratory swarms of Pantala (al- at Rutgers University, Newark. though P. flavescens has been found at least as far south as 37°S in Argentina; del Palacio et al., 2017). Our con- clusion, based on observations compiled from literature Author contributions (supplementary material), is that the extraordinary dis- JLW, GS, DT, KAH, and FS conceived the ideas; JLW, CMM, DT, KAH, persal abilities of Pantala flavescens coupled with sea- GS, FS, and HJ collected the data; JLW, MKK, FS, and KAH analysed sonal wind patterns, support the observed global pat- the data; JLW, GS, and FS led the writing; all authors contributed to tern of gene flow and migration. writing of the manuscript and figure & table preparation. International Journal of Odonatology │ Volume 25 │ pp. 43–55 52