Molecular phylogeny and diversification timing of the Nemouridae family (Insecta, Plecoptera) in the Japanese Archipelago

The generation of the high species diversity of insects in Japan was profoundly influenced by the formation of the Japanese Archipelago. We explored the species diversification and biogeographical history of the Nemouridae Billberg, 1820 family in the Japanese Archipelago using mitochondrial DNA and nuclear DNA markers. We collected 49 species among four genera: Indonemoura Baumann, 1975; Protonemura Kempny, 1898; Amphinemura, Ris 1902 and Nemoura Latreille, 1796 in Japan, China, South Korea and North America. We estimated their divergence times—based on three molecular clock node calibrations—using Bayesian phylogeography approaches. Our results suggested that Japanese Archipelago formation events resulted in diversification events in the middle of the Cretaceous (<120 Ma), speciation in the Paleogene (<50 Ma) and intra-species diversification segregated into eastern and western Japan of the Fossa Magna region at late Neogene (20 Ma). The Indonemoura samples were genetically separated into two clades—that of Mainland China and that of Japan. The Japanese clade clustered with the Nemouridae species from North America, suggesting the possibility of a colonisation event prior to the formation of the Japanese Archipelago. We believe that our results enhanced the understanding both of the origin of the species and of local species distribution in the Japanese Archipelago.


Introduction
The East Asian region-and in particular, the Japanese Archipelago-is considered to have high insect biodiversity [1], [2]. The high degree of Japanese insect biodiversity is a result of several mechanisms-in particular, the complex geological history. The Japanese Archipelago originated in the middle of the Miocene [3] as an independent formation of eastern and western Japanese landmasses. Extensive geographical changes and large-scale climatic changes throughout the islands facilitated the subsequent connection and disconnection of Japanese landmasses from the Eurasian continent, and the formation of tectonic lines (as the median tectonic line, MTL; and the Itoigawa-Shizuoka tectonic line, ISLT) [3], [4], [5]. These geological events-facilitating for the colonisation of insects from the continent and their subsequent PLOS ONE | https://doi.org/10.1371/journal.pone.0210269 January 11, 2019 1 / 17 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Study sites and sample collection
Our sampling sites in Japan comprised 32 sampling sites on Hokkaido Island, 83 on Honshu Island and 27 on Shikoku Island. The field sampling was conducted without the necessity of any special permit. None of the Nemouridae species was found on Hokkaido Island during sampling. All species reported from Hokkaido are known to occur on either Honshu or Shikoku Islands. Herein, we only reported on the sampling sites wherein specimens were found. We collected 20, 7, 8 and 1 species of the genera Nemoura, Amphinemura, Protonemura and Indonemoura, respectively, on 110 sampling sites in Japan (Fig 1, S1 Table). Additionally, 14 species collected from 8 sampling sites of Mainland China and 2 of South Korea (S1 Table, Fig  2) and 100 specimens of the three species Zapada columbiana (Claassen, 1923), Z. cinctipes (Banks, 1897), and Podmosta delicatula (Claassen, 1923) of subfamily Nemourinae collected from 15 sampling sites of North America (western United States of America and Alaska) were included in our analysis. We added these samples from outside of Japan because of their The maps were prepared in QGIS software v 2.18 under the GNU free Documentation License. We collected adult insects using hand nets around riversides. We stored samples in 80% ethanol in the field, and replaced the ethanol with fresh 99.5% ethanol after morphological identification. We identified individuals according to the taxonomical keys of [21], [23], [24], [25], [26], [27], [28], [29], [30], [31] and [32]. Undescribed species in our study were morphologically sorted based on our taxonomic expertise and inconclusive taxonomic keys.

DNA extraction, amplification and sequencing
We genetically analysed a total of 289 individuals, out of which 189 were from East Asia (males, 97; females, 92) and 100 were from North America (males, 92; females, 8). We extracted genomic DNA individually using DNeasy tissue kits (Qiagen GmbH, Hilden, Germany), following the manufacturer's instructions. We amplified a 658-bp fragment of mtDNA cox1 using LCO-1490 and HCO-2198 primers [33] with an annealing temperature of 38˚C and 40 PCR cycles. Further, we amplified a 328-bp fragment of nDNA marker histone 3 (H3) using the universal primers H3F and H3R [34] with an annealing temperature of 58˚C and 40 PCR cycles. We purified the PCR products using the QIAquick PCR Purification Kit (Qiagen GmbH, Hilden, Germany) and sequenced them in both directions using the same primers as mentioned above. Cox1 and H3 sequences were sequenced by Eurofins Operon (Tokyo, Japan). All sequence data reported here have been deposited in GenBank (COI: MK132196-MK132472; H3: MK132473-MK132742).

