Comparative Morphological, Ultrastructural, and Molecular Studies of Four Cicadinae Species Using Exuvial Legs

Previous studies have suggested that exuviae can be used for the identification of cicada species, but the precise characteristics that differ among species have not been determined. Thus, we performed the first comparative analyses of the leg morphology, ultrastructure, and mitochondrial DNA sequences of exuviae of four dominant cicada species in Korea, Hyalessa maculaticollis (Motschulsky, 1866), Meimuna opalifera (Walker, 1850), Platypleura kaempferi (Fabricius, 1794) and Cryptotympana atrata (Fabricius, 1775), the source of Cicadidae Periostracum, a well-known traditional medicine. A morphological analysis revealed that the profemur length, femoral tooth angle, and distance between the intermediate and last tooth of the femoral comb are useful characteristics for identification. We also evaluated the usefulness of the size, degree of reflex, and number of spines on the mid-legs and hind legs as diagnostic features. An ultrastructural study showed that Meimuna opalifera has a unique surface pattern on the legs. The sequences obtained using exuviae were identical to previously obtained sequences for adult tissues. Moreover, in a phylogenetic analysis using CO1 sequences, each species formed a monophyletic cluster with high bootstrap support. Accordingly, multiple methodological approaches using exuviae might provide highly reliable identification tools. The integrative data provide useful characteristics for the exuviae-based identification of closely related species and for further taxonomic and systematic studies of Cicadinae.


Introduction
Cicadas (Hemiptera, Cicadidae) are plant-sucking insects characterized by loud, species-specific calling songs produced by adult males during the summer [1]. Their immature life stage, which is much longer than the adult stage, continues for one to 17 years underground [2][3][4][5][6]. Cicadas must periodically shed their exoskeleton for growth [5,7]. Exuviae are the remains of exoskeletons after ecdysis (molting) of final instar nymphs and are associated with adult structures [8].
In traditional Chinese medicine, Cicadidae Periostracum has been used for clinical purposes owing to its anticonvulsive, sedative, and hypothermic effects [15]. Modern pharmacological studies of Cicadidae Periostracum have shown that it has fibrinolytic [16] and anti-bacterial activity [17] as

