Identifying the minimum number of microsatellite loci needed to assess population genetic structure: A case study in fly culturing.

Small, isolated populations are constantly threatened by loss of genetic diversity due to drift. Such situations are found, for instance, in laboratory culturing. In guarding against diversity loss, monitoring of potential changes in population structure is paramount; this monitoring is most often achieved using microsatellite markers, which can be costly in terms of time and money when many loci are scored in large numbers of individuals. Here, we present a case study reducing the number of microsatellites to the minimum necessary to correctly detect the population structure of two Drosophila nigrosparsa populations. The number of loci was gradually reduced from 11 to 1, using the Allelic Richness (AR) and Private Allelic Richness (PAR) as criteria for locus removal. The effect of each reduction step was evaluated by the number of genetic clusters detectable from the data and by the allocation of individuals to the clusters; in the latter, excluding ambiguous individuals was tested to reduce the rate of incorrect assignments. We demonstrate that more than 95% of the individuals can still be correctly assigned when using eight loci and that the major population structure is still visible when using two highly polymorphic loci. The differences between sorting the loci by AR and PAR were negligible. The method presented here will most efficiently reduce genotyping costs when small sets of loci ("core sets") for long-time use in large-scale population screenings are compiled.

Introduction

Culturing organisms in the laboratory is a frequently used method in many biological fields such as ecology, conservation biology, and evolutionary biology. With careful management, captive populations can maintain high genetic diversity over many generations. However, when population size is small, genetic variation may be lost due to genetic drift. A depletion of genetic variability can increase homozygosity within the population, which in turn may cause lower viability and fecundity, an effect termed inbreeding depression. Rapid identification of a loss of genetic variation within and/or among populations and generations is crucial for successful laboratory culturing.

One method to rapidly identify changes in genetic diversity is monitoring using molecular markers like microsatellites, single nucleotide polymorphisms, or whole-genome fingerprints. Such markers provide insight into genetic variation and evolutionary processes and allow identification of specimens and their populations of origin. A limiting factor for the use of genetic monitoring is that these techniques are complex, expensive, and time consuming. Here, we test a method to simplify genetic monitoring by a stepwise reduction of the number of molecular markers to the minimum needed to still detect the relevant signature of population structure.

In testing this method, we used microsatellite markers. Microsatellites were discovered in the early 1980s and are short and highly variable nucleotide tandem repeats in DNA sequences. The formerly tedious isolation of microsatellites for new species has, with the advent of next generation sequencing, become routine. Once isolated, these codominant markers are easily amplified by Polymerase Chain Reaction (PCR) and analyzed using capillary electrophoresis or, more recently, Illumina sequencing. Their mutation rate of 10−6 to 10−2 events per locus and generation qualifies them for the detection of drift effects within very few generations. While in the last years single nucleotide polymorphisms (SNPs) have gained in popularity relative to microsatellites due to their greater abundance in genomes and lower genotyping error rates, for analysis of, for instance, population size dynamics and population structure, microsatellites are on a per-locus basis two to 20 times more informative than SNPs. For many study designs, microsatellites provide an appropriate information density at substantially lower costs than genome-based approaches. Thus, microsatellites are considered to remain important genetic markers for years to come. While next generation sequencing based isolation of microsatellites now facilitates the use of arbitrarily large numbers of loci and thus very precise detection of molecular variation, in a scientific world with often limited resources the other way, that is, a cost efficient use of genetic markers, will often remain desirable. Fundamental factors for the expenses of a microsatellite-based study are the number of loci used and the resulting trade-off between accuracy and costs. In a recent study, we estimated the costs per locus and individual at ca. 3 Euro, based on multiplexing of three loci into one capillary electrophoresis run. Reducing the number of loci would thus diminish a study's costs in terms of money and time.

