Effects of Recurrent Selection on Population Structure and Allele Frequencies in the M3S Maize Population

Abstract

The effects of four cycles of recurrent selection on the allele frequencies of simple sequence repeat (SSR) markers and population structure were examined in the Maksimir 3 Synthetic (M3S) maize population (Zea mays L.). Genotyping of 32 plants from each selection cycle at 38 SSR loci revealed that the mean number of alleles per locus and the mean expected heterozygosity were preserved across selection cycles, indicating the maintenance of sufficient genetic variability in the population required for future genetic gain. The Waples test of selective neutrality revealed that genetic drift was the main force in changing allele frequencies in the population. The proportion of selectively non-neutral loci in single cycles of selection varied between 16% and 37%. Some non-neutral loci shared the same genomic locations with previously published QTLs controlling important agronomic traits. An analysis of molecular variance revealed that 5.6% of the genetic variation occurred among and 94.4% within cycle populations. Between 5% and 29% of loci were found to be in a significant Hardy–Weinberg (HW) disequilibrium, with the majority showing an excess of homozygosity. The excess of homozygosity at several loci was highly consistent across cycle populations, suggesting positive assortative mating as a possible cause of the observed HW disequilibrium. Linkage disequilibrium (LD) tests revealed that the M3S population was essentially in linkage equilibrium. The proportion of pairs of loci in significant LD varied from 0.1% to 1.8% across selection cycles, probably due to the effects of genetic drift and epistatic selection.

1. Introduction