Materials
Exuviae were collected from branch or trunk of trees (i.e., Chionanthus retusus Lindl. and Paxton (Retusa fringetree), Prunus serrulata f. spontanea (E.H. Wilson) C.S. Chang (Oriental flowering cherry), Zelkova serrata (Thunb.) Makino (Sawleaf zelkova), and other kinds of trees in Korean natural populations after molting from July to August 2018. The morphological characteristics of exuvial fore-legs, mid-legs and hind legs of these species were determined.

Morphological Analysis
Morphological characteristics were observed using a stereomicroscope (Olympus SZX16; Tokyo, Japan). Morphological images were captured using a digital camera attached to the microscope (Olympus DP21, Olympus, Tokyo, Japan).
Twenty individuals of each of the four species were used to obtain measurements, for a total of 80 measurements for the fore-, mid-, and hind legs. The following parameters were evaluated: profemur length (FL), protibia length (FTL), femoral tooth angle (FA), number of femoral comb teeth (NFC), distance between the intermediate tooth of the femur and last tooth of the femoral comb (DIF), pretarsus length (FSL), mesotibia length (MTL), mesotarsus length (MSL), number of mesotibial spines (NMS), metatibia length (HTL), metatarsus length (HSL), and number of metatibial spines (NHS).
Most of these measurements were obtained using simple measurement tools (Olympus DP21 digital photography system, Olympus, Tokyo, Japan) based on images obtained using a digital camera (Olympus DP21) for microscopes. A set of digital Vernier calipers (CD-15CP; Mitutoyo, Kawasaki, Japan) was also used for measurements.
The terminology for structural features followed Midgley et al. [26] and Hou et al. [27].

PCA
A principal component analysis (PCA) was conducted to determine whether the quantitative morphological data allowed the grouping of species. The first three principal components with eigenvalues above 1.0 were obtained. These analyses were implemented in PC-ORD version 5.31 [38]. The PCA included twelve variables (FL, FTL, FA, NFC, DIF, FSL, MTL, MSL, NMS, HTL, HSL and NHS).

Ultrastructural Analysis
The fore-legs, mid-legs, and hind legs were excised from the exuviae and washed in 70% ethanol, and all samples were allowed to air-dry. Tarsi and tibia from all three legs and femurs from fore-legs were mounted directly on aluminum stubs using double-sided adhesive carbon tape. Stubs were coated with gold using a sputter coater (208HR; Cressington Scientific Instruments Ltd., Watford, UK), and all samples were observed using a low-voltage field emission scanning electron microscope (JSM-7600F; JEOL, Tokyo, Japan) at an accelerating voltage of 5 kV with a working distance of 8-10 mm. To confirm the consistency of morphological characteristics, at least three samples from each species were evaluated.

Preparation of Genomic DNA
Samples were identified based on morphological characteristics and genomic DNA was extracted from most samples (Table S1). Genomic DNA was isolated according to the manufacturer's protocol using the DNeasy ® Blood and Tissue Kit (Qiagen, Valencia, CA, USA). The purity and concentration of DNA were assessed using a spectrophotometer (Nanodrop ND-1000; Wilmington, DE, USA) and 1.5% agarose gel electrophoresis [39]. The final DNA concentration used for PCR amplification was approximately 15 ng/µL in TE buffer. Extracted DNA samples were stored at −20 • C.
PCR amplification was performed using a DNA Engine Dyad ® PTC-0220 (Bio-Rad, Foster City, CA, USA). The following modified parameters of Kim et al. [39] were used: 95 • C for 5 min; 35 cycles of 1 min at 95 • C, 1 min at 45 • C and 1 min at 72 • C; and a final extension for 5 min at 72 • C. PCR products were separated by 1.5% agarose gel electrophoresis with a 100-bp DNA ladder (Solgent).

Nucleotide Sequence and Phylogenetic Analysis
Amplified mitochondrial CO1 DNA fragments were extracted from agarose gels using a Gel Extraction Kit (Qiagen) and sub-cloned into the pGEM-T Easy Vector (Promega, Madison, WI, USA). Inserted fragments were sequenced in both directions using an automatic DNA sequence analyzer (ABI 3730; Applied Biosystems Inc., Foster City, CA, USA) as described by Kim et al. [39]. The samples of C. atrata (14 individuals: abbreviation CA), H. maculaticollis (17 individuals: HM), M. opalifera (eight individuals: MO), and P. kaempferi (nine individuals: PK) were used in the analysis (Appendix). The following CO1 sequences for the four species deposited in NCBI GenBank were also analyzed: C. atrata, MG737717; H. maculaticollis, KY860344; M. opalifera, GQ527088; and P. kaempferi, MG737816.
Approximately 700-bp CO1 sequences were edited using BioEdit version 7.2.5 [41]. The contigs were aligned to analyze intra-and interspecific sequence variation. For the analysis of sequence identity and divergence, inter-or intraspecific genetic distances were calculated using the Kimura-2-parameter (K2P) model in MEGA 6. A phylogenetic analysis based on the full-length CO1 sequences was performed using MEGA 6 version 6.06 [42,43]. The phylogenetic tree was constructed using the NJ method with the K2P model, pairwise deletion for gaps/missing data treatment and 1000 bootstrap replicates, setting Anthocharis cardamines (Linnaeus, 1758) (MH420365) as an outgroup.

Results
We characterized morphological variation in the three leg types (fore-legs, mid-legs, and hind legs) in the four species. The major leg characteristics are summarized in Table 1  Inserted fragments were sequenced in both directions using an automatic DNA sequence analyzer (ABI 3730; Applied Biosystems Inc., Foster City, CA, USA) as described by Kim et al. [39]. The samples of C. atrata (14 individuals: abbreviation CA), H. maculaticollis (17 individuals: HM), M. opalifera (eight individuals: MO), and P. kaempferi (nine individuals: PK) were used in the analysis (Appendix). The following CO1 sequences for the four species deposited in NCBI GenBank were also analyzed: C. atrata, MG737717; H. maculaticollis, KY860344; M. opalifera, GQ527088; and P. kaempferi, MG737816.
Approximately 700-bp CO1 sequences were edited using BioEdit version 7.2.5 [41]. The contigs were aligned to analyze intra-and interspecific sequence variation. For the analysis of sequence identity and divergence, inter-or intraspecific genetic distances were calculated using the Kimura-2parameter (K2P) model in MEGA 6. A phylogenetic analysis based on the full-length CO1 sequences was performed using MEGA 6 version 6.06 [42,43]. The phylogenetic tree was constructed using the NJ method with the K2P model, pairwise deletion for gaps/missing data treatment and 1000 bootstrap replicates, setting Anthocharis cardamines (Linnaeus, 1758) (MH420365) as an outgroup.

Results
We characterized morphological variation in the three leg types (fore-legs, mid-legs, and hind legs) in the four species. The major leg characteristics are summarized in Table 1
Measurements (Table 1)  Platypleura kaempferi (Fabricius, 1794). Generally yellowish brown. Femur with posterior tooth long and sharp, curved forward very slightly, length approximately three-times longer than the width of its base; accessory tooth small and blunt; intermediate tooth present (Figure 1d); distance between itf and ltf about as long as the DIF ( Table 1); surface of profemur smooth ( Figure 3n); femoral comb usually with seven teeth (all of the individuals with seven teeth), the first tooth about two-or three-times wider than the second tooth (femoral formula 2-1-7) ( Table 1), surface smooth ( Figure  3m). Protibia slightly arched, laterally flattened; apical tooth very short, blunt; point of the blade of tibia very short, small, tooth-like, separated from the apical tooth of blade by an incision (Figure 3o). Pretarsus well-developed, folded over the inner surface of the protibia, two-segmented; apical tarsomere elongated; pretarsal claws of unequal size (Figure 3p). Mesotibia with four, very small and blunt apical spines, strongly reflexed (Figure 2d; 4m), surface smooth with sparse setae (Figure 4n). Mesotarsus smooth with sparse setae (Figure 4o), apex two-segmented; apical tarsomere elongated; mesotarsal claws of strongly unequal size (Figure 4p). Metatibia with four, very small and blunt apical spines, strongly reflexed (Figure 2h; 5m), surface of metatibia smooth with sparse setae ( Figure  5n). Metatarsus smooth with sparse setae (Figure 5o), apex two-segmented; apical tarsomere elongated; metatarsal claws of strongly unequal size (Figure 5p).
Measurements (Table 1)  Platypleura kaempferi (Fabricius, 1794). Generally yellowish brown. Femur with posterior tooth long and sharp, curved forward very slightly, length approximately three-times longer than the width of its base; accessory tooth small and blunt; intermediate tooth present (Figure 1d); distance between itf and ltf about as long as the DIF ( Table 1); surface of profemur smooth ( Figure 3n); femoral comb usually with seven teeth (all of the individuals with seven teeth), the first tooth about two-or three-times wider than the second tooth (femoral formula 2-1-7) ( Table 1), surface smooth (Figure 3m). Protibia slightly arched, laterally flattened; apical tooth very short, blunt; point of the blade of tibia very short, small, tooth-like, separated from the apical tooth of blade by an incision (Figure 3o). Pretarsus well-developed, folded over the inner surface of the protibia, two-segmented; apical tarsomere elongated; pretarsal claws of unequal size (Figure 3p). Mesotibia with four, very small and blunt apical spines, strongly reflexed (Figure 2d; Figure 4m), surface smooth with sparse setae (Figure 4n). Mesotarsus smooth with sparse setae (Figure 4o), apex two-segmented; apical tarsomere elongated; mesotarsal claws of strongly unequal size (Figure 4p). Metatibia with four, very small and blunt apical spines, strongly reflexed (Figure 2h; Figure 5m), surface of metatibia smooth with sparse setae (Figure 5n). Metatarsus smooth with sparse setae (Figure 5o), apex two-segmented; apical tarsomere elongated; metatarsal claws of strongly unequal size (Figure 5p).

Principal Component Analysis
We explored relationships among species based on quantitative leg morphological data using PCA ( Figure 6). The first two principal components (PCs) explained 84.89% of the total variance ( Table 2). PC1 explained 61.74% of the variance in leg size (FL, FTL, FSL, MTL, MSL, HTL, and HSL). PC2 accounted for 23.14% of the data variability, of which indicators of the size of the femoral tooth (FA and DIF) were the significant variables for the grouping of species (Table 2). As shown in a PCA biplot, the OTUs for C. atrata and H. maculaticollis were grouped on the negative side of the PC 1 axis, while OTUs for M. opalifera and P. kaempferi were grouped on the positive side of the PC1 axis ( Figure  6). Most OTUs for M. opalifera were positioned on the positive side of the PC2 axis using NFC. However, OTUs for P. kaempferi were positioned on the central to negative side of the PC2 axis using FA ( Figure 6).  (except b, c, f, j, k, n, o, p, scale bars = 20 µm).

Principal Component Analysis
We explored relationships among species based on quantitative leg morphological data using PCA ( Figure 6). The first two principal components (PCs) explained 84.89% of the total variance ( Table 2). PC1 explained 61.74% of the variance in leg size (FL, FTL, FSL, MTL, MSL, HTL, and HSL). PC2 accounted for 23.14% of the data variability, of which indicators of the size of the femoral tooth (FA and DIF) were the significant variables for the grouping of species (Table 2). As shown in a PCA biplot, the OTUs for C. atrata and H. maculaticollis were grouped on the negative side of the PC 1 axis, while OTUs for M. opalifera and P. kaempferi were grouped on the positive side of the PC1 axis ( Figure 6). Most OTUs for M. opalifera were positioned on the positive side of the PC2 axis using NFC. However, OTUs for P. kaempferi were positioned on the central to negative side of the PC2 axis using FA ( Figure 6).

Sequence Analysis and Phylogenetic Relationships
A 690 bp CO1 region was successfully amplified and sequenced from 48 cicada samples (see Appendix for accession numbers; Table 3). Intraspecific genetic distances were 0.11% for C. atrata, 0.16% for H. maculaticollis, 0.03% for M. opalifera, and 0.21% for P. kaempferi (Table 3). Inter-specific genetic distances ranged from 17.70% to 20.16%, and was highest for P. kaempferi (20.16%) among the three species (Table 3). Additionally, we inferred the phylogenetic relationships among the four Korean cicada species using CO1 sequences ( Figure 7). The phylogenetic tree constructed using the neighbor-joining (NJ) method revealed that all samples from each species formed a monophyletic cluster with high bootstrap support (100%; Figure 7). Samples of H. maculaticollis and M. opalifera were more closely related to each other than to other species. Overall, these data suggest that the four cicada species are identifiable on the basis of sequence variation in the CO1 region.

Sequence Analysis and Phylogenetic Relationships
A 690 bp CO1 region was successfully amplified and sequenced from 48 cicada samples (see Appendix for accession numbers; Table 3). Intraspecific genetic distances were 0.11% for C. atrata, 0.16% for H. maculaticollis, 0.03% for M. opalifera, and 0.21% for P. kaempferi (Table 3). Inter-specific genetic distances ranged from 17.70% to 20.16%, and was highest for P. kaempferi (20.16%) among the three species (Table 3). Additionally, we inferred the phylogenetic relationships among the four Korean cicada species using CO1 sequences ( Figure 7). The phylogenetic tree constructed using the neighbor-joining (NJ) method revealed that all samples from each species formed a monophyletic cluster with high bootstrap support (100%; Figure 7). Samples of H. maculaticollis and M. opalifera were more closely related to each other than to other species. Overall, these data suggest that the four cicada species are identifiable on the basis of sequence variation in the CO1 region.

Discussion
We used multiple methodological approaches in combination to obtain reliable identification information for four cicada species, including morphological and ultrastructural features of the fore-leg (femur, femoral comb, protibial, and pretarsus), mid-leg (mesotibia and mesotarsus), and hind leg (metatibia and metatarsus), as well as CO1 sequences of exuviae. Moreover, we developed an accurate taxonomic key and discussed variation in morphological characteristics in relation to molecular sequence data.
Most morphological studies of cicada exuviae have focused on size-related characteristics, mainly body length, head width, and wing length [23,24,27,28]. We found that qualitative characteristics, including ultrastructural features, as well as quantitative leg traits were valuable diagnostic characteristics for cicadas. Previously, Maccagnan and Martinelli [25] described the three leg pairs in detail and performed a comparative analysis of fifth-instar of cicadas in Brazil. Moreover, Midgley et al. [26] suggested that all three leg pairs should be included for detailed comparisons of cicadas. Thus, descriptions of all leg parts, including qualitative and quantitative characteristics, are essential for the accurate identification of the cicadas.
The fore-legs are expected to be closely related to the burrowing depth in the soil of nymphs and may provide promising characteristics for identification [2,3,[28][29][30]. We were able to infer that fore-legs characteristics, such as the protibia and profemur surface patterns, might indeed be closely related to burrowing depth. P. kaempferi, which live 10 to 30 cm underground during the nymph stage [28,44], have a blunt protibia and smooth profemur, whereas M. opalifera have a sharp protibia and echinate profemur. Hou et al. [45] reported that M. mongolica, which is closely related to M. opalifera, could extend to 60 cm underground. M. mongolica also have very sharp and long protibia [28]. Thus, a sharp protibia and echinate profemur may be effective for deep digging. To verify this, observations of burrowing depth in M. opalifera nymphs and profemur surface pattern of M. mongolica are required. Moreover, ultrastructural patterns of the profemur were a stable characteristic, with consistent patterns at the species level. To evaluate the utility of profemur surface patterns for inferring phylogenetic relationships, further studies with an expanded and targeted sampling strategy based on recent phylogenetic analyses are necessary [46]. A morphological analysis revealed that FA and DIF were useful fore-leg characteristics for distinguishing among this studied species. FA and DIF could be effective diagnostic characteristics for the identification of cicadas based on the exuviae.
The number of the metatibial spines is a phylogenetically important characteristic accordingly to Hou et al. [28], who also reported that C. atrata has five spines on its hind legs. In contrast to the results of Hou et al. [28], we found that C. atrata has six to seven metatibial spines, and the variation in spine numbers can be explained by geographical or individual differences. According to this variation within species, phylogenetic classification based on the number of the metatibial spines needs to be reconsidered. The echinate surface pattern on the mid-and hind legs, only seen in M. opalifera, is useful for distinguishing among species. Thus, this unique pattern might be of taxonomic value, but further ultrastructural studies of other cicada species and comparative analyses of the exuviae and adult traits are needed to better understand their occurrence and taxonomic implications in cicadas. Furthermore, only P. kaempferi had four very small and strongly reflexed apical spines on the tibiae. The size, degree of reflex, and number of spines might be important diagnostic characteristics.
COI sequences have been proposed as an effective barcoding tool for animal identification [47][48][49] and have been successfully used for the identification of cicadas [46,50]. We found that COI sequences enabled accurate classification at the species level. Moreover, sequences obtained from the exuviae were identical to those obtained from adults in previous investigations [46,50]. These results clearly supported the efficiency and usefulness of genetic materials obtained from cicada exuviae, which are easy to collect and could be preserved in a range of environmental conditions [8].
Based on morphological and ultrastructural analyses of exuvial legs, we generated a taxonomic key for four major species in Korea. The observed morphological patterns, which were conserved at the species level, corresponded with phylogenetic relationships. These useful morphological and ultrastructural characteristics, such as profemur length and surface, femoral tooth angle, number of apical spines, and surface of meso-and metatibia of exuviae, may be informative for taxonomic and phylogenetic analyses. The extensive divergence of morphological and ultrastructural characteristics of the legs provides a basis for a taxonomic key for distinguishing among species. Moreover, the leg morphological and exuvial molecular characteristics, which are highly consistent at the species level, might be highly reliable identification markers.

Conclusions
In conclusion, we used statistical analyses and scanning electron microscopy to explore useful characteristics and ultrastructural features of cicada legs. Although broader sampling of species is necessary, these results provide a valuable framework for future studies of evolutionary trends and variation in cicada leg morphology. Moreover, our study demonstrates the importance of integrative studies for obtaining reliable information for species descriptions and for understanding the diversity of cicadas.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/10/7/199/s1, Table S1: Exuvium information of specimen, collection locality and data, with accession number for each sample in this molecular study. Asterisked accession number are sequences obtained in this study.