In detail, the purpose of this study was to evaluate a method determining how few loci are needed to still detect the short-term changes in genetic structure in a population threatened by loss of diversity. We used individuals of two populations of the cool-climate mountain fly Drosophila nigrosparsa Strobl, 1898. Native to the Alps and threatened by climate change, the fly has been in the focus of intense research including selection experiments for which it had to be kept in the laboratory. Permanent laboratory stocks were established in 2012 and genotyped with the full set of 11 variable microsatellite loci available for this species in Generation 0 and Generation 5. Then, the performance of various subsets of the full microsatellite set was evaluated. The criteria in evaluating the performance of a set of loci were (a) the number of genetic clusters retrieved and (b) the assignment of individuals to clusters. In applying criterion (a), we considered as optimum number of loci their lowest number at which still the correct number of clusters (i. e., the number of clusters retrieved with the full set of 11 loci) was identified. In applying criterion (b), the particular aim of a research has to be set; in the extremes, this can be either a maximum accuracy in assigning individuals to clusters or a maximum number of individuals that are assigned. We explored the trade-off that applies to the two aims, in that achieving the first aim comes at the cost of excluding individuals (here termed exclusion rate) and the second at the cost of losing accuracy in assigning individuals to clusters (incorrect assignment rate). In sorting loci for being reduced, we assessed both the Allelic Richness (AR) and the Private Allelic Richness (PAR) made accessible by each locus. AR is a proxy for the number of alleles per locus, a measure of genetic diversity that indicates the potential for adaptability and persistence of a population. PAR is a proxy for the number of alleles unique to a population and is a simple measure of genetic distinctiveness. In detail, we first removed the locus with the lowest values of AR or PAR and retained to the end the locus with the highest values.

Materials and methods

Individuals of two natural populations of Drosophila nigrosparsa Strobl, 1898 were sampled at Kaserstattalm (Austria, 11.29°E 47.13°N, 2000 m above sea level, a.s.l.) and at Pfitscherjoch (Italy, 11.68°E 46.98°N, 2000 m a.s.l.) using fermented banana baits from July to August 2012. From each population, 100 males and 100 females were used to create a laboratory population, Kaserstattalm (henceforth K) and Pfitscherjoch (henceforth P). After oviposition, 31 field-caught female flies each of populations K and P were fixed in 96% ethanol and stored at −20°C for molecular analysis (Generation 0, henceforth K0 and P0). The laboratory populations K and P were kept and cultivated in quarantine for four generations to adapt flies to laboratory conditions and to eliminate potential diseases. Eggs were cultivated on malt medium (10 g agar, 1000 ml deionised water, 15 g dried yeast, 100 g ground malt, 3 ml methyl-4-hydroxybenzoate, 3.6 ml propionic acid, and 50 g ground maize), and adults were kept in inverted, transparent 0.3-l plastic cups with ventilation holes on petri-dishes with grape-juice agar (30 g agar, 1000 ml deionized water, 334 ml grape juice, 3.4 ml methyl-4-hydroxybenzoate, and 34 g sucrose). Flies were reared in environmental test chambers (MLR-352H-PE, Panasonic Healthcare Co., Japan) mimicking the diurnal temperature variation at 2000 m a.s.l. as well as the temperature at which the fly was found in the field at 2000 m a.s.l. in Tyrol in summer with a light:dark period of 16:8 hours and a humidity of 70% (Table S1). The fifth generation of laboratory populations K and P was randomly separated into eight lines. After oviposition, 31 females of each of the newly established lines were fixed in 96% ethanol and stored at −20°C for molecular analysis (Generation 5, henceforth K5 and P5).

DNA of individual flies was extracted using the Sigma GenElute Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich, St. Louis, USA) according to the manufacturer's instructions. Eleven species-specific microsatellite loci (DN16, DN31, DN35, DN36, DN37, DN39, DN40, DN41, DN45, DN48, DN49) were amplified via M13-tailed labelling. Amplifications were carried out in 5 µl reaction volume containing 1 × reaction buffer (Bioline, London, UK), 0.2 µM M13 primer, 0.02 µM M13-tailed forward primer, 0.2 µM reverse primer, 0.125 U MyTaq DNA polymerase (Bioline, London, UK), and 0.5 µl DNA extract on a UnoCycler 1200 (VWR, Radnor, USA). PCR conditions were 94°C for 2 min, followed by 35 cycles of 94°C for 30 s, Ta for 45 s, and 72°C for 1 min, followed by 72°C for 10 min. Locus-specific Ta was 48°C (DN16, DN36, DN40) and 55°C (DN31, DN35, DN37, DN39, DN41, DN45, DN49, DN55). Amplification success was checked by agarose gel electrophoresis. Capillary electrophoresis was performed by a commercial provider (CRC Sequencing Facility, Chicago, USA) using an ABI 3130 instrument (Applied Biosystems, Foster City, USA). The resulting traces were visualized using PeakScanner software v1.0 (Applied Biosystems, Foster City, USA) and scored manually. In total, microsatellite profiles of 558 individuals were generated (K0 = 31, K5 = 248, P0 = 31, P5 = 248).