A broad class of selection methods referred to as recurrent selection (RS) uses a cyclical approach to gradually increase the frequency of favorable alleles affecting quantitatively inherited traits in broad-based plant populations while maintaining genetic variability for future selection [1]. The two primary forces affecting allele frequencies in RS programs are selection, which increases the frequency of favorable alleles, and genetic drift, which is a random change in allele frequency due to small population sizes [2]. An expected effect of random genetic drift in a population is the dispersion of allele frequencies from intermediate values toward the extremes [3], a phenomenon that has been observed in maize populations under RS [4,5,6,7]. Problems associated with finite population sizes, such as the fixation of undesirable alleles due to random genetic drift, affect the response to selection [8], implying a trade-off between short-term selection favored by higher intensities and long-term selection favored by higher effective population sizes [9]. A relevant question is the extent to which selection influences allele frequencies in a population subjected to recurrent selection. A number of studies have reported significant changes in allele frequencies at isozyme, RFLP, SSR, and SNP marker loci due to selection [2,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Falke et al. [2], Coque and Gallais [12], and Wisser et al. [13] identified marker loci, whose allele frequencies were significantly changed by selection, to reside in genomic regions where QTLs affecting important maize agronomical traits of maize had previously been identified. The observed correlation between allele frequency changes at markers linked to QTLs and phenotypic performance suggested the possibility of using these markers in marker-assisted breeding programs [2,12]. In addition to information on alleles subjected to selection, the extent of genetic diversity in improved selection cycles is also important information for breeders because the future response of the population to selection depends on it. Several studies have reported a significant decrease in genetic diversity as measured by the number of marker alleles per locus and expected heterozygosity (gene diversity) in populations subjected to RS [7,14,18,19,20,21,22]. On the other hand, Kolawole et al. [17] reported that the changes in the different measures of genetic diversity due to selection in two maize composites were either small or negligible.

RS methods were developed to gradually improve the mean performance of genetically broad-based populations, and their incorporation in current breeding programs could facilitate the use of exotic germplasm [23] as well as locally adapted but non-improved germplasm [24] to increase the genetic base of maize breeding germplasm. Several researchers around the world, aware of the danger posed by the loss of genetic variability associated with the abandonment of landraces, made efforts to collect them before their total disappearance [25]. With the aim to exploit genetic variability existing in locally adapted maize germplasm from Southeast Europe, a synthetic maize population named Maksimir 3 Synthetic (M3S) was developed at the Faculty of Agriculture, University of Zagreb (Croatia) by intercrossing 12 inbred lines originating from locally adapted open-pollinated varieties and landraces from different regions of former Yugoslavia [26]. The 12 M3S progenitor lines showed equally good combining ability for grain yield with both the BSSS and Lancaster testers [27] and the same level of isoenzymic differences to both testers [28]. After its creation, M3S was subjected to four cycles of intrapopulation RS primarily for grain yield, but the resistance to leaf and stalk diseases was also considered [26,29,30,31]. After two selection cycles, grain yield increased slightly in the population per se, but the selection was more effective in reducing inbreeding depression for grain yield [26]. After the third cycle of selection, Sabljo et al. [29] observed no further improvement in grain yield, but stalk rot incidence significantly decreased. In the fourth cycle of selection, performed under low and high N conditions (0 and 150 kg N ha⁻¹, respectively), Bukan et al. [30] found an indication of the specific adaptation of the two C4 populations to contrasting N environments. The resulting C4N0 population, which has been developed by intercrossing individuals selected in a low N environment, showed a significant reduction in ASI, which is typical for the genetic material more suited to low N environments. Although no significant yield increase was observed from C3 to C4, both resulting C4 populations (C4N0 and C4N150) performed well in their target N fertilization environments. Also, after applying RS methods, the resistance to stalk rot diseases in the M3S population appeared to be maintained [31].

Šarčević et al. [6] examined changes in allele frequencies at nine SSR loci in the M3S maize population after two cycles of RS and found significant changes in allele frequencies at four loci. In some previous studies, including Šarčević et al. [6], genetic changes in populations under RS were assessed by comparing the starting (base) populations and the advanced populations developed through a certain number of selection cycles. However, according to Coque and Gallais [12], the main problem to solve when studying more than one selection cycle is to test whether allelic frequency changes are due to selection or genetic drift. To obtain a better insight into the genetic changes of a population under RS, it would be useful to monitor the effects of selection at the molecular level from cycle to cycle. In the present study, we examined the M3S maize population at the molecular level after four cycles of RS by monitoring the selection response in each single cycle of selection using 38 SSR markers. The objectives of the study were (1) to investigate changes in allele frequencies in the population due to the effects of random genetic drift and selection and (2) to investigate population structure concerning the partitioning of genetic variation among and within cycle populations and concerning the Hardy–Weinberg (HW) disequilibrium at individual loci and linkage disequilibrium (LD) between pairs of loci within cycle populations.

2. Materials and Methods

2.1. Development of M3S Cycle Populations

The Maksimir 3 Synthetic (M3S) cycle 0 population (C0) was developed by intermating 12 maize inbred lines, whose origin traces back to several landraces and open-pollinated varieties from different regions of former Yugoslavia [26]. Four cycles of intrapopulation recurrent selection (RS), primarily for grain yield, were conducted in the M3S (Figure 1, Table S1).

In the first cycle, two-stage S1-S2 selection was performed with emphasis on disease resistance (Setosphaeria turcica and Colletotrichum graminicola (Ces) GW Wils) among S1 progenies and grain yield among S2 progenies resulting in the C1 population. In the second cycle, selection for higher grain yield was conducted among S1 progenies, resulting in the C2 population. In the third cycle, two selection methods, including evaluation of S1 and full-sib (FS) progenies, were conducted simultaneously to improve grain yield, resulting in two C3 populations, C3S1 and C3FS, respectively. The low response to selection observed in the C3S1 population [29] led to the fourth cycle of selection starting from the C3FS population, in which S1 progenies were evaluated simultaneously at low nitrogen fertilization (0 kg N ha⁻¹) and high nitrogen fertilization (150 kg N ha⁻¹), resulting in two C4 populations, C4N0 and C4N150, respectively. Details of the experimental procedures used in population synthesis and the four selection cycles have been described previously [26,29,30,31].

2.2. Simple Sequence Repeat (SSR) Genotyping

Thirty-two plants randomly selected from each of the C0, C1, C2, C3S1, C3FS, C4N0, and C4N150 cycle populations of M3S were grown in a growth chamber. After about three weeks, DNA was extracted from each plant using a GenElute Plant Genomic DNA Miniprep Kit (Sigma-Aldrich, St. Louis, MO, USA). Forty SSR primer pairs were chosen for analysis on the basis of previous polymorphism identification in M3S [6] and based on genomic location in order to provide a uniform coverage of all ten maize chromosomes. After an initial analysis, two markers were discarded because of poor amplification. Primer pairs were fluorescently labeled prior to polymerase chain reactions (PCRs). PCRs were performed in 15 μL final volumes containing 25 ng of template DNA, a 1× PCR buffer with added 1.5 mM MgCl₂, 0.2 μM of each of the forward and reverse primers, 0.2 mM of dNTPs, and 0.5 U of Taq polymerase (Sigma-Aldrich, St. Louis, MO, USA). Reactions were carried out on a Veriti 96-Well Fast Thermal Cycler (Applied Biosystems, Foster City, CA, USA). Reactions were denatured at 95 °C for 2 min, followed by 30 cycles of 92 °C for 1 min, 55 °C for 1 min, and 72 °C for 2 min. This was followed by a final extension step at 72 °C for 10 min. Diluted PCR products (in ddH20) were mixed with 8.75 μL HiDi (Applied Biosystems, Foster City, CA, USA) and 0.25 μL Genescan 500 LIZ size standard (Applied Biosystems, Foster City, CA, USA), centrifuged, denatured at 95 °C for 5 min, and placed on ice. Fragment analysis was performed using an ABI 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). GeneMapper 4.0 software (Applied Biosystems, Foster City, CA, USA) with the sizing algorithm “2nd order least square” was used to detect allele sizes.

2.3. Statistical Analysis of the SSR Data

2.3.1. Diversity Statistics and Allele Frequency