Sequence analysis
We assembled and edited forward and reverse sequences using CodonCode Aligner v 3.5 (Codon Code Corporation, Dedham, USA). All sequences were aligned using ClustalW (https://www.genome.jp/tools-bin/clustalw) [35]. Alignment of H3 was based on the conservation of the amino acid reading frame across the data. Only a few nucleotide sites in the sequence showed heterozygous nucleotides (i.e. double peaks). These nucleotide sites were treated as unambiguous and they were omitted from our phylogenetic analysis. All alignment data have been deposited in Dryad repository.
We calculated the genetic diversity by the number of polymorphic sites, number of haplotypes and both mean nucleotide substitution rate (i.e. individuals within species) and pairwise nucleotide substitution rate (i.e. between species), with the Kimura 2-parameter model. We performed all analyses using DnaSp v5.10 [36]. All analyses were performed for cox1 and h3 separately.
All sequences of the mtDNA and nDNA markers were compared with the NCBI nucleotide database using blastn queries (http://blast.ncbi.nlm.nih.gov) to corroborate species identification (DNA barcoding, similarity > 98%) and to discard possible sequence errors.

DNA species delimitation
To corroborate the morphological species identification match with our molecular data, we implemented a DNA species delimitation analysis. Putative DNA species were delineated using the General Mixed Yule Coalescent model (GMYC; [37]). An ultrametric gene tree of cox1 gene was constructed using BEST v1.8.3 [38], and the GMYC analysis was performed using the splits package [39] in R ver. 3.3 (R Core Team). We used a single version of the GMYC model. The maximum likelihood of the GMYC model was tested using the likelihood ratio test against a one-species null model (i.e. entire tree is considered as a single coalescent).

Molecular clock analysis
We estimated the evolutionary history of the family in the Japanese Archipelago according to the timing of the divergence of the lineages. For this estimation, we implemented a Bayesian phylogenetic analysis in combination with a molecular clock analysis using BEAST v.2.4.4 [40] with Zwicknia bifrons (Newman, 1838) (Capniidae Banks, 1990) as a outgroup for cox1 and H3 (own sequences) separately. This outgroup was selected owing to their close phylogenetic relationship with Nemouridae [41], [42]. To observe the divergence time, we adopted a relax clock model [43] following a log normal distribution, and calibrated the phylogenetic tree nodes using three types of molecular clock analysis. The first calibration was based on fossil records of the Nemouridae family [44]. We calibrated the nodes at 180 million years ago (Ma) and adjusted the parameters with a standard deviation of 20 Ma, as suggested in a previous study [45] for a 95% highest posterior density (HPD). For this analysis, we implemented a fossilised birth death model [46] for tree prior parameter. The second calibration was based on the Japanese Archipelago formation events dated from 15 to 30 Ma [3]. We applied several calibrations from 15, 20, 25 and 30 Ma at all nodes representing taxonomic species. All calibrations were adjusted to 5 Ma as a standard deviation for a 95% HPD and a fossilised birth death model [46] for tree prior parameter. Lastly, the third calibration was the time to the most recent common ancestor (TMRCA) to observe species diversification patterns based on the mean substitution rate of cox1. Using a Yule model tree prior parameter [47], we applied the substitution rate for insect cox1 of 1.5% [48] and 3.54% [49] per million years for a 95% HPD.
For all branch age calibrations (namely, fossil, biogeographic and mtDNA substitution rate), we performed MCMC for 50 million generations, and log dating trees (BEAST parameters) for every 5000 generations. We tested the output files for convergence after removing a 10% burn-in by examining the effective sampling size using Tracer v1.5 [50]. We pooled the four resulting output trees from biogeographical calibration analysis into a single tree. We then pooled the resulting single tree from biogeographical calibration and the single tree from fossil calibration analyses into a single tree. We performed all pooling analyses using Log Combiner v1.6.1 (BEAST package) summarised with Tree Annotator (BEAST package) and visualised using FigTree v1.3.1 [51]. We performed the analyses for cox1 and H3 separately. In summary, three molecular clock calibration trees were obtained, two for cox1 (fossil + biogeography and cox1 rates) and one for H3 (fossil + biogeography). The incongruence length difference test (ILD) [52] was conducted to test the congruence of tree topologies (fossil + biogeography) between cox1 and H3 using Tree Analysis Technology (TNT) [53]. ILD test revealed no significant differences in terms of the Bayesian tree topologies between cox1 and H3 (P = 0.8); therefore, both markers were polled into a single tree for further analysis. The tree files have been deposited in the Dryad repository.

Phylogenetic analysis between Nemoura from Japan and North America
To observe the phylogenetic relationship between Nemouridae from Japan and North America (Zapada columbiana, Z. cinctipes and Podmosta delicatula), we analysed the maximum likelihood (ML) phylogenetic trees of cox1 and H3 separately using PhyML 3.1 [54]. The General Time-Reversible (GTR) model and gamma distribution were selected for both markers (cox1 and H3) based on a separate test performed with jModel Test v.3 [55] and using Zwicknia bifrons (Capniidae) as an outgroup as described above. The trees were bootstrapped using 10,000 replications.

Genetic diversity and DNA phylogeny
For studying the phylogeny of Nemouridae in the Japanese Archipelago, we analysed two molecular markers. Cox1 sequences were of 658 bp length, with 247 polymorphic sites, 237 parsimony-informative sites, 10 singletons and a mean nucleotide substitution rate of 0.151. H3 sequences were of 328 bp length, with 67 polymorphic sites, 54 parsimony-informative sites, 13 singletons and a mean nucleotide substitution rate of 0.051. No gaps were detected for either cox1 or H3 sequences. In total, for cox1 and H3, we identified 128 and 68 haplotypes, respectively.
The log-likelihood of the GMYC model at the optimal threshold (111.02) was significantly better than the null model of a single coalescent (logL = 56.99) in the likelihood ratio test (p < 0.001). Most clades have GMYC-support values higher than 0.9, implying that the probability of the clades being delimited as separate GMYC-species among the alternative models of delimitation (within a 95% confidence set) is higher than 0.9. The single-threshold model delimited 61 putative species (S1 Table) (Kawai, 1960), N. uenoi Kawai, 1954, N. chinonis Okamoto, 1922, and N. cf. cercispinosa Kawai, 1960) showed two putative DNA-species. While A. decemseta showed multiple putative DNA-species (three putative DNA-species), N. sanbena Shimizu (1993) and P. kohnoae , showed two putative DNA-species in the same sampling site suggesting the presence of cryptic species. The congruence of H3 phylogenetic groups provided confirmation of DNA-based groups detected by GMYC.
We observed the genetic diversity of the species per island (Table 1). Honshu had the highest number of species (26 species), haplotype richness (63) and mean nucleotide substitution rate (average 0.027). Five species were found throughout the three Japanese islands (Honshu, Shikoku and Kyushu), i.e. A. decemseta, A. zonata, N. cf. cercispinosa, N. chinonis and I. nohirae, with a mean nucleotide substitution rate ranging from 0.011 to 0.126 and a total of 23 haplotypes. N. sanbena haplotypes were observed in two different branches in the phylogenetic tree, both within N. cf. cercispinosa and as an isolated branch.

Divergence dates
The Bayesian phylogenetic trees for cox1 and H3 showed tree topology similarity (ILD test, P = 0.8). Three clades corresponded to the three families-Protonemura, Amphinemura and Nemoura-whereas Indonemoura was divided into two clades-the Mainland China clade, clustered with Protonemura, and the Japanese clade (Fig 2).
The evolutionary divergence between the Nemouridae and Capniidae families was settled at 180 Ma, with a 95% HPD interval of 160 to 198 Ma, in the Jurassic geological period (Fig 2). Genus-level diversifications within Nemouridae occurred in the early and middle Cretaceous.

Phylogeographic pattern between Nemoura from Japan and North America
DNA sequences in the Japanese clade of Indonemoura (single species, I. nohirae) showed a high homology with those in the Alaskan species of Z. columbiana (COI: KM874174; >93% sequence similarity) and Z. cinctipes (H3: EF622600; >98% sequence similarity) based on blastn results. Based on the Blast results and the previous reported geographical distribution of these species in Asia and North America, we decided to observe the phylogenetic relationships of these Alaskan specimens (Z. columbiana, Z. cinctipes and P. delicatula) with stoneflies taxa from Japan. The ML phylogenetic trees for both cox1 and H3 (Fig 3) showed that the Indonemoura Japanese clade clustered with three North American species (Z. columbiana, Z. cinctipes and P. delicatula) and the Indonemoura Mainland China clade clustered with the East Asian Nemouridae genera (Nemoura, Protonemura and Amphinemura). The pairwise nucleotide substitution rate based on cox1 between the Indonemoura Japanese clade and Zapada spp. or P. delicatula from North America ranged from 0.13 to 0.15, whereas a higher pairwise nucleotide substitution rate based on cox1 of 0.26 was observed between the Indonemoura Japanese and Mainland China clades (Table 2).

Discussion
We studied mitochondrial cox1 and nuclear H3 gene sequences to determine the patterns of diversification and phylogenetic relationships of species belonging to four genera of stoneflies of the Nemouridae family in the Japanese Archipelago. We estimated the divergence among Nemoura, Amphinemura, Indonemoura and Protonemura to have occurred in the early and mid-Cretaceous (around 100 Ma), which is compatible with previous studies based on fossil records [56], [57]. Our results suggested that these four genera might have dispersed and colonised different areas of the Eurasian continent-including the Japanese landmasses-when they were still connected to the Eurasian continent. Among the four genera, the diversification of Indonemoura occurred earlier (120 Ma) than that of the other genera (<100 Ma), suggesting that it is an ancient genus. The geological isolation of colonised areas [58], the long evolutionary time [59] and poor dispersal ability of Indonemoura [23], [60] might have accounted for their ancient diversification.
Based on the phylogenetic relationships of both molecular markers (cox1 and H3), we observed that the three genera Nemoura, Amphinemura and Protonemura were monophyletic and clustered as three independent groups, as previously observed by morphological systematics [16]. However, Indonemoura was paraphyletic. This genus was divided into two clades corresponding to the Mainland China and the Japanese clades. Surprisingly, the Japanese clade of Indonemoura (single species, I. nohirae) clustered together with North American species (Z. columbiana, Z. cinctipes and P. delicatula), with a low pairwise nucleotide substitution rate (<0.15). The distribution range of these two North America genera covers North America and Eastern Asia. Previous studies suggested that their distribution could be related to the land connection (i.e. the islands) between Alaska and Eastern Asia [15]. Dispersal by island connectivity between Alaska, the Aleutian Islands, the Kamchatka peninsula and the Kuril Islands has been observed in other stonefly families (for instance, Arcynopteryx dichroa (McLachlan, 1872), Capnia nearctica Banks, 1919, Mesocapnia variabilis (Klapálek, 1920) and Nemoura arctica Esben-Petersen, 1910) [61]. However, the distribution of Indonemoura on these islands is unknown.
The complex history of the geological formation of the Japanese Archipelago may provide a possible alternative explanation. The ancestral Japanese landmasses were located on the borders of four major tectonic plates, of which two are continental plates-the Eurasian plate and the North American plate [4] (S2 Fig). The eastern Japanese landmass was located on the North American plate, whereas the western Japanese landmass was located on the Eurasian plate [5]. The dispersal and colonisation of Indonemoura might have occurred from the North American plate to the Eurasian continent or vice versa (from the Eurasian continent to the North American plate) before their geographic separation in an ancient time (around 70 to 80 Ma) [62]. Dispersal events between Eurasian and Japanese landmasses are commonly reported for aquatic insects [10], [63]. Particularly, a dispersal event between North America and the Japanese Archipelago was detected by the phylogenetic relationship of the monophyletic group of caddisflies, Palaeagapetus spp. [9]. However, no prior studies have observed speciation events of aquatic insects associated with geological events that occurred in ancient times (>12 Ma). Our result suggests an ancient divergence time and a distribution pattern of Indonemoura, consistent with a hypothesis of an ancient colonisation influenced by the connection of the Japanese landmass with the North American plate in the Eurasian continent. Nemouridae species diversification, as has been observed in other species of aquatic insects, such as beetles [10], caddisflies [9], water bugs [14] and mayflies [8], [13], was also observed to be affected by the geological formation of the Japanese Archipelago. The diversification of the Nemouridae species occurred during the Paleogene period (<50 Ma). This geological period is consistent with the movement of landmasses (S2 Fig) about 70 Ma ago [4] and the active geological formation of the Japanese Archipelago around 20 Ma ago [5], which could be the cause of the Nemouridae diversification, as previously reported for the mayfly Dipteromimus flavipterus Tojo and Matsukawa, 2003 (35 Ma) [2].
Indonemoura nohirae is the single species of Indonemoura on the Japanese Archipelago [25], [26]. The morphology of their terminalia resembles that of Protonemura rather than of Indonemoura, but the characteristic gill formula justifies their taxonomical classification in Indonemoura [25], [26]. To date, there are 24 Indonemoura species from China [16], [24] and 30 species belonging to the Himalayan and Oriental regions in East Asia [15], [20]. These species are morphologically different from I. nohirae in Japan [15], [16], [20], [24], [25], [26]. We hypothesise that the Indonemoura species of East Asia could be forming separate phylogenetic clades clustered by geographical regions. For the hypothesis testing, further collection of molecular data on Indonemoura from wider areas such as Northeast China, Southeast China, Mongolia, Russia and other countries in Asia is needed in future studies.
Eight species (I. nohirae, A. decemseta, A. zonata, A. longispina, A. megaloba, N. chinonis, N. uenoi and N. cercispinosa) showed interesting patterns of intra-species separation into two genetic groups corresponding to eastern and western areas of the Fossa Magna region of Honshu Island (S1 Fig). Honshu is the centre of insect biodiversity [10]; apart from its extensive territorial space, it is the main island with a geological history [3], [4], [5]. We found supporting evidence on the genetic diversity of these eight species. We found a larger mean nucleotide substitution rate and haplotype number in the Honshu region than in other islands ( Table 1). The mean nucleotide substitution rate and haplotype diversity are indications of biodiversity [64], which could lead to evidence of speciation [65]. Out of eight species, the diversification of six species (A. decemseta, A. zonata, A. longispina, A. megaloba, N. uenoi and N [5], [66]. The speciation of aquatic insects was often observed to be influenced by these two geological events [2]. Additionally, species diversification-from eastern or western Japan of the Fossa Magna regionshowed recent diversification events (3 to 5 Ma) corresponding with the formation of the small islands in northeastern or southwestern edge areas of Japan (Fig 1). The northeastern islands created land bridges between the Japanese Archipelago and China or Korea, whereas the southwestern islands connected Taiwan or the Philippines with the Japanese Archipelago [5], [66]. This connectivity promoted immigration events in Japan that might have contributed to the formation of the current genetic diversity, as previously observed in mayflies [13] and beetles [10].
The evolutionary divergence of the Nemouridae family was promoted by the complex geological formation of the Japanese Archipelago. Despite the different evolutionary rate of both molecular markers, Bayesian analysis found congruence between them; however, failed to find congruence with their morphological taxonomy. The main morphological character used for identification of adult stoneflies species is its genital morphology. The evolution of genital morphology is; however, governed by within-population sexual selection rather than environmental or geological history of the locations [67]. Conversely, the genetic variation of natural populations has been observed to be directly associated with environmental [68] and geological variations [2]. Therefore, the genetic variation could reflect an independent course in the evolutionary history of Indonemoura than do the morphological characters used for their taxonomy. However, we detected that N. sanbena shared haplotypes from different lineages, revealing a possible introgression or incomplete sorting of ancestral polymorphisms [10]. This is an often reported phenomenon in stoneflies [40], [69], [70], which remains as unresolved species. Resolving the problems between the process of evolution of morphological characters and the genetic variation within species will improve our future understanding of the origin of the species and the local species distribution.
Finally, our inference of divergence time was based on the coalescent simulation approach. Despite the frequent use of this approach, a biased sampling of lineages and extreme statedependent molecular substitutions rate heterogeneity are known to potentially cause erroneous inference of divergence time [71]. One of the causes of biased sampling is under-sampling (i.e. the incomplete-coverage samples included in the phylogenetic tree). Previous studies demonstrated that under-sampling increases the estimates of the posterior probabilities (i.e. variance of age estimates become less precise) [72] and led to bias and low precision of the divergence time of shallow nodes (i.e. evolutionary recently divergent taxa) [73], [74]. However, under-sampling tends not to affect deep nodes (i.e. internal nodes close to the root) time scales estimations or the tree shape [75], [76], [77]. Therefore, in order to confirm the accuracy of the time scales estimations of the recently divergent taxa in our study, additional taxonomic samples would be required to be included in the molecular clock analysis. Conversely, the state-dependent molecular substitution rate heterogeneity could be address by combining node calibrations generated by more than one calibration analyses, as recommended by [71], [78]. A cautious method such as the combined uses of fossil records and biogeographic ages as employed in our analysis may minimize the risk of such erroneous inference.  I. nohirae, A. decemseta, A. zonata, A. longispina,  A. megaloba, N. chinonis, N. uenoi and N. cf. cercispinosa (GMYC = 2 species). (TIFF)

S2 Fig. Putative formation of the Japanese Archipelago [2], [3], [4]
, [5]. (A) Around 30 to 130 Ma, the Japanese landmasses were located in two major tectonic plates from the Eurasian continent. (B) Around 15 to 30 Ma, the Japanese landmasses began to separate from Eurasia and the North American Plates began to separate from the Eurasian continent, and remained separated by a sea zone called Fossa Magna-a geological event called double-door. (C) Current map of the Japanese Archipelago in East Asia, where the names of the four main Japanese islands and the two tectonic lines are shown. The maps was prepared using QGIS v 2.18 under the GNU free Documentation License with political boundaries from the Global Database of Administrative Areas (https://gadm.org/). (TIFF)