MAIN
A dramatic decline in the biomass and biodiversity of species and populations of insects has been reported in recent years, and climate change is thought to be one of the contributing factors (Hallmann et al., 2017; Lister and Garcia, 2018; Sánchez-Bayo and Wyckhuys, 2019). This decline is expected to have severe impacts on global biodiversity due to the important role that insects play in various ecosystems (Thomas et al., 2004). In particular, for species adapted to cool climates, it is expected that their effective population size will expand and/or shrink, and their distribution will shift to higher latitudes and altitudes. These factors increase the risk of extinction of these cool climate species (Menéndez et al., 2006; Colwell et al., 2008). To evaluate the impact of climate change on species in subarctic and cool temperate forests and ecosystems, it is important to examine population genetic diversity, structure and demographic history of the species that inhabit them. However, although population genetic studies have been conducted for cool temperate forest tree species across Japan (e.g., Tsuda et al., 2015), insect species are yet to be well examined and there is less information about the population demographic history of insects distributed in cool temperate forests.
Yezoterpnosia nigricosta (Motschulsky, 1866) is a cicada distributed in subarctic and cool temperate forests in Japan, China and the Russian Far East. In the northern area of Japan, Y. nigricosta is widespread and occurs continuously from low altitudes to subalpine forests, while in the central and southern parts of Japan, the species’ distribution is restricted to the higher altitudes of cool temperate forests, which are mainly dominated by beech (Fagus crenata) and oak (Quercus crispula) trees. Thus, this species is listed as a threatened or near-threatened species on the Red List of 12 prefectures in Japan (Association of Wildlife Research and EnVision, 2007, http://jpnrdb.com/search.php?mode=map&q=07150110701). The species is characterized by its emergence period from spring to early summer (e.g., May–July), much earlier than other cicada species such as Cryptotympana facialis (Walker, 1858) and Graptopsaltria nigrofuscata (Motschulsky, 1866) (Hayashi and Saisho, 2011). According to a recent taxonomic study by Lee (2012), this species was transferred from the genus Terpnosia Distant, 1892 to Yezoterpnosia Matsumura, 1917. In the present study we focused on Y. nigricosta as an environmental indicator species in Japan, representative of cool temperate forest insects. A population genetic study was therefore performed to elucidate the genetic structure and population demographic history of this species and to evaluate the impacts of climate change.
Recent next-generation sequencing (NGS)-based technology has enabled the efficient development of not only single-nucleotide polymorphisms but also microsatellite markers (Drechsler et al., 2013; Yamakawa et al., 2019). Thus, 17 polymorphic microsatellites were isolated from Y. nigricosta through restriction site-associated DNA sequencing (RAD-seq) and their genetic variation was evaluated using two populations in Japan.
An adult individual was collected in Sapporo, Hokkaido for genome sequencing. In addition, eight individuals were collected from Hokkaido to Kyushu, covering the species’ range in Japan, for screening polymorphic loci. To evaluate the genetic diversity of the developed loci, 17 individuals were collected from the Sugadaira Research Station, Mountain Science Center, University of Tsukuba (36°31′ N/138°20′ E) in Nagano Prefecture; and 15 individuals were collected from the Hokkaido Research Center, Forestry and Forest Products Research Institute (42°59′ N/141°23′ E) in Hokkaido. The whole body or a foreleg of the samples was preserved with 99.5% ethanol and stored at –20 ℃ until DNA extraction.
Total genomic DNA for NGS library preparation was extracted from flight muscle of the individual collected in Sapporo, Hokkaido, using a commercial DNA extraction kit, NucleoSpin Tissue (Macherey-Nagel, Düren, Germany). Construction of the RAD-seq library and sequencing with a high-throughput DNA sequencing platform was conducted by Novogene (Beijing, China). Briefly, the genomic DNA was digested by the restriction enzyme EcoRI. After ligating adapters and random shearing, 300–500-bp fragments with both adapters were selected and amplified by PCR. Paired-end sequencing with 150 bp length at each end was performed on HiSeq X Ten (Illumina, San Diego, CA, USA). Raw sequence data were assembled with a cloud computing service of Read Annotation Pipeline (Nagasaki et al., 2013), provided by the DNA Data Bank of Japan (DDBJ). After preprocessing by trimming low-quality bases, de novo assembly was conducted with Velvet Version 1.2.10 (Zerbino and Birney, 2008). From assembled contigs, sequences including microsatellite(s) were extracted and primers for each locus were designed using the software QDD version 3.1 (Meglécz et al., 2014). From the list outputted, sequences which included one microsatellite composed of at least six repetitions of a 2–6-bp motif, but not any homopolymers or nanosatellites, were chosen. Because some of the contigs could be misassembled by Velvet, the accuracy of target sequences was confirmed by mapping raw sequence reads to the sequences extracted from the contigs with Geneious R10 (Biomatters, Auckland, New Zealand).
After preliminary PCR amplification tests for all 51 designed primer pairs, the level of polymorphism of amplified loci was checked using eight individuals collected from Hokkaido to Kyushu, covering the species’ range in Japan. PCR amplification was carried out using the Type-it Microsatellite PCR Kit (Qiagen, Valencia, CA, USA), based on two methods of fluorescent labeling of primers, depending on loci (Table 1). First, in the method using fluorescently labeled universal primers (Blacket et al., 2012), each 5-μl reaction contained 2.5 μl of 2× Type-it Multiplex PCR Master Mix, 0.5 μl of 0.01 μM forward tailed primer and 0.2 μM each of fluorescently labeled universal primer and reverse primer, 1.2 μl of RNase-free water and 2–10 ng/μl of genomic DNA. Second, the traditional approach in which forward primers have fluorescent dyes added was employed, with each 5-μl reaction containing 2.5 μl of 2× Type-it Multiplex PCR Master Mix, 0.5 μl of 0.2 μM each of forward and reverse primer, 1.0 μl of RNase-free water and 2–10 ng/μl of genomic DNA. In both methods, the amplification process consisted of an initial denaturation at 95 ℃ for 5 min; 30 cycles of denaturation at 95 ℃ for 30 s, annealing at 57 ℃ for 90 s and extension at 72 ℃ for 30 s; and a final extension at 60 ℃ for 30 min. Fragment sizes were determined using an ABI PRISM 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) and GeneMarker software (SoftGenetics, State College, PA, USA) with GeneScan 500 LIZ dye Size Standard v2.0 (Applied Biosystems).
Table 1. Characteristics of 17 microsatellite primers developed for
Yezoterpnosia nigricosta| Locus | Primer sequences (5′-3′) | Repeat motif | Ta (℃) | Fluorescent label* | Directly labeling locus-specific | Size range (bp) | DDBJ accession no. |
|---|
| TeNi002 | F: TGCTGAAGAGGTATCCCACT | (AT)6 | 57 | | PET | 241–255 | LC537600 |
| R: CCATTCCCATTGCCTGCAAG |
| TeNi003 | F: ACGTAGCAGAAAGACTAGTGGA | (GC)6 | 57 | Tail A | | 101–103 | LC537601 |
| R: ACTAGTGTGATAACCGCCACT |
| TeNi012 | F: TTCTCCGCCTGAATATAACACA | (AC)6 | 57 | Tail A | | 185–187 | LC537602 |
| R: CGCATTGCAAAGGAATGTGC |
| TeNi019 | F: TGATCGAATTCGGAACAGGCA | (GAT)6 | 57 | | FAM | 103–127 | LC537603 |
| R: TCCAAAGGATCACAGTGGCA |
| TeNi028 | F: CCATGGTTATGTTGTTTATGCGC | (TG)6 | 57 | Tail D | | 152–158 | LC537604 |
| R: ACATCGTATCCAGACATGCAGT |
| TeNi029 | F: TGGAAGGTAGTGGCTTGAGTG | (GC)6 | 57 | Tail D | | 158–160 | LC537605 |
| R: TGGGCACAATTTCTTTCACGC |
| TeNi032 | F: TATCATCTTGGCAGCGGGTG | (CA)6 | 57 | | HEX | 130–142 | LC537606 |
| R: AAGCAACATCTGTGCATGCG |
| TeNi038 | F: TCACACGTGATGGAAACCGA | (CA)7 | 57 | Tail A | | 139–153 | LC537607 |
| R: CGCAGTAAATACACGAGTGCG |
| TeNi040 | F: GCGTACGAGTAAGTTATGCGC | (GT)7 | 57 | | NED | 109–122 | LC537608 |
| R: AACAGAAATCGGACTGATCATTAAA |
| TeNi041 | F: TAAGTACAGGGTCGGCCGAT | (GT)7 | 57 | Tail C | | 134–140 | LC537609 |
| R: AAGGCGCGTTGTACAGAGAT |
| TeNi045 | F: CGTTCACACTCTTAGTTTCCCA | (TG)8 | 57 | Tail D | | 109–131 | LC537610 |
| R: AAGGTAACATATGGAACAGTGGA |
| TeNi047 | F: GGCGTGTACGTACGATTTGC | (AG)8 | 57 | | HEX | 228–288 | LC537611 |
| R: AAGCGAGAATTCAGACGGGA |
| TeNi048 | F: AAATACAGCCTTGCGTGGAA | (TG)8 | 57 | | NED | 160–168 | LC537612 |
| R: TACTGACACTCGCTTCAGGC |
| TeNi050 | F: GATGTATTGACGTGCGCGAG | (GT)9 | 57 | | HEX | 80–86 | LC537613 |
| R: CACTTGTATGTTTACCGTCAAACG |
| TeNi051 | F: CGAGATGGCTAGACAGAAAGGT | (AG)9 | 57 | | FAM | 92–122 | LC537614 |
| R: CTGTGGCGGCTTAGATTTGC |
| TeNi052 | F: GTAGAGGTGCTCAGCTGGTG | (GA)9 | 57 | | FAM | 163–169 | LC537615 |
| R: ACAGCCAGCTCTTTACTGCA |
| TeNi053 | F: TTATCGAAGCGAGGAAAAGC | (AC)10 | 57 | Tail B | | 165–195 | LC537616 |
| R: CGGTCTGGAACGCTTTAAGT |
Ta, annealing temperature; * sequence of the fluorescent labels: Tail A = 5′-GCCTCCCTCGCGCCA-3′, Tail B = 5′-GCCTTGCCAGCCCGC-3′, Tail C = 5′-CGGAGAGCCGAGAGGTG-3′, Tail D = 5′-CAGGACCAGGCTACCGTG-3′.
The genetic diversity of the loci that showed amplification and polymorphism was evaluated using 32 individuals from two populations of Y. nigricosta as mentioned above. Total genomic DNA of these samples was extracted using the DNeasy Blood & Tissue kit (Qiagen), and PCR and genotyping were conducted by the same methods as used in the polymorphism check mentioned above.
For each population, the number of alleles (A), expected heterozygosity (HE) and fixation index (FIS) at each locus were evaluated using FSTAT version 2.9.3 (Goudet, 1995). The deviation of FIS value from zero was tested by 1,000 randomizations at each locus and population using FSTAT to check Hardy–Weinberg equilibrium. In addition, genetic differentiation between the two populations was evaluated by calculating FST (Weir and Cockerham, 1984) and its standardized value, F’ST, which ranges from 0 to 1 (Meirmans and Hedrick, 2011), using GenAlEx version 6.5 (Peakall and Smouse, 2012).
A total of 22,233,102 read pairs (approximately 6.7 Gb) obtained from RAD-seq were assembled into 2,401,270 contigs (maximum contig size = 7,699 bp, minimum contig size = 45 bp, N50 contig size = 218 bp). From the contigs, 793 sequences containing microsatellites were extracted by the QDD pipeline. Next, 135 sequences with a pure microsatellite composed of six or more repetitions of a motif were selected. After excluding misassembled sequences, 51 loci remained. Finally, 17 polymorphic loci with clear fragment patterns that could be easily genotyped were obtained (Table 1). The number of alleles across the two populations was 1–15 (mean = 4.294) and three loci were not polymorphic in the Sapporo population (Table 2). HE and FIS per locus were 0.000–0.914 (mean = 0.519) and −0.225–0.456 (mean = 0.108), respectively (Table 2). Although FIS values showed negative and positive values, these values did not significantly deviate from 0, suggesting they were in Hardy–Weinberg equilibrium. FST and F’ST between the two populations were 0.111 and 0.237, respectively, suggesting moderate genetic differentiation, probably due to a distribution shift and demographic history of cool temperate forests in Japan, in relation to the Quaternary ice ages (Tsuda et al., 2015).
Table 2. Genetic variation of the 17 microsatellite loci for two populations of
Yezoterpnosia nigricosta| Locus | Sugadaira population (N = 17) | Sapporo population (N = 15) |
|---|
| A | HE | FIS | A | HE | FIS |
|---|
| TeNi002 | 7 | 0.673 | 0.301 | 3 | 0.657 | 0.188 |
| TeNi003 | 2 | 0.500 | 0.125 | 2 | 0.514 | −0.167 |
| TeNi012 | 2 | 0.257 | −0.143 | 2 | 0.238 | −0.12 |
| TeNi019 | 5 | 0.665 | −0.061 | 5 | 0.664 | −0.004 |
| TeNi028 | 4 | 0.285 | −0.095 | 1 | 0.000 | NA |
| TeNi029 | 2 | 0.165 | −0.067 | 1 | 0.000 | NA |
| TeNi032 | 7 | 0.743 | 0.287 | 4 | 0.574 | 0.303 |
| TeNi038 | 5 | 0.590 | 0.152 | 3 | 0.355 | 0.436 |
| TeNi040 | 6 | 0.629 | 0.251 | 5 | 0.710 | 0.154 |
| TeNi041 | 4 | 0.324 | 0.273 | 1 | 0.000 | NA |
| TeNi045 | 9 | 0.840 | 0.090 | 3 | 0.381 | −0.225 |
| TeNi047 | 15 | 0.914 | 0.292 | 8 | 0.843 | 0.051 |
| TeNi048 | 5 | 0.798 | 0.263 | 5 | 0.693 | −0.058 |
| TeNi050 | 3 | 0.540 | 0.456 | 3 | 0.481 | 0.03 |
| TeNi051 | 5 | 0.526 | 0.441 | 3 | 0.605 | 0.008 |
| TeNi052 | 3 | 0.507 | −0.043 | 3 | 0.555 | 0.279 |
| TeNi053 | 4 | 0.627 | −0.196 | 6 | 0.795 | 0.162 |
N, number of analyzed individuals; A, number of alleles; HE, expected heterozygosity; FIS, fixation index; NA, not available.
Overall, 17 polymorphic microsatellite markers were developed for Y. nigricosta, which will be useful for population genetics studies of this species. Additional samples covering the species’ range will be collected and examined in the future, and the genetic structure and past demographic history of the species will be inferred using the loci developed here, in order to evaluate the impacts of climate change and understand the life history strategy and evolution of cicadas.