For each cycle population, the total number of SSR alleles (At), the mean number of alleles per locus (Am), the observed heterozygosity (Ho), and the expected heterozygosity (He) were calculated. Allele frequencies were estimated using the software package Genealex 6 [32]. The Waples test of temporal variation in allele frequency [33] was used to test the selective neutrality of alleles at 38 SSR loci from cycle to cycle, as well as after four cycles of selection. The Waples method tests the hypothesis that observed differences in allele frequency can be explained entirely by sampling processes, both in choosing gametes to form the next generations (genetic drift) and in choosing the sample for genetic analysis (sampling error). Loci with one or more non-neutral alleles (rejecting the null hypothesis of the Waples test) were designated as non-neutral loci. The effective population size (Ne) for the particular cycle of selection was assumed to be equal to the number of intermated progenies (N) in that cycle. Because S2 progenies were intermated in the first cycle of selection, Ne for the first cycle of selection was corrected by multiplying N with the term 1/(1 + Fp), where Fp = 0.5 is the inbreeding coefficient of the parental generation (S1) from which the S2 progenies were produced [34]. The Ne after four selection cycles was calculated as the harmonic mean of the Ne values from individual selection cycles [3].

2.3.2. Population Structure

Wright’s fixation index was estimated using the formula F_IS = 1 − (Ho/He) to quantify the lack or excess of heterozygosity. Significant deviations from the Hardy–Weinberg (HW) equilibrium at individual loci were determined using the likelihood ratio G test [35]. The significance of linkage disequilibrium (LD) between pairs of alleles from different loci was determined using the χ² test outlined by Weir [36]. The calculations of Ho, He, and F_IS, as well as the HW equilibrium and LD tests, were performed using the software package Popgene 1.31 [37].

To compare the observed proportions of pairs of loci in a significant LD involving linked vs. unlinked loci, as well as non-neutral vs. neutral loci (based on the Waples neutrality test) in a particular cycle population, we derived two equations (Equations (1) and (2)). The probability that the two randomly selected loci are linked $\left(P_{ll}\right)$ was calculated as the number of possible pairs involving linked loci relative to the number of possible pairs involving all loci (Equation (1)):

$$P_{ll}=\frac{{{\displaystyle \sum }}_{i=1}^{10} n_i\left(n_i-1\right)}{N\left(N-1\right)}$$

(1)

where $n_i$ is the number of loci at the i-th chromosome, and $N$ is the total number of loci. The probability that the two randomly selected loci are non-neutral $P_{nnl}$ was calculated as the ratio of the number of possible pairs of non-neutral loci and the number of possible pairs of all loci (Equation (2)):

$$P_{nnl}=\frac{n_{nnl} \left(n_{nnl}-1\right)}{N\left(N-1\right)}$$

(2)

where $n_{nnl}$ is the number of non-neutral loci, and N is the total number of loci. The derivation of Equations (1) and (2) is shown in Appendix A.

In order to detect differences in the distribution of multilocus genotypes among and within the cycles of selection, we performed an analysis of molecular variance (AMOVA) [38] using Arlequin 3.5 [39]. The sources of variation included seven cycle populations (C0, C1, C2, C3S1, C3FS, C4N0, and C4N150) and 224 plants representing the entire sample. The pairwise fixation index, Φ_ST, provides an indication of the genetic distance between cycle populations [12]. The significance of the average fixation index, as well as of each pairwise Φ_ST value, was obtained after 1023 permutations.

3. Results

3.1. Diversity Statistics and Allele Frequency

The 38 SSR primer pairs generated a total of 133 different alleles (an average of 3.5 alleles per locus) across seven M3S cycle populations. The number of alleles per locus varied from two to seven. The total number of alleles (At) for individual cycle populations varied from 123 (mean of 3.2 alleles per locus) found for C3FS to 130 (mean of 3.4 alleles per locus) found for C1 and C2, although the differences were not significant (

Table 1. Genetic diversity of the seven Maksimir 3 Synthetic (M3S) cycle populations: total number of alleles (At), mean number of alleles per locus (Am), and mean expected heterozygosity (He) with their standard errors.

Cycle Population	At	Am	He
C0	127	3.34 ± 0.242	0.5170 ± 0.031
C1	130	3.42 ± 0.231	0.5318 ± 0.025
C2	130	3.42 ± 0.234	0.5171 ± 0.029
C3S1	126	3.32 ± 0.207	0.5135 ± 0.030
C3FS	123	3.24 ± 0.218	0.5012 ± 0.030
C4N0	129	3.39 ± 0.240	0.5266 ± 0.024
C4N150	125	3.29 ± 0.237	0.4924 ± 0.028

). Mean expected heterozygosity (He) ranged from 0.49 in C4N150 to 0.53 in C1 but observed differences among cycle populations were not found to be significant. Six alleles found in advanced cycle populations were not found in the original base population (C0). Seventeen alleles from 13 loci that were present in C0 were absent from one or more improved cycle populations. Most of these alleles in the base population were found in low frequencies (less than 0.10).

