Taxonomic revision of Aegista subchinensis (Möllendorff, 1884) (Stylommatophora, Bradybaenidae) and a description of a new species of Aegista from eastern Taiwan based on multilocus phylogeny and comparative morphology

Abstract Aegista subchinensis (Möllendorff, 1884) is a widely distributed land snail species with morphological variation and endemic to Taiwan. Three genetic markers (partial sequence of the mitochondrial cytochrome c oxidase subunit I [COI], the 16S rDNA and the nuclear internal transcribed spacer 2 [ITS2]) were analysed to infer phylogenetic relationships and genetic divergence of closely related species of the genus Aegista, Aegista vermis (Reeve, 1852) and Aegista oculus (Pfeiffer, 1850). A new species from Aegista subchinensis has been recognized on the basis of phylogenetic and morphological evidences. The nominal new species, Aegista diversifamilia sp. n. is distinguished from Aegista subchinensis (Möllendorff, 1884) by its larger shell size, aperture and apex angle; wider umbilicus and flatter shell shape. The northernmost distribution of Aegista diversifamilia sp. n. is limited by the Lanyang River, which is presumed to mark the geographic barrier between Aegista diversifamilia sp. n. and Aegista subchinensis.


Introduction
Traditional morphology-based taxonomy provides a window to explore biodiversity and evolutionary history. The application of molecular genetic markers opens new avenues to discover biodiversity. In recent years, it was found that the species richness of land snails bearing comparably few morphological characteristics and exhibiting limited abilities of dispersal had been underestimated once molecular tools were applied (Hirano et al. 2014, Nantarat et al. 2014, Nishi and Sota 2007, Prevot et al. 2013, Wu et al. 2008. In contrast, taxonomic revision based merely on a single molecular locus may lead to an overestimation of number of taxa. Integrative taxonomy was therefore proposed to integrate multiple independent lines of evidence for objective taxonomic treatment (Dayrat 2005, Padial et al. 2010, Schlick-Steiner et al. 2010. Taiwan is a continental island that was formed through the collision of the Philippine Sea plate and the Eurasian plate. This collision uplifted the Central Mountain Range (CMR), forming a major physical barrier for animals inhabit lowland areas. The CMR has contributed to evolutionary divergences between organisms on either side of the CMR both on interspecific and intraspecific levels (Huang and Lin 2010, Jang-Liaw et al. 2008, Kuo et al. 2014, Lin et al. 2012, Wang et al. 2007, Yu et al. 2014. Aegista subchinensis (Möllendorff, 1884) is one of the most widely distributed species of the genus Aegista in Taiwan and is commonly found in lowland hardwood forests near the CMR (Hsieh 2003, Hsieh et al. 2006, Hsieh et al. 2013, Lee and Chen 2003, Lee and Wu 2004. The morphological differences observed between western and eastern populations indicate that the evolutionary diversification of this species complex may be underestimated and requiring further investigation Chen 2003, Lee andWu 2004). Based on multilocus sequence analyses and comparative morphology, we demonstrate that Aegista snails from eastern Taiwan, originally identified as A. subchinensis, represent a new species which is herein described as A. diversifamilia sp. n.

Sample collection and molecular techniques
Live snails identified as A. subchinensis were collected from ten localities in Taiwan. Similar species, A. vermis (Reeve, 1852) of Ishigaki Island and A. oculus (Pfeiffer, 1850) of Miyako Island, were collected from two and four localities, respectively, on the southern Ryukyu Islands. Four congenic species, A. mackensii (Adams & Reeve, 1850), A. granti (Pfeiffer, 1865), A. inrinensis (Pilsbry & Hirase, 1905), and A. shermani (Pfeiffer, 1865), distributed in Taiwan were used as outgroups to root the phylogenetic tree. Global positioning system (GPS) coordinates of sampling sites (including latitude, longitude and altitude) were recorded using Garmin GPSmap 60CSx with an uncertainty of less than 10 metres ( Figure 2 and Table 1). Vouchers and type specimens of A. subchinensis and A. diversifamilia sp. n. were deposited in the National Museum of Natural Science, Taiwan (NMNS, NMNS-7276) and the Natural History Museum, United Kingdom (NHMUK 20140070-20140074). Snails were relaxed in water for at least 6 hours, quickly fixed by submersion in boiling water and then preserved in 95% ethanol. DNA was extracted from 10 mg of foot tissue using AxyPrep™ Multisource Genomic DNA Miniprep Kit (Axygen Bioscience, USA) following the manufacturer's protocol. A partial sequence of mitochondrial cytochrome c oxidase subunit I (COI) was amplified using the LCO1490 and HCO2198 universal primers (Folmer et al. 1994). Partial 16S ribosomal DNA was amplified using the 16Sar and 16Sbr universal primers (Palumbi et al. 1991). Complete nuclear internal transcribed spacer 2 (ITS2) was amplified using the LSU1 and LSU3 primers (Wade et al. 2006). The PCR mixture was composed of 10 ng DNA template, 1 μM primers, 1X Taq DNA polymerase 2.0 Master mix kit (Ampliqon, Denmark) and water. The total volume of the PCR mixture was 20 μl. PCR was performed under the following conditions: initial denaturation at 94 °C for 1 min followed by 36 cycles of denaturing at 94 °C for 30 s, annealing at 48 °C or 52 °C for 30 s and a final extension at 72 °C for 30 s. Primer annealing temperature was 48 °C for COI and 52 °C for 16S and ITS2. The size of the PCR products was checked under ultraviolet light after gel electrophoresis. The PCR mixture was purified using Genomics Universal DNA Purification kit (Genomics BioSci and Tech, Taiwan). Sanger sequencing was performed on an ABI PRISM 3730 DNA Analyzer at Institute of Cellular and Organismic Biology, Academia Sinica.

Phylogenetic reconstruction
Sequences were visually checked using BIOEDIT version 7.2.5 (Hall 1999) and deposited in GenBank (KJ574281-KJ574400, Table 1). The sequence of ingroups made available by Hirano et al. (2014) were incorporated into the phylogenetic reconstruction (Table 1)  Maximum likelihood phylogeny and sampling sites of Aegista spp. Reconstructed phylogeny was based on concatenated sequences of mitochondrial COI, 16S and nuclear ITS2 genes. Branch support confidences of clades are shown in bootstrap, approximate likelihood-ratio test and Bayesian posterior probability, respectively. The log likelihood of maximum likelihood tree = -6584.1713. The numbered sampling sites are detailed in Table 1.  (Darriba et al. 2012, Guindon andGascuel 2003). The model filtering threshold of heuristic was set to 1.0 to increase the efficiency for evaluating 56 substitution models. The best-fit model was evaluated using BIC and DT criteria because their accuracy has been shown to outperform hLRT and AIC (Luo et al. 2010). Phylogenetic relationship was inferred from each gene and concatenated sequences using maximum likelihood and the Bayesian inference method. Phylogeny of maximum likelihood was inferred using PHYML 3.0  implemented in SEAVIEW version 4.4.2 (Gouy et al. 2010). Parameters were set to empirical nucleotide frequencies with optimized invariable sites; the number of rate categories equalled four; tree searching operation as the best of nearest neighbour interchange (NNI) and subtree pruning and regrafting (SPR); BioNJ with optimized tree topology as the starting tree. Branch support confidences were estimated via traditional bootstrap with 100 replicates (Felsenstein 1985) and an fast and accurate alternative method, approximate likelihood-ratio test (aLRT) (Anisimova and Gascuel 2006). Phylogeny of Bayesian inference was performed using MRBAYES version 3.2 (Ronquist et al. 2012). Two parallel runs of three heated chains and one cold chain were performed for 2-5 million generations and sampled every 1000 generations. Sampling was stopped when all chains reached stationarity and the runs converged (split frequencies < 0.01). The first 25% of sampling trees were discarded as burn-in and a 50% majority-rule consensus tree was computed. Branch support confidences were inferred from Bayesian posterior probability. Network is an alternative way to infer the phylogenetic relationship among haplotypes (Bandelt et al. 1999, Bapteste et al. 2013, Huson and Bryant 2006. Haplotype network was reconstructed by medianjoining method (Bandelt et al. 1999) via NETWORK version 4.612 (available at http://www.fluxus-engineering.com/sharenet.htm). Interspecific and intraspecific genetic distances were calculated by using PAUP* version 4.0 (Swofford 2003) employing the best-fit substitution model but excluding sequences from Hirano et al. (2014).

Morphological analyses
For genital morphological comparison, we dissected two samples of one adult snail each from the western (Zhishanyan, Taipei City) and the eastern (Heren 1, Xiulin Township, Hualien County) population, respectively (Table 1). Snails were relaxed in water for more than 6 hours, fixed by submersion in boiling water for 20 seconds, and preserved in 70% ethanol. Shell was crushed before dissection. All soft tissues were preserved in 70% ethanol. Terminology of genital morphology follows Gómez (2001). Images of genitalia were obtained via Canon PowerShot A650IS digital camera. QPcard 201 was included in each image for colour and scale correction with PhotoImpact X3 (Corel Corp., Taipei, Taiwan). The following organs were measured in cm and scaled by their shell width: length of albumen gland (AG), hermaphroditic duct (HD), spermoviduct (SOD), free oviduct (FO), vagina (V), dart-sac of auxiliary copulatory organ (DS), penis (P), epiphallus (E), and epiphallus flagellum (EF). The following ratios were calculated: HD/AG, AG/SOD, HD/SOD, FO/V, EF/E, E/P, and P/V. Shell morphological comparisons were based on examination of 36 specimens of A. subchinensis from the west of the CMR and 43 specimens from the east (Table 1). Morphological differences between the western and the eastern A. subchinensis were analysed on the basis of measurements and shape outline coordinates. Shell images were obtained as mentioned above. The image's background was removed using CLIPPING MAGIC (available at https://clippingmagic.com/). Images were converted into thinplate splines (TPS) format using TPSUTIL version 1.56 (Rohlf 2013). The Janssen's method of counting the number of whorls was used (Janssen 2007). Fifteen characteristics were measured: the number of whorls; width of the shell (SW), aperture (AW), umbilicus (UW), first whorl (FW), and 2 nd -6 th whorl (2W-6W); height of the shell (SH), aperture (AH), body whorl (BH), and secondary body whorl (SBH); and angle of apex (AA) (Figure 1). Measurements of length (in cm) and angle were conducted using TPSDIG version 1.4 (Rohlf 2004) and IMAGEJ version 1.47 (Abramoff et al. 2004), respectively. The shell flatness and aperture shape were expressed by calculating ratios of shell height divided by shell width (SH/SW), aperture height divided by aperture width (AH/AW) and umbilicus width divided by shell width (UW/SW), and additional ratios were calculated: AW/SW, BH/SW, SBH/SW, FW/SW, AH/SH, SBH/SH, SBH/ SH, SBH/BH, SH/UW, AW/UW, AH/UW, BH/UW, SBH/UW, 2W/3W, 3W/4W, 4W/5W, 5W/6W. Summary statistics and 95% confidence intervals of these characteristics were calculated for the western and the eastern population of A. subchinensis. Statistical differences of each characteristic were analysed using Welch's t test or the Mann-Whitney U test when the distribution of characteristics did or did not follow a normal distribution, respectively. Multivariate analysis of variance (MANOVA) was conducted with the species as independent variables and the characteristics mentioned above as dependent variables. A Bonferroni correction was applied for p-values. All characteristics were log-transformed for principle component analysis (PCA) and discriminant analysis with Jackknifed classification via PAST version 3.01 (Hammer et al. 2001). Outline coordinates of each shell were digitalized using TPSDIG for geometric morphometric analysis of shell shape. Outline coordinates were approximated with 16 control points via MORPHOMATICA version 1.6 (Linhart et al. 2006). The outline coordinates of the western and the eastern A. subchinensis were used to calculate mean shell shape in MORPHOMATICA. A relative warps PCA was conducted using PAST.

Molecular phylogeny
A total of 50 individuals were sequenced from the Aegista ingroup. The lengths of the COI, 16S and ITS2 after alignment were 655 bp, 280 bp and 750 bp, respectively. The best maximum likelihood tree and Bayesian consensus phylogram had similar topologies: Aegista oculus from Miyako Island formed a monophyly outside the other three species (Figure 2). Aegista vermis and A. caerulea from Ishigaki Island formed a sister pair with high branch support except for one individual of A. vermis from Tozato, Ishigaki Island (sampling site 32), which was found in an unresolved position in Bayesian phylogeny. The cluster of A. subchinensis revealed a basal bifurcation between eastern and western populations. The western and eastern A. subchinensis formed reciprocal monophyletic clades with highly supported by the bootstrap, aLRT and Bayesian posterior probability. However, the monophyly of the A. subchinensis clade was only weakly supported. Low statistical branch support for the monophyly of the A. subchinensis clade probably resulted from conflicting phylogenetic relationships between the COI and the 16S gene tree and unresolved phylogenetic relationships among species of the ITS2 gene tree (Suppl. material 1, Figures  S1-6). Maximum likelihood and Bayesian COI gene trees suggested that the eastern A. subchinensis was the sister clade of A. vermis and the western A. subchinensis (Suppl. mate- rial 1, Figures S1 and S4). Maximum likelihood 16S gene tree showed sister relationship between the western and the eastern A. subchinensis with low branch support (Suppl. material 1, Figure S2), but the eastern A. subchinensis was sister species to A. vermis with moderate branch support in Bayesian tree (Suppl. material 1, Figure S5). COI gene tree estimated by maximum likelihood and Bayesian inference and the 16S Bayesian gene tree suggested that the western and the eastern A. subchinensis were not sister clades.
No haplotypes were shared among species for COI, 16S or ITS2 genes (Figure 3). For COI haplotype network, the eastern A. subchinensis was nested with the other three species for at least 16 mutations, no clear sister species relationship could be inferred. Haplotype network of 16S gene suggested the western and the eastern A. subchinensis were distant related. The western A. subchinensis was nested with A. oculus at least 9 mutations. The eastern A. subchinensis was nested with A. vermis at least 8 mutations. For ITS2 haplotype network, the eastern A. subchinensis was more closely related to A. vermis that separated by at least 3 mutations. The eastern and western A. subchinensis were diverged at least 5 mutations.
Gene trees and haplotype networks suggested the western and the eastern A. subchinensis were not sister clades. The eastern A. subchinensis was more closely related to A. vermis that distributed in Ishigaki Island, Japan. The absence of shared haplotype between the eastern and the western A. subchinensis might suggest that they were diverged and currently no gene flow between them.
The mean genetic distance between the western and the eastern A. subchinensis clades was 5.9% for COI, 4.2% for 16S, and 0.8% for ITS2. This divergence corresponded to the divergence between other closely related congeneric species (Table 2). Based on phylogenetic trees, haplotype networks and genetic distance analyses, the western and the eastern A. subchinensis were diverged and might not be sister clades. More genetic markers are needed to resolve the low statistical support of phylogenetic relationships among species.

Morphological analyses
The major differences of genital morphology of the western and the eastern A. subchinensis were the length of AG, DS, EF, the shape of DS, and the number of lobes of mucus gland in the auxiliary copulatory organ (M) (Figure 4, Table 3). The length of AG of the eastern A. subchinensis (0.88, scaled by shell width) was three times longer than the western (0.24). The length of DS and EF were nearly two times longer in the eastern A. subchinensis (DS=0.45, EF=0.57) than the western (DS=0.23, EF=0.29). The shape of DS was more rounded and larger in the eastern A. subchinensis than the western. The eastern A. subchinensis had three lobes of M and the western A. subchinensis had two (Figure 4). The HD/AG ratio was the most different characteristics between the eastern and the western A. subchinensis that showed HD was three times longer than AG in the western A. subchinensis but HD was shorter than AG in the eastern (Table 3). The eastern A. subchinensis showed larger values than the western in AG/SOD, EF/E and E/P ratio.   The eastern and western populations of A. subchinensis differed significantly from each other in all studied shell parameters (p < 0.001) except the number of whorls and the height of the secondary body whorl (Table 4, Suppl. material 1, Table S1). The eastern A. subchinensis had a similar number of whorls (p > 0.05) but shells were significantly larger than those of the western A. subchinensis and had a wider apex angle. The eastern and western populations also differed significantly in all morphometric ratios except UW/SW, AW/UW, 2W/3W, 3W/4W and 5W/6W. MANOVA suggested that the morphological difference between the western and the eastern A. subchinensis was statistical significant (Bonferroni-corrected p-value=1.12E-10). Results of the PCA of morphological characteristics suggested that principle component axis 1 (PC1) explained 46.32% of the total variation and had the highest loading scores for aperture height (0.284). PC2 explained 20.03% of the total variation and had the highest loading scores for height of the secondary body whorl (0.476). The PCA scores plot of PC1 and PC2 showed that morphology between the eastern and the western A. subchinensis were well distinguished with very limited overlap ( Figure 5A). Discriminant analysis showed that specimens could be correctly classified (100%) into eastern and western clades (93.67% using Jackknifed analysis). Relative warps PCA of shell outline coordinates suggested that PC1 and PC2 represented 94.43% and 3.07% of the total variation, respectively. The scores plot of relative warps PC1 and PC2 also showed prominent morphological differences between the eastern and the western clade ( Figure 5B). The mean shape of the western and the eastern A. subchinensis were presented in Figure 6. The eastern A. subchinensis ( Figure 6B) was more flat than the western A. subchinensis ( Figure 6A). The type locality of A. subchinensis was Tamsui in northern Taiwan (Figure 2). The illustration of A. subchinensis presented in plate 7 figure 8 of the original description (Möllendorff 1884) showed higher conic shape that was similar to our analysed western populations of A. subchinensis. Despite we did not have sample from Tamsui, the type locality of A. subchinensis, we sampled much wider geographical region around Tamsui that included more variation of molecular and morphological characteristics. Based on the observed amounts of morphological and genetic differentiation, we conclude that the eastern and western populations assigned to A. subchinensis have diverged into separate species. Phylogeny reconstructed from concatenated sequences supports monophyly of both clades corresponding to their allopatric distributional pattern that separated by the Lanyang River. The Lanyang River was a biogeographic barrier for a high elevation mammal, Formosan wood mouse Apodemos semotus (Hsu et al. 2001). Some researchers identified the Xueshan Mountain Range, located in the northern area of Lanyang River, as a biogeographic barrier for lowland animals (Lin et al. 2012, Shih et al. 2011, Yu et al. 2014. To our knowledge, this is the first study revealed that the Lanyang River as a barrier for lowland terrestrial animals. We suggested that the eastern A. subchinensis might be diverged from the western A. subchinensis by vicariance event.  Table 4, Suppl. material 1, Table S1 Aegista subchinensis Hsieh, 2003: 200, figs;Lee and Chen 2003: 234, figs above text, figs 1-2; Lee and Wu 2004: 13-14, figures 2A, 3D;Hsieh et al. 2006: 250, figs;Wu and Jian 2006: fig 33;Hsieh et al. 2013: 335, figs. Aegista (Aegista) subchinensis Hemmen and Niederhöfer, 2007:   Other material examined. Anpingkeng, Dongshan Township, I-Lan County, 24°36'52.5"N, 121°46'38.1"E (3 adult dry shells); Wushibi, Su'ao Township, 24°29'13.5"N, 121°50'02.9"E (1 juvenile in EtOH); Chaoyang Trail, Nan'ao Township, 24°27'35.9"N, 121°48'53.9"E (2 adult dry shells); Heren 1, Xiulin Township, Hualien County, 24°14'49.1"N, 121°43'06.4"E (1 adult dry shells); Heren 2, 24°14'54.8"N, 121°42'51.4"E (7 adult dry shells); Heren Trail, 24°13'58.5"N, 121°42'27.73"E (1 adult and 4 juvenile in EtOH); Southern Chongde Tunnel, 24°11'22.0"N, 121°39'36.8"E (2 juvenile in EtOH); Sanjianwu, 24°10'55.3"N, 121°37'34.3"E (6 adult dry shells); Taroko Service Center, 24°09'31.9"N, 121°37'20.7"E (6 adult dry shells); Badagang, 24°10'36.8"N, 121°33'43.6"E (1 adult and 4 juvenile in EtOH, 6 adult dry shells) (materials mentioned above were deposited in NMNS, NMNH-7276); Hoping Forest Road, (1 adult in EtOH, NMNS-004875-00015 and 1 adult dry shell, NMNS-004962-00038); Sanzhan Northern Stream, (1 adult dry shell, NMNS-003348-00023).

Systematics
Description. Shell Morphology. Shell depressed globosed, dextral, medium sized, shell width range 1.98-3.24 cm, shell height range 0.97-1.68 cm, shell height/ shell width ratio range 0.43-0.55. Shell thin but solid, glossy with chestnut brown or yellowish-brown, usually with narrow and light brown spiral band on periphery. Shell surface with distinct oblique and curved growth lines. Apex obtuse, angle range 148.56°-165.02°. Spire depressed conic, slightly convex, suture depressed. Whorl range 6.6-8.2 in number, earlier whorl narrow then slowly increases regularly, and last whorl shouldered. Body whorl height range 0.53-0.88 cm. Aperture little descending, ovate or nearly circular, width range 0.78-1.32 cm, height range 0.48-1.05 cm. Peristome white, expanded and reflected. Umbilicus widely open, width range 0.77-1.59 cm. Mean and standard errors of each characteristics were provided in Table 4. Morphological measurements of all specimens were presented in Suppl. material 1, Table S1.
Genital morphology. Atrium thick and short. Penis slender and long. Epiphallus slender, longer than penis. Penis retractor muscle thin and long, attached to one-third part of epiphallus. Epiphallial flagellum thick and long, logner than epiphallus, wider than penis and epiphallus. Dart-sac of auxiliary copulatory organ thick and large, inserted into the base of vagina, with one small accessory-sac of auxiliary copulatory organ. Three mucus glands of the auxiliary copulatory organ. Vagina slender at the base of dart-sac, gradual wider and thick toward free oviduct, inflated at the connected region of free oviduct, about equal length of penis. Free oviduct thick, short,inflated. Pedunculus of bursa copulatory thin and long. Sac of bursa copulatory large and oval. Vas deferens thin and long, wider than penis retractor muscle. Spermoviduct long, about four times longer than penis and oviduct. Hermaphroditic duct slender and long, about half length of spermoviduct. Albumen gland thick and long, longer than hermaphroditic duct.
Etymology. Named after the recent efforts supporting equal marriage rights in Taiwan and around the world. Derived from "diversus" (Latin for different) and "familia" (Latin for family), adjective of feminine gender.
Distribution. Endemic to Taiwan and is currently known from I-Lan and Hualian Counties. Aegista diversifamilia sp. n. is absent from Gueishan Island based on our Table 4. Measurements (in cm) of Aegista subchinensis and A. diversifamilia sp. n. Mean, standard error, statistical method and the p-value were provided.  Table 1).

Ecology.
Live snails are generally found on the ground or under leaf litter in shady, moist environments in lowland hardwood forests (Figure 8). Eggs white and round, approximately 3 mm in diameter with 20-30 eggs in each spawn (personal observation of reared snail in laboratory).
Remarks. Aegista diversifamilia sp. n. can be distinguished from A. subchinensis by its overall larger shell width (1.98-3.24 cm), whorl width and aperture, more depressed shell, and wider umbilicus (0.77-1.59 cm) and larger apex angle (148.56°-165.02°) (see Suppl. material 1, Table S1). For the genital morphology, A. diversifamilia sp. n. was distinguished from A. subchinensis by thicker and about three times logner albumen gland, larger and about two times longer dart-sac of auxiliary copulatory organ and epiphallial flagellum. The length of hermaphroditic duct/ albumen gland ratio was three times larger in A. diversifamilia sp. n. than in A. subchinensis.
The morphological divergence between the eastern and the western A. subchinensis was firstly noticed by Lee and Chen (2003), who found that the shells from the western population were roughly one third smaller than those from the eastern population. When newly describing A. caperata, Lee and Wu (2004) suggested the presence of cryptic species within A. subchinensis from different sides of the CMR. Wen and Hwang (2014) compared reproductive system between subgenus Aegista and Plectotropis. Aegista (Aegista) subchinensis and A. (Plectotropis) mackensii were dissected as representative species. Wen and Hwang (2014) mentioned there were two lobes of mucus glands of A. subchinensis. According to the shell masurements, sampling locality (Xiulin Township, Hualien County) and the illustration of genital morphology figure 1 of A. subchinensis was actually the nominal new species presented here, A. diversifamilia sp. n. It might suggested that the number of lobes of mucus glands is a variable characteristic in A. diversifamilia sp. n. .