- Research Article
- Open Access
Population structure of Apodemus flavicollis and comparison to Apodemus sylvaticus in northern Poland based on RAD-seq
BMC Genomics volume 21, Article number: 241 (2020)
Mice of the genus Apodemus are one the most common mammals in the Palaearctic region. Despite their broad range and long history of ecological observations, there are no whole-genome data available for Apodemus, hindering our ability to further exploit the genus in evolutionary and ecological genomics context.
Here we present results from the double-digest restriction site-associated DNA sequencing (ddRAD-seq) on 72 individuals of A. flavicollis and 10 A. sylvaticus from four populations, sampled across 500 km distance in northern Poland. Our data present clear genetic divergence of the two species, with average p-distance, based on 21377 common loci, of 1.51% and a mutation rate of 0.0011 - 0.0019 substitutions per site per million years. We provide a catalogue of 117 highly divergent loci that enable genetic differentiation of the two species in Poland and to a large degree of 20 unrelated samples from several European countries and Tunisia. We also show evidence of admixture between the three A. flavicollis populations but demonstrate that they have negligible average population structure, with largest pairwise FST<0.086.
Our study demonstrates the feasibility of ddRAD-seq in Apodemus and provides the first insights into the population genomics of the species.
Mice of the genus Apodemus (Kaup, 1829) (Rodentia: Muridae) are one the most common mammals in the Palaearctic region . The genus comprises of three subgenera (Sylvaemus, Apodemus and Karstomys) , however the systematic classification of the 20 species belonging to the genus  is not fully settled . In the Western Palearctic, the yellow-necked mice A. flavicollis (Melchior, 1934) and the woodmice A. sylvaticus (Linnaeus, 1758) are widespread, sympatric and occasionally syntopic species. They are often difficult to distinguish morphologically in their southern range , but in the Central and Northern Europe both are easily recognisable by the full yellow collar around the neck of A. flavicollis, which only forms a narrow elongated spot on the breast in A. sylvaticus .
Their prevalence in Western Palearctic and common status in Western and Central Europe made them one of the model organisms to study post-glacial movement of mammals [22, 41]. Both species have traditionally been studied in a parasitological context, as one of the vectors of Borellia-carrying ticks Ixodes ricinus, who often feed on Apodemus [43, 58], tick-borne encephalitis virus  and hantaviruses [31, 46] and have been used as markers for environmental quality [36, 63]. Lastly, they have extra-autosomal chromosomes, called B chromosomes, with varied distribution among the populations  and suggested involvement in a variety of physiological phenomena, from cell division and development to immune response .
Previous studies on Apodemus typically employed a small number of microsatellite  and mtDNA markers [22, 38, 40, 41], which are insufficient to learn about the species’ population structure and admixture patterns in detail, or to identify loci under selection. In the absence of high-quality reference genome, which remains cost-prohibitive for complex genomes, whole-genome marker discovery enabled by restriction site-associated DNA sequencing presents a cost-effective method to study species on a population scale even with no previous genetic and genomic resources available .
Here we employ the double-digest restriction site-associated DNA sequencing (ddRAD-seq) to elucidate the genetic structure and connectivity of three populations of A. flavicollis and compare it to a population of A. sylvaticus in Poland. We demonstrate clear divergence between the two species and very low differentiation between populations of A. flavicollis. Our results provide the first estimates of population parameters in A. flavicollis based on thousands of loci, calculation of p-distance between the two Apodemus species, as well as a selection of loci enabling their accurate identification.
Sequencing and variant calling
The sequencing produced a total of 92741120 reads. The number of reads per individual varied from 346810 to 4157586, with an average of 1078385 reads per individual and median of 905786,5 (Supplementary Table S2). The best parameters for calling the stacks and variants for the entire dataset were: minimum number of identical, raw reads required to create a stack m = 2, number of mismatches allowed between loci for each individual M = 4 and number of mismatches allowed between loci when building the catalogue n = 5 (Supplementary Figure S1). The best parameters calculated for A. flavicollis samples only were: m = 2, M = 4 and n = 3 (Supplementary Figure S3). The coverage per sample ranged from 4.95x to 26.20x with an average of 10.13x and median of 9.32x for the entire dataset (Supplementary Figures S2 and S4).
SNPs and loci co-identification rates
Analysis of the duplicated samples showed that loci and allele misassignment rates were of similar magnitude, on average, between all pairs of duplicates. The duplicate pair F06-B02 showed the highest discrepancy between loci, of 10%, and also between alleles, of 8%. When only shared loci were included in the comparisons, all four sets of duplicates showed on average 0.5% ±0.2% SNPs called differently (Table 1).
Comparison of A. flavicollis and A. sylvaticus
The number of assembled loci per individual ranged from 46286 to 117366 (mean: 73711, median: 71395, standard deviation: 29917). 52494 loci passed the population filters established for species differentiation (see Methods, section "Variant calling and filtering"), representing 8,3% of the total 632063 loci included in the catalogue. Out of 158144 SNPs called, 60366 (38.1%) were removed after filtering for minor allele frequency (MAF) and 52298 (33%) were removed after failing the HWE test at p<0.05; further 35302 (22.3%) were removed due to a minimum mean depth lower than 20, leaving 10178 SNPs (6.6%) to be used in the downstream analyses (Fig. 1). PCA plot of the first two components (Fig. 2), accounting for 13.13% of the total variance, shows differentiation of the two species but also distinguish different populations of A. flavicollis.
Similarly, the phylogenetic tree shows A. sylvaticus as a separate clade to the three populations of A. flavicollis, with A. flavicollis from geographically closer regions (Białowieża and Haćki, 50 km) grouped closer than a population from Bory Tucholskie, 450 km away from Białowieża (Fig. 3). The A. sylvaticus and A. flavicollis clusters have high bootstrap value support (100% and 99% respectively).
We then investigated the suitability of the loci we identified on Polish populations to distinguish A. sylvaticus and A. flavicollis from other European populations. The genotyping of the extra 10 samples from each species (see Methods) produced 179763 SNPs. 62158 (34.58%) were removed after filtering for MAF and 69125 (38.45%) were removed after failing the HWE test at p<0.05; further 42054 (23.39%) were removed due to a minimum mean depth lower than 20 and 5203 (2.89%) were removed due to more than 5% missing data, leaving 1223 SNPs (0.68%) to be used in the downstream analyses.
The first axis of the PCA plot (Fig. 4) constructed from this data accounts for the 65.73% of the total variance and shows clear differentiation between the two species. All the A. flavicollis samples cluster with the Polish A. flavicollis samples, while all but Tunisian samples of A. sylvaticus cluster with the Polish samples of the same species. Tunisian A. sylvaticus appear as a separate cluster but still closer to the A. sylvaticus group. The catalogue of loci used for species identification is included in the Supplementary Materials, Section 6.
Genetic diversity and population structure of A. flavicollis
The number of assembled loci per individual in the Polish populations ranged from 46286 to 117366 (mean: 72738, median: 70592, stdev: 12575). 30722 loci passed the population filters established for population differentiation, representing and 4,43% of the total 691960 loci included in the catalog. Out of 63742 SNPs called, 31401 (49.26%) were removed after filtering for MAF and 10034 (15.74%) were removed after failing the HWE test at p<0.05. Further 9653 (15.14%) were removed due to a minimum mean depth lower than 20, leaving 12654 (19.85%) SNPs to be used in the downstream analyses (Fig. 1).
PCA plot (Fig. 5) shows differentiation between the three Polish A. flavicollis populations, with PC1 and PC2 cumulatively explaining 10.47% of the total variance. Haćki population shows larger diversity than the other populations, with some Haćki individuals closer to Białowieża individuals than to others from this location. Phylogenetic tree (Fig. 6) supports this pattern of differentiation. Bory Tucholskie and Haćki populations each form a cluster with a 100% of bootstrap support value, whereas Białowieża forms a third cluster with an 95% of bootstrap support. Białowieża and Bory Tucholskie population together form a large cluster with a 100% bootstrap support.
In the ADMIXTURE analysis, the lowest cross-validation errors  were always found for K = 3, indicating contribution of three ancestral populations (Fig. 7). Majority of samples from each of the populations show a single dominant component of ancestry with little contribution from other populations, with the exception of four individuals from Haćki, which show clear admixture of the Białowieża population.
Recognising that STRUCTURE-type analyses (on which ADMIXTURE is based) may be sensitive to the effects of uneven number of samples in compared groups , we repeated the ADMIXTURE analysis 10 times, each time randomly drawing the same number of individuals (n = 15) from each population. In all cases, the lowest cross-validation errors were found for K = 2, followed by K = 3 (Supplementary Figure S5). At even sampling, ADMIXTURE pattern found for K = 3 was the closest to the observed ecological and geographical distribution of the samples and closely matched our results when all samples were included (Supplementary Figure S6).
The patterns of heterozygosity highlight Haćki as the only population where the values of Ho is higher than He, where the FIS is negative (Table 2). As parameters such as number of private alleles, nucleotide diversity and heterozygosity can vary with sample size, we performed 100 calculations of the above parameters using random sampling of the same number of individuals (n = 15) from each population. The parameters showed similar relationships except for the number of private alleles (data not shown).
Fst values are consistently very low between all the populations, even though populations from Haćki and Bory Tucholskie show three-fold higher Fst values that for the other two pairs of populations (Table 3).
Finally, we calculated that the average p-distance between A. flavicollis and A. sylvaticus, based on 21377 shared loci, is 1.51% (standard deviation = 1.11%).
We then identified the top 117 most divergent loci between the species, which all had the divergence larger than 4.9% (The loci ID are provided in the Supplementary Table S3), and checked whether these loci alone allow for accurate assignment of samples to the two species. We constructed PCA plots from the Polish samples only and from the Polish, other European and Tunisian samples together. They demonstrate that while the 117 loci are sufficient to clearly assign Polish samples to the two species (Supplementary Figure S8), some uncertainty remains when we use these loci for the broader set of samples. Whereas all A. flavicollis samples do cluster together, A. sylvaticus samples do not form a clearly differentiated group (Supplementary Figure S9).
We also identified fixed loci, where all individuals within each species have identical sequences. There were 3526 such fixed loci for A. flavicollis and 5843 for A. sylvaticus. We then used 1273 of those loci that were shared among the two species and calculated that the average p-distance based on fixed differences is 0.97% (standard deviation = 0.94%).
RAD-sequencing approaches, including double-digest RAD-seq and its variants [6, 19, 42, 49, 50], have allowed a cost-effective discovery of thousands of genetic markers in both model and non-model organisms [21, 60], proving to be a transformative research tool in population genetics [8, 13, 24], phylogeography and phylogenetics [4, 23, 27, 57], marker development , linkage mapping studies , species differentiation  and detecting selection . However, despite the widespread use of this approach to marker discovery, only few studies have used RAD-seq in mammals [18, 30, 32, 44, 61]. Here, we have identified over 10000 markers in two closely related and common species of Apodemus in Western Palearctic, characterised the population structure of A. flavicollis and compared it to A. sylvaticus, for the first time providing estimates of the species divergence and population genetic parameters based on thousands of SNPs.
We have used four pairs of technical duplicates to check the accuracy of the RAD-seq genotyping based on the Poland protocol . The largest source of discrepancy in SNP calls between the duplicates is caused by unequal identification of loci: the difference in our case averaged approximately 10% (Table 1) and was similar to allele misindentification rates. However, when considering only shared loci between the duplicates, the discrepancy in SNP calls fell by over an order of magnitude to an average of 0.5%, indicating high accuracy and reliability of calls in once-defined shared loci. Our finding of loci calls being the major source of genotyping variability agrees with Mastretta et al. (2015), although our discrepancies are almost an order of magnitude smaller. Moreover, despite the differences in number of loci included in the analysis, each duplicated pair of samples clustered together with a 100% bootstrap values support and branch length equal to 0 on the phylogenetic tree (Fig. 6), indicating that the samples were identical. Overall, our finding reiterates the importance of the influence of stochastic events and imprecise size selection in the library preparations on genotyping calls . We note that some of these variables could be better controlled with automated size-selection approaches . Our findings also illustrate the usefulness of including technical replicates during library preparation.
Effect of group size
Permutations performed for the calculations of genetic diversity parameters (Table 4) have shown that with the exception of the number of private alleles, the results are comparable, regardless of the number of samples included per each population.
In the ADMIXTURE, we observed different optimal K depending on whether all samples were included in the analysis (K = 3) or a set of 15 randomly-chosen set of samples from each population (K = 2, although closely followed by K = 3). While the previously reported tendency of STRUCTURE-like analyses to produce ΔK = 2 does not apply in our case due to different method to select optimal number of clusters , we chose to use K = 3 for our analyses due to close match to the spatial and ecological locations from which our populations were sampled. The results obtained for K = 3 in the evenly-sampled dataset were similar to the clusters obtained for K = 3 with the complete dataset. ADMIXTURE plots for K = 2 to K = 5 are shown as Supplementary Figure S7.
The FST values calculated in this study between all three pairs of populations of A. flavicollis, based on 12654 SNPs, are consistently low and are not affected when we randomly draw the same number of individuals from each population to compute pairwise FST (Table 5). Previous studies of A. flavicollis populations in north-eastern Poland based on a small number of microsatellites showed similarly and consistently low values [15, 20], even though  suggested some population structure based on statistically significant differences between very low pairwise FST values.  also suggest large, broadly geographically defined clusters of A. flavicollis in north-eastern Poland that are separated by highly admixed individuals, but, again, FST between those clusters are as low as those reported by  and this study.
We would argue, based on a much larger set of markers reported here, that A. flavicollis has a negligible population structure across the entire area studied. Large number of markers nevertheless allows us to discover evidence for admixture of Białowieża population and Haćki (Fig. 7), further indicated by relatively high heterozygosity and negative FIS in this population. It is therefore intriguing that such a low differentiation occurs across hundreds of kilometres of varying landscape in a species that typically has a limited range of about 4 km2 and that suffers up to 86% winter mortality rate , which would lead to multiple bottlenecks and drift-driven population differentiation. With this in mind, our data suggests a much larger dispersal ability of the species, a much better connectivity between populations, or both.
Both low overall FST and moderate heterozygosity suggest it would be worthwhile to conduct a genome-wide scan for selection using FST as a metrics of local genomic differentiation to identify geographically local regions under selection. This, however, is not yet possible given the lack of high-quality reference genome for Apodemus and unknown synteny to the available genome of Mus musculus.
Divergence and differentiation of A. flavicollis and A. sylvaticus
Given that accurate identification of the two species using morphological characters is problematic, especially in their southern range , a large collection of markers identified in this study allowed us to create a catalogue of 632060 loci that allow clear differentiation between species. This identification is somewhat biased, as the catalogue was built using many more samples of A. flavicollis than A. sylvaticus (72 vs 10) and both from a relatively limited geographical range. Nevertheless, it allowed for accurate assignment of A. flavicollis samples and to a large degree of A. sylvaticus, as we demonstrated on a set of 20 independent samples from other European countries and Tunisia (Figure 4). The imperfect clustering of some European A. sylvaticus samples on Figure 4 may reflect the higher genetic differentiation of population from the edge of the range of the species: samples near the centre of the PC1 axis on Figure 4 (green dots in a circle) all come from Tunisia. Given the wide distribution of both species in Western Palearctic, a more representative sample from both species from a broader geographic range would likely provide more accurate set of markers for their identification.
Finally, we calculated the nucleotide divergence between the two species, based on 21377 shared loci, which is 1.51%. Considering a divergence time between A. flavicollis and A. sylvaticus estimated from archeological data of 4 Mya , the evolutionary rate is 0.0019 substitutions per site per million of years. While this estimate of sequence divergence level is in broad agreement with calculations based on mitochondrial 12S rRNA, IRBP and Cytochrome b genes , we note that Aghova et al. (2018)  provided an updated timing of a split between A. sylvaemus and A. mystacinus, which, at 9.6Mya, is 2Mya older than previously suggested by Michaux et al. (2002) . If we use that date as a reference and move the presumed split between A. flavicollis and A. sylvaticus by 2.6Mya, the estimated evolutionary rate from our data would be 0.0011. In both cases, our calculation is likely an underestimate, as we only used shared loci to calculate divergence and did not include the potential impact of insertion/deletion events, which can significantly affect the total genomic divergence between species [9, 35]. Highly divergent sequences would have been identified as different loci, and would not be compared to their true homologous sequences.
We have successfully applied the ddRAD-seq approach to discover tens of thousands of SNPs in wide-spread and common mammalian species of A. flavicollis and A. sylvaticus. The high resolution data obtained here allowed us to delineate geographically close populations, including identifying admixture between them, but suggest that A. flavicollis effectively forms a single population in an entire sampling area that spans 500 km in the W-E direction. Comparing A. flavicollis and A. sylvaticus, we have calculated their divergence and identified a set of genomic loci that enable effective molecular identification of the species. We anticipate that with the development of further whole-genome resources, Apodemus, thanks to its common status, broad geographic range and long history of ecological observations, will become an excellent model species for evolutionary and ecological research in the genomic era.
Sample collection and dNA extraction
Eighty two individuals (10 Apodemus sylvaticus and 72 Apodemus flavicollis) from four locations in northern Poland spanning 500 km were trapped in 2015 (Fig. 8). A. flavicollis were collected in Białowieża (E23.8345814, N52.7231935), an oak-lime-hornbeam forest (n = 35), Bory Tucholskie (E17.5160265, N53.7797608), in an oak-lime-hornbeam and pine forest (n = 23) and Haćki (E23.1793284, N52.834369), in a xerothermic meadow (n = 14). A. sylvaticus were trapped in Białowieża (E21.3778496, N53.2089113) in a dry pine forest (n = 5) and in Bory Tucholskie, mainly in a pine forest (n = 5) (Supplementary Table S1). While A. flavicollis are present in all sampled locations, there have been no trappings of A. sylvaticus in Białowieża for the last 20 years, despite Białowieża being within the European range of this species (Dr Karol Zub, personal communication). The sampling procedures were approved by the Local Ethical Commission on Experimentation on Animals in Białystok, Poland, under permission number 2015/99. All animals were released after sample collection.
Tail clippings were collected, preserved in ≥95% ethanol and stored at −20∘C until DNA extraction. The tissues were digested by incubating at 55∘C overnight with lysis buffer (10mM Tris, 100mM NaCl, 10mM EDTA, 0.5% SDS) and proteinase K (20mg/ml). Subsequently, potassium acetate and RNAse A were used to remove protein and RNA contamination. Three ethanol washes were performed using Sera-Mag SpeedBeads solution (GElifesciences, Marlborough, MA, USA). The quality and integrity of the DNA was tested in a 2% agarose gel. Twenty-fold dilutions of the samples were used to measure the DNA concentration using Quant-iT PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA, USA) and concentration of each sample was then normalised to 10 ng/ μl in 20 μl volume. Four samples were used as technical duplicates (F06-B02, G02-D01, H11-G06, F12-A12). Technical duplicates had the same DNA but were digested and ligated to barcodes independently.
ddRAD-seq library preparation
ddRAD-seq library was prepared following the protocol from , adapted to a different combination of enzymes. Briefly, genomic DNA was digested in a 20 μl reaction with CutSmart® buffer, 8 units of SbfI and 8 units of HF-MseI (New England Biolabs, Frankfurt am Main, Germany). Digestion was performed at 37∘C for 2 hours. Enzymes were inactivated at 65∘C for 20 minutes and the reactions were kept at 8∘C. Adapter ligation was performed at 22∘C for 2 hours and the ligase was inactivated by incubating the samples at 65∘C for 20 minutes. Samples were cooled down to 8∘C and multiplexed by combining 5 μl of each sample. P1 adapters contained barcodes with a length between 5 and 10 bp.
PCR amplification was conducted in 25 μl with 1 μl of each primer (IlluminaF_PE: AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT and IlluminaR_PE: CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAA) at 10mM, 0.5 μl of 10 mM dNTPs, 13.25 μl of PCR-grade water, 5 μl of 5x Phusion HF Buffer, 0.25 μl of Phusion DNA Polymerase (New England Biolabs, Frankfurt am Main, Germany) and 4 μl of the multiplexed DNA. After an initial denaturation step of 30s at 98∘C, PCR reaction was carried out for 12 cycles (10s at 98∘C, 20s at 58∘C and 15s at 72∘C). Final elongation step was performed at 72∘C for 5 minutes.
PCR products were loaded into a single lane on a 1% agarose gel with 100 bp DNA ladder (New England Biolabs, Frankfurt am Main, Germany). Fragments between 200 and 500 bp were cut from the gel with a scalpel and purified using the QIAquick gel extraction kit (QIAgen, Hilden, Germany), followed by the second cleanup step with Sera-Mag SpeedBeads (GElifesciences, Marlborough, MA. USA). Sizing, quantification and quality control of the DNA was performed using Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) before paired-end sequencing on an Illumina HiSeq 3500 with 150 bp read length.
Processing of RAD-tags
Sequences were analysed with Stacks version 1.48 . Samples were demultiplexed using process_radtags allowing no mismatches in barcodes and cutting sites. Sequences with uncalled bases and low quality scores were removed and all reads were trimmed to 141 bp. The four files generated per sample by process_radtags were concatenated using a custom bash script. The best parameters for building and calling SNPs de novo, using denovo_map, were calculated following  approach, using either samples from both species or only from A. flavicollis. Secondary reads were not used to call haplotypes in denovo_map (option -H).
SNPs and loci co-identification rates
We estimated the loci and SNP co-identification rates by analysing a set of four samples that were prepared and sequenced in duplicates. Sequences for 52494 loci from both species, were extracted using –fasta_samples option from the population package in Stacks. We extracted sequences for each of the duplicated samples with a custom script and calculated co-identification rates as described by . Briefly, the locus misassignment rate is the percentage of unidentified loci, calculated by dividing the number of loci found only in one of the duplicates by the total number of loci in each sample. The allele misassignment rate is the percentage of mismmatches between the IUPAC consensus sequences between homologous loci from each pair of duplicates. Finally, the two SNP error rates: the percentage of different SNPs called in each of the duplicated samples using either all 10178 SNPs or using the SNPs called without missing data between duplicate samples excluded (see Table 1).
Variant calling and filtering
We combined the data from A. sylvaticus and A. flavicollis to establish species differentiation and then filtered the SNPs using the population package from Stacks  and VCFtools . We kept SNPs from the loci present in the 80% of the individuals in each species (p=1, r=0.8) and excluded SNPs with minor allele frequencies MAF<0.05 and which deviated from the Hardy-Weinberg equilibrium (HWE) at P<0.05. We also removed sites with mean depth values lower than 20. We manually modified the chromosome numbers in the vcf file to input it into SNPhylo , which we used to build the tree. We set a missing rate (-M) of 1, minor allele frequencies (-m) of 0, linkage disequilibrium threshold (-l) of 1 and the -r option to skip the step of removing low quality data. Confidence values were estimated using 1000 bootstrap replicates. The root was manually fixed to separate both species. Principal Component Analysis (PCA) was performed using the R package Adegenet  (Fig. 3).
The set of divergent loci identified between the two species in Polish samples was tested for its ability to differentiate an extra set of samples from other locations in Europe and Tunisia. Ten A. flavicollis (2 samples from Austria, 5 from Lithuania and 3 from Romania) and 10 samples of A. sylvaticus (4 samples from Wales, 3 from Tunisia and 3 from Scotland) were kindly provided by Dr Jeremy Herman, National Museums Scotland, Dr Johan Michaux, University of Liege and Dr Karol Zub, Mammal Research Institute of the Polish Academy of Sciences (MRI) (Supplementary Table S4). We considered all 20 test samples as a different group from Polish A. sylvaticus and A. flavicollis for SNP calling. We kept SNPs from the loci present in the 80% of the individuals in each group (p=1, r=0.8) and excluded SNPs with minor allele frequencies MAF<0.05, SNPs which deviated from the Hardy-Weinberg equilibrium (HWE) at P<0.05, sites with mean depth values lower than 20 and with more than 5% of missing data.
To analyse genetic diversity and population connectivity within A. flavicollis, we analysed the three populations (Bory Tucholskie, Białowieża and Haćki) separately (p=3, r=0.8), while keeping the other parameters as described above. We note that absence of filtering out SNPs in Hardy-Weinberg disequlibrium did not markedly affect the reported pairwise Fst values (data not shown). Due to the lack of outgroup, a mid-point root was chosen in the phylogenetic tree. Individual ancestries were estimated following a maximum likelihood approach with ADMIXTURE , after conversion of the VCF file to ped with plink version 1.9 [12, 55]. ADMIXTURE analysis was run for each of K=1 to K=5, each using 10 different seeds. Weighted (Weir-Cockerham) Fst was calculated with VCFtools v0.1.13. Heterozygosity, Pi and Fis were calculated with the population package from Stacks .
To calculate the p-distance between the two species, first a set of common loci was extracted with a custom script and the strict consensus sequences for each species were calculated with Consensus.pl script  with no threshold parameter set, such that the most common nucleotide was set as a consensus. p-distance was then calculated using a custom R script that counts the number of differences between pairs of sequences (one from each species, for each of the 21377 loci) (Supplementary Materials, Section 9). The same set of scripts was also used to select loci with fixed differences within each species and then to calculate the p-distance between them.
Availability of data and materials
The sequence data, in the format of the output from the Stacks 1.48 proces_radtags (four files per sample), have been deposited in EBI SRA repository under accession number PRJNA554851. The scripts used in the analysis are deposited at GitHub.com, with links available in the Supplementary Materials.
double-digest restriction site-associated DNA sequencing
minor allele frequency
principal component analysis
single nucleotide polymorphism
variant call format
Aghová T, Kimura Y, Bryja J, Dobigny G, Granjon L, Kergoat GJ. Fossils know it best: Using a new set of fossil calibrations to improve the temporal phylogenetic framework of murid rodents (Rodentia: Muridae). Mol Phylogenet Evol. 2018; 128:98–111.
Alexander DH, Lange K. Enhancements to the admixture algorithm for individual ancestry estimation. BMC bioinformatics. 2011; 12(1):246.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009; 19(9):1655–1664.
Alter SE, Munshi-South J, Stiassny MLJ. Genome wide snp data reveal cryptic phylogeographic structure and microallopatric divergence in a rapids-adapted clade of cichlids from the congo river. Mol Ecol. 2017; 26(5):1401–1419.
Andrews KR, Good JM, Miller MR, Luikart G, Hohenlohe PA. Harnessing the power of radseq for ecological and evolutionary genomics. Nat Rev Genet. 2016; 17(2):81.
Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA. Rapid snp discovery and genetic mapping using sequenced rad markers. PloS one. 2008; 3(10):e3376.
Baxter SW, Davey JW, Johnston JS, Shelton AM, Heckel DG, Jiggins CD, Blaxter. ML. Linkage mapping and comparative genomics using next-generation rad sequencing of a non-model organism. PloS one. 2011; 6(4):e19315.
Blanco-Bercial L, Bucklin A. New view of population genetics of zooplankton: Rad-seq analysis reveals population structure of the north atlantic planktonic copepod centropages typicus. Mol Ecol. 2016; 25(7):1566–1580.
Britten RJ. Divergence between samples of chimpanzee and human dna sequences is 5%, counting indels. Proc Natl Acad Sci. 2002; 99(21):13633–13635.
Bugarski-Stanojević V, Blagojević J, Adnadjević T, Jovanović V, Vujošević M. Identification of the sibling species apodemus sylvaticus and apodemus flavicollis (rodentia, muridae)—comparison of molecular methods. Zool Anz J Comp Zool. 2013; 252(4):579–587.
Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH. Stacks: building and genotyping loci de novo from short-read sequences. G3: Genes, genomes, genetics. 2011; 1(3):171–182.
Chang CC, Chow CC, Tellier LC, Vattikuti s, Purcell sm. Second-generation plink: rising to the challenge of larger and richer datasets. Gigascience. 2015; 4(1):7.
Cromie GA, Hyma KE, Ludlow CL, Garmendia-Torres C, Gilbert TL, May P, Huang AA, Dudley AM, Fay JC. Genomic sequence diversity and population structure of saccharomyces cerevisiae assessed by rad-seq. G3: Genes, Genomes, Genet. 2013; 3(12):2163–2171.
Cull B, Vaux AGC, Ottowell LJ, Gillingham EL, Medlock JM. Tick infestation of small mammals in an english woodland. J Vector Ecol. 2017; 42(1):74–83.
Czarnomska SD, Niedziałkowska M, Borowik T, Jędrzejewska B. Regional and local patterns of genetic variation and structure in yellow-necked mice-the roles of geographic distance, population abundance, and winter severity. Ecol Evol. 2018; 8:8171–8186.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al.The variant call format and vcftools. Bioinformatics. 2011; 27(15):2156–2158.
Darvish J, Mohammadi Z, Ghorbani F, Mahmoudi A, Dubey S. Phylogenetic relationships of apodemus kaup, 1829 (rodentia: Muridae) species in the eastern mediterranean inferred from mitochondrial dna, with emphasis on iranian species. J Mamm Evol. 2015; 22(4):583–595.
Fernández R, Schubert M, Vargas-Velázquez AM, Brownlow A, Víkingsson GA, Siebert U, Jensen LF, Øien N, Wall D, Rogan E, Mikkelsen B. A genomewide catalogue of single nucleotide polymorphisms in white-beaked and Atlantic white-sided dolphins. Mol Ecol Resour. 2016; 16(1):266–276.
Franchini P, Parera DM, Kautt AF, Meyer A. quaddrad: a new high-multiplexing and pcr duplicate removal ddrad protocol produces novel evolutionary insights in a nonradiating cichlid lineage. Mol Ecol. 2017; 26(10):2783–2795.
Gortat Tt, Gryczyńska-Siemiątkowska A, Rutkowski R, Kozakiewicz A, Mikoszewski A, Kozakiewicz M. Landscape pattern and genetic structure of a yellow-necked mouse apodemus flavicollis population in north-eastern poland. Acta Theriol. 2010; 55(2):109–121.
Hammerman NM, Rivera-Vicens RE, Galaska MP, Weil E, Appledoorn RS, Alfaro M, Schizas NV. Population connectivity of the plating coral agaricia lamarcki from southwest puerto rico. Coral Reefs. 2018; 37(1):183–191.
Herman JS, Jóhannesdóttir FJ, Jones EP, McDevitt AD, Michaux JR, White TA, Wójcik JN, Searle JB. Post-glacial colonization of europe by the wood mouse, apodemus sylvaticus: evidence of a northern refugium and dispersal with humans. Biol J Linn Soc. 2017; 120(2):313–332.
Hipp AL, Eaton DAR, Cavender-Bares J, Fitzek E, Nipper R, Manos PS. A framework phylogeny of the american oak clade based on sequenced rad data. PLoS One. 2014; 9(4):e93975.
Hohenlohe PA, Day MD, Amish SJ, Miller MR, Kamps-Hughes N, Boyer MC, Muhlfeld CC, Allendorf FW, Johnson EA, Luikart G. Genomic patterns of introgression in rainbow and westslope cutthroat trout illuminated by overlapping paired-end rad sequencing. Mol Ecol. 2013; 22(11):3002–3013.
Hughes J. Sequence-manipulation. GitHub Repository. 2011. https://github.com/josephhughes/Sequencemanipulation.
Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, Andrew RL. The k= 2 conundrum. Mol Ecol. 2017; 26(14):3594–3602.
Jeffries DL, Copp GH, Handley LL, Olsén KH, Sayer CD, Hänfling B. Comparing radseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the crucian carp, carassius carassius, l. Mol Ecol. 2016; 25(13):2997–3018.
Jojić V, Bugarski-Stanojević V, Blagojević J, Vujošević M. Discrimination of the sibling species apodemus flavicollis and a. sylvaticus (rodentia, muridae). Zool Anz J Comp Zool. 2014; 253(4):261–269.
Jombart T. adegenet: a r package for the multivariate analysis of genetic markers. Bioinformatics. 2008; 24(11):1403–1405.
Knowles LL, Massatti R, He Q, Olson LE, Lanier HC. Quantifying the similarity between genes and geography across alaska’s alpine small mammals. J Biogeogr. 2016; 43(7):1464–1476.
Kolodziej M, Melgies A, Joniec-Wiechetek J, Michalski A, Nowakowska A, Pitucha G, Niemcewicz M. First molecular characterization of dobrava-belgrade virus found in apodemus flavicollis in poland. Ann Agric Environ Med. 2018; 25(2):368–373.
Lanier HC, Massatti R, He Q, Olson LE, Knowles. LL. Colonization from divergent ancestors: glaciation signatures on contemporary patterns of genomic variation in collared pikas (ochotona collaris). Mol Ecol. 2015; 24(14):3688–3705.
Lapinski AG, Pavlenko MV, Solovenchuk LL, Gorbachev VV. Some limitations in the use of the mitochondrial dna cytb gene as a molecular marker for phylogenetic and population-genetic studies by the example of the apodemus genus. Russ J Genet Appl Res. 2016; 6(1):84–90.
Lee T-H, Guo H, Wang X, Kim C, Paterson AH. Snphylo: a pipeline to construct a phylogenetic tree from huge snp data. BMC Genomics. 2014; 15(1):162.
Li W-H, Tanimura M, Sharp PM. An evaluation of the molecular clock hypothesis using mammalian dna sequences. J Mol Evol. 1987; 25(4):330–342.
Martiniaková M, Omelka R, Jancová A, Stawarz R, Formicki G. Heavy metal content in the femora of yellow-necked mouse (apodemus flavicollis) and wood mouse (apodemus sylvaticus) from different types of polluted environment in slovakia. Environ Monit Assess. 2010; 171(1-4):651–660.
Mastretta-Yanes A, Arrigo N, Alvarez N, Jorgensen TH, Piñero D, Emerson BC. Restriction site-associated dna sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Mol Ecol Resour. 2015; 15(1):28–41.
René Michaux J, Magnanou E, Paradis E, Nieberding C, Libois R. Mitochondrial phylogeography of the woodmouse (apodemus sylvaticus) in the western palearctic region. Mol Ecol. 2003; 12(3):685–697.
Michaux JR, Chevret P, Filippucci M-G, Macholan M. Phylogeny of the genus apodemus with a special emphasis on the subgenus sylvaemus using the nuclear irbp gene and two mitochondrial markers: cytochrome b and 12s rrna. Mol Phylogenet Evol. 2002; 23(2):123–136.
Michaux JR, Libois R, Paradis E, Filippucci M-G. Phylogeographic history of the yellow-necked fieldmouse (apodemus flavicollis) in europe and in the near and middle east. Mol Phylogenet Evol. 2004; 32(3):788–798.
Michaux JR, Libois R, Filippucci MG. So close and so different: comparative phylogeography of two small mammal species, the yellow-necked fieldmouse (apodemus flavicollis) and the woodmouse (apodemus sylvaticus) in the western palearctic region. Heredity. 2005; 94(1):52–63.
Miller MR, Dunham JP, Amores A, Cresko WA, Johnson EA. Rapid and cost-effective polymorphism identification and genotyping using restriction site associated dna (rad) markers. Genome Res. 2007; 17(2):240–248.
Mlera L, Bloom ME. The role of mammalian reservoir hosts in tick-borne flavivirus biology. Front Cell Infect Microbiol. 2018; 8:298.
Moura AE, Kenny JG, Chaudhuri R, Hughes MA, Welch AJ, Reisinger RR, de Bruyn PJN, Dahlheim ME, Hall N, Hoelzel AR. Population genomics of the killer whale indicates ecotype evolution in sympatry involving both selection and drift. Mol Ecol. 2014; 23(21):5179–5192.
Pante E, Abdelkrim J, Viricel A, Gey D, France SC, Boisselier M-C, Samadi S. Use of rad sequencing for delimiting species. Heredity. 2015; 114(5):450.
Papa A, Rogozi E, Velo E, Papadimitriou E, Bino S. Genetic detection of hantaviruses in rodents, albania. J Med Virol. 2016; 88(8):1309–1313.
Paris JR, Stevens JR, Catchen JM. Lost in parameter space: a road map for stacks. Methods Ecol Evol. 2017; 8(10):1360–1373. https://0-doi-org.brum.beds.ac.uk/10.1111/2041-210x.12775.
Pegadaraju V, Nipper R, Hulke B, Qi L, Schultz Q. De novo sequencing of sunflower genome for snp discovery using rad (restriction site associated dna) approach. BMC Genomics. 2013; 14(1):556.
Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest radseq: an inexpensive method for de novo snp discovery and genotyping in model and non-model species. PloS One. 2012; 7(5):e37135.
Poland JA, Rife TW. Genotyping-by-sequencing for plant breeding and genetics. Plant Genome. 2012; 5(3):92–102.
Poland JA, Brown PJ, Sorrells ME, Jannink J-L. Development of high-density genetic maps for barley and wheat using a novel two-enzyme genotyping-by-sequencing approach. PloS One. 2012; 7(2):e32253.
Pucek Z, Stachurska J, Bibrich M, et al.Keys to vertebrates of poland: mammals. 1981.
Pucek Z, Jędrzejewski W, Jędrzejewska B, Pucek M. Rodent population dynamics in a primeval deciduous forest (białowieża national park) in relation to weather, seed crop, and predation. Acta Theriol. 1993; 38(2):199–232.
Puechmaille SJ. The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Mol Ecol Resour. 2016; 16(3):608–627.
Purcell S, Chang CC. Plink 1.9 package. 2018. https://www.cog-genomics.org/plink/1.9/. Accessed 27 Feb 2019.
Rajičić M, Romanenko SA, Karamysheva TV, Blagojević J, Adnađević T, Budinski I, Bogdanov AS, Trifonov VA, Rubtsov NB, Vujošević M. The origin of b chromosomes in yellow-necked mice (apodemus flavicollis)—break rules but keep playing the game. PloS One. 2017; 12(3):e0172704.
Reitzel AM, Herrera S, Layden MJ, Martindale MQ, Shank TM. Going where traditional markers have not gone before: utility of and promise for rad sequencing in marine invertebrate phylogeography and population genomics. Mol Ecol. 2013; 22(11):2953–2970.
Richter D, Schlee DB, Matuschka F-R. Reservoir competence of various rodents for the lyme disease spirochete borrelia spielmanii. Appl Environ Microbiol. 2011; 77(11):3565–3570. https://0-doi-org.brum.beds.ac.uk/10.1128/aem.00022-11.
Rico A, Kindlmann P, Sedláček F. Can the barrier effect of highways cause genetic subdivision in small mammals?. Acta Theriol. 2009; 54(4):297–310.
Rodriguez-Ezpeleta N, Álvarez P, Irigoien X. Genetic diversity and connectivity in maurolicus muelleri in the bay of biscay inferred from thousands of snp markers. Front Genet. 2017; 8:195.
Shafer A, Peart CR, Tusso S, Maayan I, Brelsford A, Wheat CW, Wolf JBW. Bioinformatic processing of rad-seq data dramatically impacts downstream population genetic inference. Methods Ecol Evol. 2017; 8(8):907–917.
Shultz AJ, Baker AJ, Hill GE, Nolan PM, Edwards SV. Snp s across time and space: population genomic signatures of founder events and epizootics in the house finch (haemorhous mexicanus). Ecol Evol. 2016; 6(20):7475–7489.
Velickovic M. Measures of the developmental stability, body size and body condition in the black-striped mouse (apodemus agrarius) as indicators of a disturbed environment in northern serbia. Belg J Zool. 2007; 137(2):147.
Vujošević M, Rajičić M, Blagojević J. B chromosomes in populations of mammals revisited. Genes. 2018; 9(10):487.
The authors wish to thank Jeremy Herman (National Museums Scotland) and Johan Michaux (Université de Liège) for providing the European/Tunesian samples for testing our catalog of SNPs. MLMC and JB also wish to acknowledge the use of the Orion High Performance Computing cluster at the School of Applied Sciences, University of Huddersfield.
This project was supported by the University of Huddersfield, the Friedrich Miescher Laboratory of the Max Planck Society and Microsoft Azure Research Award CRM:0518338.
Ethics approval and consent to participate
The sampling procedures were approved by the Local Ethical Commission on Experimentation on Animals in Białystok, Poland, under permission number 2015/99.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Martin Cerezo, M.L., Kucka, M., Zub, K. et al. Population structure of Apodemus flavicollis and comparison to Apodemus sylvaticus in northern Poland based on RAD-seq. BMC Genomics 21, 241 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-020-6603-3
- population structure
- Apodemus flavicollis
- Apodemus sylvaticus