Mean allele frequencies remained unchanged (0.29 ± 0.02) after four cycles of selection, but the shape of the frequency distribution changed slightly (

Table 2. Allele frequency distribution (%) in seven Maksimir 3 Synthetic (M3S) cycle populations.

Allele Frequency Class	Cycle Population
	C0	C1	C2	C3S1	C3FS	C4N0	C4N150
	%
0.00–0.10	23	28	28	26	29	27	29
0.11–0.20	25	12	21	22	20	15	20
0.21–0.30	13	21	15	14	9	17	12
0.31–0.40	11	11	8	10	13	14	12
0.41–0.50	8	10	8	11	10	8	4
0.51–0.60	8	6	8	6	6	4	7
0.61–0.70	5	2	4	3	6	5	7
0.71–0.80	4	8	5	3	2	8	4
0.81–1.00	3	2	3	5	5	2	5
Mean allele frequency ± SE	0.29 ± 0.021	0.29 ± 0.020	0.29 ± 0.021	0.29 ± 0.021	0.29 ± 0.021	0.29 ± 0.020	0.29 ± 0.022

). Generally, slightly higher percentages of alleles with low (≤0.10) and high (>0.80) frequencies were observed in improved cycle populations relative to the base population.

The number of non-neutral loci detected by the Waples test varied among single cycles of selection from six (as found between C2 and C3S1) to 14 (as found between C2 and C3FS) (

Table 3. Chromosomal position of non-neutral loci in single cycles and across four cycles of selection in the Maksimir 3 Synthetic (M3S) population.

Bin	Locus	SCP 1	C0	C1	C2	C2	C3FS	C3FS	C0	C0
		RCP 2	C1	C2	C3S1	C3FS	C4N0	C4N150	C4N0	C4N150
1.01	phi056
1.03	phi339017
1.08	dupssr12
1.11	phi064
2.00	phi96100
2.02	phi098
2.04	phi083
3.02	phi036
3.05	phi053
3.07	bnlg197
4.01	phi213984
4.04	phi308090
4.05	phi438301
4.11	phi076
5.00	nc130
5.01	bnlg143
5.03	phi109188
5.07	phi128
6.00	phi126
6.03	phi389203
6.04	phi452693
6.07	phi123
7.02	phi034
7.03	bnlg572
7.05	phi082
7.06	phi116
8.00	umc1359
8.03	phi115
8.08	phi015
8.09	phi233376
9.01	phi033
9.02	umc1033
9.05	phi236654
9.07	bnlg128
10.02	phi96342
10.03	phi050
10.04	phi084
10.3	phi059

¹ SCP: starting cycle population. ² RCP: resulting cycle population. Significant Waples neutrality test at the 0.05 probability level.

). In the fourth cycle of selection, more non-neutral loci were found for C4N150 than for C4N0 (nine vs. seven). However, the cumulative changes in allele frequency over four selection cycles resulted in a higher number of non-neutral loci for C4N0 (six) than for C4N150 (three). Considering single cycles of selection, non-neutral loci were found on all ten maize chromosomes, with the lowest incidence on chromosome 6 (four cases) and the highest incidence on chromosome 1 (nine cases). The occurrence of non-neutral loci was inconsistent across cycles of selection, but discrepancies were also observed between the neutrality status of loci after four cycles of selection and their neutrality status across individual cycles of selection. Of the 38 loci, 28 were non-neutral in at least one single cycle, 5 were non-neutral across three cycles, and none of them were non-neutral across all four cycles of selection. Even in cases where a particular marker locus was recognized as non-neutral over several cycles of selection, there was inconsistency in the neutrality status of different alleles at these loci (Table S2).

3.2. Analysis of Molecular Variance (AMOVA)

Analysis of molecular variance revealed that 5.6% of the genetic variation occurred among and 94.4% within M3S cycle populations. The average fixation index, Φ_ST, commonly used to estimate the extent of differentiation in population subdivisions [19], was 0.0559. The pairwise Φ_ST values between cycle populations, interpreted here as the genetic distance between them, were all significant (

Table 4. Pairwise Φ_ST values between Maksimir 3 Synthetic (M3S) cycle populations (bellow diagonal) and their probability values after 1023 permutations (above diagonal).

Cycle Population	C0	C1	C2	C3S1	C3FS	C4N0	C4N150
C0		<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
C1	0.0503		<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
C2	0.0395	0.0535		<0.0001	<0.0001	<0.0001	<0.0001
C3S1	0.0437	0.0575	0.0143		<0.0001	<0.0001	<0.0001
C3FS	0.0702	0.0872	0.0542	0.0750		<0.0001	<0.0001
C4N0	0.0598	0.0664	0.0537	0.0728	0.0247		<0.0001
C4N150	0.0705	0.0843	0.0595	0.0739	0.0320	0.0143
Mean 1	0.0560	0.0670	0.0460	0.0560	0.057	0.0490	0.0560

¹ Mean Φst value of a cycle population to the remaining cycle populations.

).

The highest Φ_ST was observed between C1 and C3FS (0.0872), followed by Φ_ST observed between C1 and C4N150 (0.0843). The lowest Φ_ST values were determined between C2 and C3S1 and between C4N0 and C4N150 (0.0143). The highest mean Φ_ST value of a cycle population to the remaining cycle populations was observed for C1.

3.3. Hardy–Weinberg (HW) Disequilibrium

The number of loci that were not in the HW equilibrium varied across M3S cycle populations from 2 (5%), observed in C3FS, to 11 (29%), observed in C3S1 (Figure 2). Most of them showed an excess of homozygotes.