Deviations from Hardy-Weinberg-Equilibrium (HWE), pairwise FST, and Analyses of Molecular Variance (AMOVA) were computed in GenAlEx v6.4140. In the AMOVA, the generations (G0, G5) and the populations (K, P) were used as hierarchy levels in multiple combinations. Linkage disequilibrium (LD) was calculated using Arlequin v3.5 with 10,000 iterations. For HWE and LD, Bonferroni-Holm corrections for multiple testing were performed.

To detect the genetic variation of K and P and to asses a potential effect of husbandry (loss of genetic variation from Generation 0 to Generation 5), AR and PAR were calculated for each population and generation using HP-Rare v1.0 assuming 52 genes, and significance (α = 0.05) was assessed by f- and two-sided t-tests in Excel 2013 (Microsoft Corp., Redmond, USA). HP-Rare uses rarefaction of alleles to compensate for differences in sample size, as larger samples are expected to display higher allele numbers.

STRUCTURE v2.3.3 was used to identify the population structure of the complete data set. The admixture model with correlated allele frequencies was chosen as recommended for faint population structures. The number of clusters (K) assumed was set to [1, 8], and each value of K was run 10 times. The Markov Chain was run for 20,000 generations burnin and 180,000 generations data collection. The optimum K was calculated using the method described by Evanno, Regnaut & Goudet (2005).

For further assessment of the minimum number of loci needed to distinguish between K and P flies, the individuals of Generation 5 were used. STRUCTURE analysis and search for best K were performed as described above with 496 individuals (K5 = 248, P5 = 248). The loci were then sorted by the mean AR value of the two populations; the locus with the lowest AR was removed and STRUCTURE analysis/best K search re-performed. Locus removal and data analysis were repeated until only the locus with the highest AR remained. The same procedure was applied on the dataset using PAR, resulting in each 10 reduced datasets for AR and PAR.

STRUCTURE analysis with 11 loci resulted in a best K of 2 (see Results and Discussion). From the 10 repetitions of this K, the individual cluster assignment of the run with the highest LnP(D) was used as benchmark for comparison with the STRUCTURE assignments of the reduced datasets. In detail, at K = 2, STRUCTURE provides for each individual a probability p to belong to cluster A and a probability p = 1 – p to belong to cluster B. We considered an individual assigned to cluster A when its p was larger than 0.5 plus a variable exclusion threshold value x. Accordingly, an individual was assigned to cluster B when p > 0.5 + x. Individuals with p or p larger than 0.5 but smaller than or equal to 0.5 + x were considered as unassignable (U). Thresholds from x = 0.00 to x = 0.50 were used in 0.05-steps. As this algorithm for cluster assignment handles just K = 2, which is also the true number of clusters as discernible from the analysis using all data, also in the two reduced datasets where another value than two was suggested as best K (see Results), the STRUCTURE results for K = 2 were used for computation. On the individual level, this algorithm allowed three outcomes: (i) Correct assignment: An individual was assigned to the same cluster it had been assigned to when using all 11 loci. (ii) Incorrect assignment: An individual was assigned to the cluster opposite to the one it had been assigned to when using 11 loci. (iii) Exclusion: Due to an admixture value close to 0.5, no assignment was achieved. At x = 0.00, no exclusion occurred, while at x = 0.50 all individuals were excluded. Higher levels of x can be expected to reduce cases of incorrect assignment but also reduce the number of individuals assigned at all.

The assignment of each individual to A, B, and U in all loci-reduced datasets was compared with its assignment when using all data, and incorrect assignment rate (rate of individuals assigned to different clusters in benchmark and reduced dataset) and exclusion rate (rate of individuals unassignable in loci-reduced datasets) were calculated in Excel 2013.

Results and discussion

Scoreable alleles were found in 89.6% (DN48) to 98.2% (DN40) of the individual amplicons (Table S2). The number of alleles per locus ranged from 12 to 25, with a mean of 19.36. AR per locus ranged from 7.63 (DN40 in K5) to 20.00 (DN39 in K0), PAR from 0.05 (DN40 in K5) to 3.95 (DN45 in P0). Over all loci, average AR was 13.38 for K0, 13.30 for P0, 11.77 for K5, and 11.18 for P5; average PAR was 0.99 for K0, 1.28 for P0, 0.45 for K5, and 0.51 for P5 (Fig. 1, Table S3). The loss of allelic diversity from Generation 0 to Generation 5 was significant for population P concerning PAR. This loss of diversity can likely be explained by drift acting on a relatively small laboratory population. Thus, the fly lines here are a suitable example for genetic pauperization, a situation that should be avoided in, for instance, laboratory culturing.

