History of the terrestrial isopod genus Ligidium in Japan based on phylogeographic analysis

Background
Phylogeographical approaches explain the genetic diversity of local organisms in the context of their geological and geographic environments. Thus, genetic diversity can be a proxy for geological history. Here we propose a genus of woodland isopod, Ligidium, as a marker of geological history in relation to orogeny and the Quaternary glacial cycle.

Results
Mitochondrial analysis of 721 individuals from 97 sites across Japan revealed phylogenetic divergence between the northeastern and southwestern Japan arcs from 7 to 3.5 million years ago. It also showed repeated population expansions in northeastern Japan in response to Quaternary glacial and interglacial cycles. Genome-wide analysis of 83 selected individuals revealed multiple genetic nuclear clusters. The genomic groupings were consistent with the local geographic distribution, indicating that the Ligidium phylogeny reflects its migration history.

Conclusion
Ligidium DNA sequence analysis can provide insight into the geological, geographical, and paleoenvironmental history of the studied region.

Methods
Sample collection
We surveyed Ligidium populations in Japan. We collected 828 Ligidium specimens from 97 sites and sequenced 721 samples (Fig. 1, Table 1). Samples were preserved in 70–99.5% ethanol in 2-mL microtubes at 4°C or room temperature.

DNA extraction, amplification, and sequencing
Genomic DNA was isolated from the muscles of the abdomen and legs with a DNA Mini Kit (Qiagen, Germantown, MD, USA). PCR amplification was conducted using the primer pair LCO-1490 and HCO-2198. Amplification and cycling conditions were as in a previous work. Details of the experimental conditions are provided in the Supplemental Data. Sequencing was conducted with an ABI 3130 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA). Sequences were checked and assembled using MEGA7. Mitochondrial gene locus (CO1) sequences for 721 individuals were determined and aligned using ClustalW.

RAD-seq
We performed RAD-seq analysis to search for SNPs in individuals from Niigata and Hokkaido obtained in previous studies, two regions in northern Japan (Aomori and Sendai), and two regions in western Japan (Shizuoka and Sendai). Tables 1 and S1 list the samples used, and the sampling sites are shown in Fig. 1c. Genomic DNA was isolated from almost whole-body tissues with a DNA Mini Kit (Qiagen). Libraries for RAD-seq were prepared with EcoRI and BglII restriction enzymes. The library was sequenced with 150 + 150 bp paired-end reads in one lane of an Illumina HiSeqX instrument (Illumina, San Diego, CA, USA) by Macrogen (Seoul, South Korea). Raw reads were trimmed using Trimmomatic-0.39 with the following parameters: ILLUMINACLIP: adapter.fasta:2:30:10:keepBothReads, SLIDINGWINDOW: 4:15, CROP: 132, HEADCROP: 2, and MINLEN: 130. Sequences are available at the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (DRA014204). We used two pipeline programs to call the SNPs: denovo_map.pl provided by Stacks and ipyrad. Following a previous work, we varied the combinations of the denovo_map.pl parameters as follows: (n, M) = (2, 1), (3, 2), (4, 3), (5, 6), and selected (n, M) = (2, 1), which called the most SNPs. We used Stacks’ populations program to analyze populations of individual samples, calculate population genetics statistics, and export data in various output formats for analysis. PLINK v1.90b6.18 was used for data handling. Alleles with a frequency of < 1% and sites with > 50% heterozygosity were removed. Only SNPs shared by ~80% of the individuals were retained. With ipyrad, loci with frequencies of > 50% heterozygosity were removed, and SNPs shared by ~70% of the local populations were retained. We retained SNPs shared by at least two individuals and filtered out individuals that did not have 80% of all SNPs using TASSEL 5. After filtering with TASSEL 5, we used PGDSpider to convert the VCF files for other analyses.

Genetic structure analysis
We tested the ability of Structure v. 2.3.4 to determine the genetic structure of populations using Bayesian cluster analysis. Ten simulations were run, with the burn-in period and Markov chain Monte Carlo iterations set to 105 and 106, respectively. The maximum value of K was determined based on the mtDNA results and geographical distribution. For the Structure analysis, one SNP was randomly sampled from each locus to avoid the effect of linkage disequilibrium. The python script vcf_single_snp.py (radseq/vcf_single_snp.py at master · pimbongaerts/radseq · GitHub) was used to obtain the one SNP datum from ipyrad, and drawings were created using the R package pophelper. In addition, PCA was performed to visualize the genetic differences among populations using the adegenet package in R. We obtained pairwise Fst values for the RAD-seq dataset using Arlequin 3.5.1.2. Fst values were used to test population structure, supported by cluster analysis, and statistical significance was based on 1000 restored extractions.