In total, 26 SSR loci (68%) were identified to be in significant HW disequilibrium in at least one cycle population, and such loci were found on all chromosomes except chromosome 3 (Table S3). Fixation indices (F_IS) for six loci, which were not in HW equilibrium in three or more of the cycle populations, are shown in

Table 5. Fixation index for six loci across seven Maksimir 3 Synthetic (M3S) cycle populations. Only those loci with significant departure from the HW equilibrium in three or more cycle populations are shown.

Locus	Bin	Cycle Population
		C0	C1	C2	C3S1	C3FS	C4N0	C4N150	Mean
		Fixation index (FIS)
phi098	2.02	−0.12	−0.38 *	ml 1	1.00 **	−0.21	−0.27	0.53 *	0.09
phi076	4.11	−0.07	0.43 **	0.93 **	−0.07	−0.00	1.00 **	−0.02	0.31
bnlg143	5.01	0.09	0.49 **	0.21	−0.03	0.19 *	0.23 *	0.26 *	0.21
phi452693	6.04	0.34 *	0.50 **	−0.21	−0.03	−0.19	0.30 *	0.15	0.13
bnlg572	7.03	0.30 **	0.86 **	0.23	0.91 **	0.59 **	0.90 **	0.77 **	0.65
phi050	10.03	0.28 **	0.20	0.57 **	0.10 **	−0.02	0.19	1.00 **	0.33

*, ** Significant at the 0.05 and 0.01 probability levels, respectively. ¹ ml monomorphic locus.

. Two of these loci (bnlg143 and phi050) showed a significant excess of homozygosity (positive F_IS value) in five and one locus (bnlg572) in six of seven cycle populations. All six loci had a positive F_IS, averaged over populations, with the highest value observed for the locus bnlg572 (0.65).

3.4. Linkage Disequilibrium (LD)

In the base population (C0), 3 (0.4%) out of 703 possible pairs of loci were in significant LD (

Table 6. Pairs of loci in significant (p < 0.05) linkage disequilibrium (LD) in seven Maksimir 3 Synthetic (M3S) cycle populations (C0, C1, C2, C3S1, C3FS, C4N0, and C4N150).

Bin	Locus	Bin	Locus	Bin	Locus	Bin	Locus
C0				C3S1
3.07	bnlg197	4.04	phi308090	9.01	phi033	9.02	umc1033
3.07	bnlg197	8.00	umc1359	C3FS
8.08	phi015	8.09	phi233376	9.07	bnlg128	9.01	phi033
C1				3.07	bnlg197	1.01	phi056
9.07	bnlg128 1	9.01	phi033	3.07	bnlg197	6.00	phi126
9.07	bnlg128	5.03	phi109188	3.07	bnlg197	6.04	phi452693
3.07	bnlg197	1.08	dupssr12	8.08	phi015	8.09	phi233376
3.07	bnlg197	3.02	phi036	1.01	phi056	6.00	phi126
3.07	bnlg197	4.11	phi076	7.05	phi082	7.06	phi116
3.07	bnlg197	2.00	phi96100	C4N0
1.08	dupssr12	8.08	phi015	1.08	dupssr12	8.08	phi015
1.08	dupssr12	4.11	phi076	1.08	dupssr12	9.02	umc1033
1.08	dupssr12	4.05	phi438301	8.08	phi015	8.00	umc1359
1.08	dupssr12	2.00	phi96100	1.01	phi056	9.02	umc1033
8.08	phi015	4.05	phi438301	7.05	phi082	7.06	phi116
3.02	phi036	7.06	phi116	C4N150
10.03	phi050	10.04	phi084	3.07	bnlg197	1.08	dupssr12
C2				3.07	bnlg197	9.02	umc1033
3.07	bnlg197	1.08	dupssr12	9.01	phi033	8.03	phi115
3.07	bnlg197	10.03	phi050	4.04	phi308090	4.05	phi438301
3.07	bnlg197	6.03	phi389203
3.02	phi036	7.05	phi082

¹ Selectively non-neutral loci (according to the Waples neutrality test at the 0.05 probability level) are indicated in bold.

). The number of significant LD tests in the improved cycle populations varied from only 1 (0.1%) in C3S1 to 13 (1.8%) in C1. Out of three pairs of loci found to be in significant LD in C0, only one pair (which included two linked loci on chromosome 8) was found to be later in LD (in C3FS). Besides this pair, eight cases of newly generated LD between linked loci were found in the improved cycle populations. However, the majority of the newly generated LD involved unlinked loci. The proportions of marker pairs in LD across cycle populations involving linked loci were 0.33 (C0), 0.23 (C1), 0.00 (C2), 1.00 (C3S1), 0.43 (C3FS), 0.40 (C4N0), and 0.25 (C4N150), which is in most cases more than three times higher than the probability that the two randomly selected loci were linked (0.08, according to Equation (1)).