Figure 1.

Allelic and Private Allelic Richness of populations Kaserstattalm (K) and Pfitscherjoch (P) in Generation 0 and Generation 5. Pairwise comparisons by two-sided t-tests are indicated by squared brackets. Significant differences (α = 0.05) are indicated by an asterisk.

Significant deviations from HWE after Bonferroni-Holm correction were found in four and one loci in K0 and P0, respectively, and in seven and eight loci in K5 and P5, respectively (Table S4). Thus, the number of loci significantly deviating from HWE had increased from K0 to K5 and from P0 to P5. A similar trend emerged in LD, where the number of significantly linked locus pairs increased from K0 to K5 (3 to 30) and from P0 to P5 (0 to 21) (Table S5). Significant deviations from HWE and significant LD in wild populations like K0 and P0, which were collected in the field to establish the laboratory stocks, would indicate that non-neutral evolutionary forces like inbreeding, non-random mating, and/or selection are acting on the populations. However, both K0 and P0 samples were comparatively small, and more data are needed for reliable population genetic estimates of these wild populations of D. nigrosparsa. More importantly here, already after five generations in the laboratory, clear changes in HWE, LD, and diversity metrics were visible.

The AMOVA revealed significant but very low FST values ranging from 0.01 to 0.03 among the various combinations of populations and generations except between the two Generation 0 populations, for which the values were not significant (Table 1). The FIS value ranged from 0.11 to 0.19 and was always significant in the various combinations. The percentage of variation between populations increased from 0.01% between K0 and P0 to 2.71% between K5 and P5. These results indicate a weak population differentiation, slightly increasing during five generations of laboratory culturing, and are in line with the loss of Allelic Richness observed between Generation 0 and 5.

Table 1.

Results from Analyses of Molecular Variance and F-Statistics from various combinations of populations. df … degrees of freedom; reg … regions; pop … populations; ind … individuals; FST … fixation index of subpopulation compared with total population; FIT … inbreeding coefficient of individuals relative to total population; FIS … inbreeding coefficient of individuals relative to subpopulation; K0, K5 … populations from Kaserstattalm at Generation 0 and Generation 5; P0, P5 … populations from Pfitscherjoch at Generation 0 and Generation 5. All … all hierarchy levels (populations and generations).

AMOVA				Fixation indices
Sample	Source of variation	df	% Variation	Index	Value	p
All	Among reg	1	1.37	FST	0.02	0.01
	Among pop	2	1.07	FIS	0.15	0.01
	Among ind	554	14.82	FIT	0.17	0.01
	Within ind	558	82.74
	Total	1115	100.00
K0 & K5	Among pop	1	0.87	FST	0.01	0.01
	Among ind	277	18.60	FIS	0.19	0.01
	Within ind	279	80.53	FIT	0.19	0.01
	Total	557	100.00
P0 & P5	Among pop	1	1.27	FST	0.01	0.01
	Among ind	277	11.38	FIS	0.12	0.01
	Within ind	279	87.35	FIT	0.13	0.01
	Total	557	100.00
K0 & P0	Among pop	1	0.01	FST	0.00	0.48
	Among ind	60	11.19	FIS	0.11	0.01
	Within ind	62	88.80	FIT	0.11	0.01
	Total	123	100.00
K5 & P5	Among pop	1	2.71	FST	0.03	0.01
	Among ind	494	14.56	FIS	0.15	0.01
	Within ind	496	82.73	FIT	0.17	0.01
	Total	991	100.00

STRUCTURE analysis of K0, P0, K5, and P5 resulted in a best K = 2. In Generation 0, single individuals assigned mainly to one of the two clusters, or intermediate, existed in both populations K and P. After five generations, some sorting had occurred, with one cluster dominated by individuals from K and the other by individuals from P (Fig. 2). This sorting was likely driven by genetic drift, which leads to loss of some (as seen in the analysis of AR) but not the same alleles in the two populations, as seen in the STRUCTURE plot. This sorting is far from being complete; still there are several ‘K-type’ individuals in P5 and vice versa, and many intermediates. By systematically reducing the number of microsatellite loci, we then tested whether an individual was still assigned to the same cluster as with the full set of 11 loci.

Figure 2.