Some loci were found more frequently to be in significant LD, like bnlg197 and phi015 (in five and four out of seven cycle populations, respectively), but there were only three pairs of loci in significant LD, which were found at the same time in two or more improved cycle populations (dupssr12 and phi015 in C1 and C4N0; bnlg197 and dupssr12 in C1, C2, and C4N150; and bnlg128 and phi033 in C1 and C3FS). Most pairs of loci found to be in significant LD in improved cycle populations included one or both non-neutral loci (according to the Waples neutrality test). The observed proportion of LD pairs, including both non-neutral loci, was at least three times higher than expected (according to Equation (2)) in four out of five improved cycle populations (C3S1 was not considered because only one test was significant) with the observed to expected ratio being 0.31:0.11, 0.25:0.08, 0.43:0.13, 0.00:0.03, and 0.50:0.05 for C1, C2, C3FS, C4N0, and C4N150, respectively.

4. Discussion

4.1. Diversity Statistics and Allele Frequency

In the present study, the effects of four cycles of recurrent selection on the allele frequencies of simple sequence repeat (SSR) markers and population structure of the Maksimir 3 Synthetic (M3S) maize population were examined. The variability of the C0 population, in terms of the mean number of alleles per locus and expected heterozygosity, was in the order of magnitude of the starting population of other RS programs [5,7,9,14,19].

The mean number of alleles per locus and the mean expected heterozygosity did not change significantly in the M3S population after four cycles of RS. Similarly, in the study of Wisser et al. [40], neither of the two diversity measures, determined by SSR markers, changed significantly after four cycles of RS for quantitative disease resistance in a complex maize population from CIMMYT. Kolawole et al. [17] and Wisser et al. [40] reported that the changes in different measures of SNP diversity were either small or negligible in maize populations subjected to recurrent selection. In the study by Daas et al. [41], the genetic diversity of the two maize populations also did not change significantly after two cycles of genomic selection. In contrast, Labate et al. [18] and Hinze et al. [19] observed a significant decrease in the mean number of alleles per locus and expected heterozygosity in the BSSS and BSCB1 populations after 12 and 15 cycles of reciprocal RS, respectively. A decrease in marker diversity (in terms of the number of polymorphic loci, mean number of alleles per locus, expected heterozygosity, or observed heterozygosity) within maize populations that underwent various numbers of RS cycles has also been reported in several previous studies [7,14,20,21,22,42,43,44,45,46,47].

In the present study, some alleles found in the base population (C0) were absent from subsequent cycle populations, while some alleles absent from the base population were detected in one or more subsequent cycle populations. Most of the missing alleles were generally found at low frequencies (less than 0.10), similar to the study of McLean-Rodríguez et al. [48], where most of the alleles lost or gained over time in 13 Mexican landraces had rare or low frequencies. Such rare alleles may not have been detected in the particular cycle population in the present study because of the relatively small sample size (32 plants per cycle population genotyped), which was also the case in the study by McLean-Rodriguez et al. [48], who sampled only 10 plants per population. Nevertheless, our sample size is comparable to the sample size of 30 plants used for the SSR analysis of the BSSS and BSCB1 populations studied by Hinze et al. [19]. In an earlier study, Labate et al. [18] genotyped 100 individuals from the same two maize populations using RFLP markers and reported higher estimates of average number of alleles per locus, expected heterozygosity, heterozygous plants, and number of unique alleles, which reflects not only the differences between the two types of markers used in the two studies (RFLPs versus SSRs), but also the power of larger sample sizes in detecting less frequent alleles [19]. Although the sample size of 32 genotyped plants per population in the present study was relatively small, it is comparable to sample sizes reported in some previous studies in maize using SSR markers [17,18] as well as SNP markers [19,20,48], which ranged from 10 to 36 individuals per population. However, the possible role of pollen or seed contamination during the development of the M3S cycle populations cannot be excluded either.

In improved M3S cycle populations compared to C0, the mean allele frequency did not change, although slightly higher proportions of alleles with low and high frequencies were found. Similar changes in allele distribution after various cycles of RS in maize have been earlier found by Labate et al. [4], Pinto et al. [5], and Šarčević et al. [6], whereas Kolawole et al. [17] reported the opposite, a decrease in the proportion of alleles at both low and high frequencies and an increase in those at intermediate frequencies.

The changes in allele frequencies from cycle to cycle as well as after four cycles of selection in M3S as determined by the Waples test, were mainly attributable to the effects of random genetic drift. Similar results were reported in previous studies examining changes in allele frequencies in maize populations subjected to RS [2,4,5,43]. Assuming the value of Ne = N, the Waples test identified six (16%) and three (8%) non-neutral loci after four cycles of selection from C0 to C4N0 and C4N150, respectively (Table 3). The number of non-neutral loci in single cycles varied between 6 (16%), as found between C2 and C3S1, and 14 (37%), as found between C2 and C3FS. Labate et al. [4] observed 17% non-neutral loci in BSSS(R) and BSCB1(R) populations after 12 cycles of reciprocal RS. Pinto et al. [5] reported significant changes in allele frequency due to selection at 13% and 7% of SSR loci in two tropical maize synthetics subjected to a single cycle of high-intensity reciprocal RS. Falke et al. [2] detected 20.13% non-neutral loci in one and 12.87% in a second biparental maize population after four and seven cycles of intrapopulation RS, respectively. In these studies, loci with significant changes in allele frequency due to selection were not restricted to particular chromosomes or genomic regions but were scattered throughout the genome. In our study, selectively non-neutral loci were also found on all ten maize chromosomes, but their number and chromosomal position varied among cycles of selection. The occurrence of non-neutral loci was inconsistent across the four selection cycles, but discrepancies were also observed between the neutrality status of loci after four cycles of selection and their neutrality status across individual cycles of selection (Table S2). Even in cases where a particular marker locus was recognized as non-neutral over multiple cycles of selection in the M3S population, there were discrepancies in the neutrality status of different alleles at these loci. In addition to selfed progeny RS (used through all four cycles of selection), FS RS was implemented in the third cycle, and, in addition to yield, other traits such as disease resistance in the first and N use efficiency in the fourth selection cycle were also considered. These factors may have contributed to the selection pressure on different QTLs during the four selection cycles in the M3S population. According to Wisser et al. [13], the most important drawback to selection mapping of an individual trait arises if selection is exerted for multiple traits, which is typically the case in breeding populations used for production. In such cases, selection mapping cannot distinguish the loci that respond to a particular selection pressure.

It has also been shown that the effect of QTLs on trait values can vary in different environments [49,50], leading to significant QTL × environment interaction (QEI). Because the selection of progenies for the recombination in the different cycles of RS in the M3S population was based on data collected in different environments, QEI may have influenced the inconsistency of neutrality test results between individual cycles of selection. The reason for this, besides QEI, could be the fact that more than two alleles (up to seven) were found in the population at these marker loci. Thus, it can be assumed that more than one marker allele per locus was initially (in the base population) linked to favorable as well as unfavorable alleles at a particular QTL, leading to random changes in frequencies within the two groups of marker alleles as a result of selection pressure at that QTL.

The C1 cycle population of M3S was created by intermating the highest-yielding S2 progenies after stringent selection for disease resistance among the preceding S1 progenies. The applied two-stage selection method resulted in a population with the highest mean Φ_ST value (Table 4) between a single cycle population and all other cycle populations (mean Φ_ST = 0.067). Selection for two generally negatively correlated traits might increase the genetic distance of C1 to other cycle populations developed through selection for yield only. The observed differentiation of C1 based upon molecular data was also observed on the phenotypic level reported by Bukan et al. [30] (decreased yield of C1 at both N fertilization levels investigated). Yield decreases after primary selection for pest resistance were also reported by Devey and Russell [51] and Klenke et al. [52]. Butrón et al. [43] found a significant linear trend for the departure from the random genetic drift model for some allelic versions of the two SSR markers, umc1329, and phi076, in their study of molecular changes in the maize composite during the selection for resistance to pink stem borer. In the C1 population of M3S, a significant non-neutral SSR marker was also phi076. In the third cycle of selection, a large difference in the number of non-neutral markers was observed between the S1 and FS methods of selection (14 vs. 6 from C2 to C3FS and from C2 to C3S1, respectively). The higher number of non-neutral markers found for the FS method is in accordance with the higher yield and disease resistance observed for C3FS compared to C3S1 [29]. The pairwise Φ_ST values between C3S1 and C3FS (Table 4) showed a divergence of the two populations from each other, confirming the different effects of the two methods of selection applied in the third cycle. In the fourth cycle of selection, we observed a higher number of non-neutral SSR loci in C4N150 than in C4N0 (nine vs. seven from C3FS to C4N150 and from C3FS to C4N0, respectively). Coque and Gallais [12] also found more SSR loci under selection in high N fertility environments. The same authors found that the two genomic regions responding to selection were common to both high N and low N conditions, which, according to them, corroborates the observation of Bertin and Gallais [53] that grain yield QTLs detected in low N conditions were very often a subset of QTLs detected in high N conditions, but probably differentially expressed. Three SSR markers used by Coque and Gallais [12] were located in genomic regions found to be associated with grain yield, N uptake, and kernel number under both high and low N conditions (bnlg1643); grain yield and kernel weight under low N conditions (umc1653); and N utilization efficiency under both high and low N conditions (bnlg1402). These three SSR markers share the same bin location (1.08, 6.07, and 9.02, respectively) as the three selectively non-neutral SSR markers (dupssr12, phi123, and umc1033) in the fourth cycle of selection of the present study, which was conducted under contrasting N fertilization regimes. The C4N0 cycle population, besides exhibiting possible adaptation to low N conditions, also exhibited a significant reduction of anthesis–silking interval (ASI) compared to earlier cycle populations [30]. Two SSR loci that were selectively non-neutral in the C4N0 population (dupssr12 and phi438301) had the same bin location (1.08 and 4.05, respectively) as the RFLP and SSR markers previously found to be associated with QTL affecting ASI in diverse sets of environments [54]. Recent studies [55,56] also reported that significant SNP bases and QTLs for ASI delay due to drought or high-density stress were located on chromosomes 1 and 4.

4.2. Hardy–Weinberg (HW) Disequilibrium