STRUCTURE results of the populations Kaserstattalm (K) and Pfitscherjoch (P) in Generation 0 and Generation 5. K = 2 was the most probable number of clusters. Loss of different alleles due to genetic drift during five generations laboratory rearing lead to increasing population differentiation.

The sequence of locus removal is given in Table 2, and the resulting STRUCTURE plots for K = 2 are given in Fig. 3. The method of Evanno et al. suggested K = 2 for all reduced datasets except for those where just one locus remained. The detailed results of assignment accuracy by stepwise exclusion of microsatellite loci are given in Fig. 4 and Supplemental Tables S6-S9. As expected, the number of correctly assigned individuals decreased both with fewer loci and higher thresholds. A trade-off between accepting incorrect assignments versus accepting high numbers of excluded individuals is inherent in this approach. Considering the often-used error level of 5%, a correct assignment was still possible with eight loci both when sorted by AR and PAR. With fewer than eight loci, also the rate of incorrect assignments exceeded 5% when low thresholds were applied. Remarkably, a correct assignment of more than 80% of the individuals was still possible with as few as three loci. Correct assignment dropped to 76.6% and 74.6% for AR and PAR, respectively, when data of only two loci were used, although these were the most polymorphic markers. Any population structure vanished when just the single locus with the highest AR and PAR was used (Fig. 3). The differences between sorting the loci by AR and PAR were negligible.

Table 2.

Sequence of locus removal for microsatellite data on population Kaserstattalm in Generation 5. Loci were sorted from lowest to highest Allelic Richness (AR) and Private Allelic Richness (PAR) and sequentially removed. The number of alleles remaining and the value suggested for the best K using the method of Evanno et al. (2005) are given. In AR, locus DN16 (AR = 15.96) remained as last marker, in PAR, locus DN41 (PAR = 1.38).

Removal by AR value					Removal by PAR value
n loci	removed locus	AR of removed locus	n alleles in dataset	best K	n loci	removed locus	PAR of removed locus	n alleles in dataset	best K
11	none	—	193	2	11	none	—	193	2
10	DN40	7.73	174	2	10	DN40	0.19	174	2
9	DN45	9.13	159	2	9	DN31	0.27	158	2
8	DN35	9.83	135	2	8	DN45	0.32	143	2
7	DN49	10.24	121	2	7	DN48	0.33	122	2
6	DN48	10.64	100	2	6	DN36	0.59	110	2
5	DN31	11.19	84	2	5	DN39	0.59	94	2
4	DN41	12.12	69	2	4	DN35	0.74	70	2
3	DN39	12.90	53	2	3	DN16	0.75	47	2
2	DN36	12.94	41	2	2	DN49	0.92	33	2
1	DN37	13.54	23	7	1	DN37	1.25	15	4

Figure 3.

STRUCTURE results after sequential removal of loci, sorted according to Allelic (AR) and Private Allelic Richness (PAR). Plots in which the estimate for the best K deviated from 2 are printed in grey. While the separation quality gradually deteriorated, the major population structure was still visible when only the two loci with highest AR or PAR are used.

Figure 4.

Based on the STRUCTURE results, individuals genotyped at one to ten loci were assigned to two clusters. Individuals with intermediate probabilities to belong to one cluster were excluded using a variable threshold x; x = 0 results in no exclusions, x = 0.5 excludes all individuals. The cluster assignment with 11 loci was used as a benchmark. Each individual could be either correctly assigned (i. e., in the same cluster as with 11 loci; upper row of plots), incorrectly assigned (i.e., in the other cluster than with 11 loci; middle row), or excluded from assignment based on the threshold value (lower row). High thresholds minimize incorrect assignments at the cost of many excluded individuals.

Screening for loss of genetic diversity is of high importance for laboratory culturing and will often require the screening of many individuals. While the costs of genotyping a single microsatellite locus are around 3 Euro, costs pile up with adding loci, especially when the population to be screened is large. In this study, we have provided the proof of concept for a method of reducing the number of loci while still detecting genetic structure and its changes across generations. This method is especially promising when core sets, that is, smaller sets of markers selected from a large set of microsatellite markers, are to be established for broad-scale screening. The particular number of loci one removes will depend on whether maximum accuracy or information on a maximum number of individuals is desired. Using AR and PAR as a criterion for the removal of loci yields similar results. In adopting this method, long-time and large-scale population genetic investigations will benefit from a substantial reduction of costs in terms of money and processing time.