In the present study, the loci that were in a significant HW disequilibrium, in most cases, exhibited an excess of homozygosity, which is consistent with the results reported in previous studies [9,14,19,20,43,57]. Factors like sample size used during random mating and sample size used to estimate the HW equilibrium [57] or genotyping errors [19] might affect departure from HW equilibrium in the M3S cycle populations in both directions, i.e., towards an excess of homozygosity as well as heterozygosity. On the other hand, positive assortative mating within a population due to genotypic differences in flowering time is expected to increase homozygosity in the population [3]. The excess of homozygosity at several loci (bnlg143, bnlg572, and phi050) was highly consistent across M3S cycle populations, suggesting a possible role of positive assortative mating in the observed departure from HW equilibrium at these loci. One of these loci (bnlg572), which showed an excess of homozygosity (positive F_IS value) in seven out of seven M3S cycle populations, shared the same bin location (7.03) as the SSR locus phi114, which had a positive fixation index in nine out of nine maize populations studied by Ordas et al. [9]. Similarly, the SSR loci bnlg572 and phi050, which showed the excess of homozygosity over cycles of selection in the M3S population, had the same bin locations (7.03 and 10.03, respectively) as SNP markers previously found to be associated with QTLs for flowering time-related traits [55].

4.3. Linkage Disequilibrium (LD)

The LD test revealed that the M3S population was essentially in linkage equilibrium across selection cycles, with the number of significant LD tests varying from only one (0.14%) in C3S1 to 13 (1.85%) in C1. For the three pairs of loci found to be in LD in the base population (C0), we assumed that they originated either from parental LD or that they were created during population maintenance by chain sib-mating. In all but one case of observed LD in improved cycle populations (from C1 to C4), the instances of LD were not found in the C0 population and must have been generated over the course of the RS program. The total number of pairs of markers in LD generally increased with selection, which is consistent with the results reported for other populations improved through RS [7,14,17,43]. Theoretically, LD in a population can arise from the intermixture of populations with different allele frequencies, by chance in small populations (random genetic drift), from selection favoring one combination of alleles over another (epistatic selection), or assortative mating [3,58]. On the other hand, hitchhiking can lead to an increase but also to a decrease of LD between two neutral loci linked to a locus experiencing positive directional selection, depending on the position of this locus relative to two neutral loci [59]. In several previous studies, the LD generated during the course of recurrent selection in maize synthetic populations was suggested to result mainly from genetic drift [60], from natural selection for epistatic effects [57], or from selection for epistatic effects [7,9,14,57]. All the above-mentioned evolutionary forces could also be involved in creating LD between loci in the M3S population. In a single selection cycle, genetic drift is expected to generate new LD between different loci regardless of whether they are linked or unlinked. According to Equation (1) (given in Section 2), the generation of drift-related LD for each single cycle of selection is in favor of unlinked loci with the probability of the two randomly selected loci being linked of only 0.08. However, due to positive correlations between the rate of decay of LD and the recombination rate between the two loci [3], we can assume that the rate of LD decay in the M3S population due to the intercrossing of selected progenies and the seed multiplication of cycle populations was lower for linked than for unlinked loci. This can possibly explain the observed surplus of LD pairs, including linked loci, in the present study. Selection for favorable epistatic interactions may have also been involved in generating LD in the M3S population because of the observed overrepresentation of non-neutral pairs among pairs of loci detected to be in significant LD (based on Equation (2) given in Section 2).

5. Conclusions

The present study was undertaken to investigate changes in allele frequencies and population structure in the M3S maize population that was subjected to four cycles of phenotypic recurrent selection for grain yield and disease resistance. The proportion of non-neutral SSR loci detected using the Waples test of selective neutrality varied among single cycles of selection from 16% to 37%. Multiple trait selection and changing methods of selection applied in the M3S population may have caused some discrepancies in the neutrality status of loci among cycles of selection. In addition, competition among multiple “positive” and multiple “negative” alleles at loci under selection (random changes of allele frequencies within selectively non-neutral loci) may have also affected the consistency of the results. In this sense, the multiple allelomorphism of SSR markers can be a constraint, which decreases the power of the neutrality test to detect genomic regions controlling quantitative traits in multiparental synthetic populations used in RS programs. Most previous studies reported changes in maize populations under recurrent selection at the molecular level using SSR markers, whereas only a few recent studies used SNP-based assays. Due to their biallelic nature, SNP markers can overcome the inadequacy of multiallelic SSR markers observed in the present study in terms of unambiguous detection of non-neutral loci/alleles. In addition, due to their high abundance in plant genomes, SNPs are the marker of choice for future studies on the molecular basis of population response to selection. On the other hand, multiallelic SSR markers are much more informative than biallelic SNP markers [61] and are convenient for studying genetic diversity within and among populations, including synthetics, landraces, and inbred lines. The results of the present study showed that the mean number of alleles per locus and average gene diversity (expected heterozygosity) were preserved over cycles of selection, indicating maintenance of sufficient genetic variability in the population required for future genetic gain. Furthermore, several SSR loci declared as non-neutral in the present study have been previously reported to be under selection or share the same genomic locations with previously published QTLs controlling important agronomic traits and can be implemented in marker-assisted breeding programs.