Population genetic structure of Puccinia striiformis f. sp. tritici at the junction of Gansu, Sichuan and Shaanxi Provinces in China

Abstract

Wheat stripe rust, caused by Puccinia striiformis f. sp. tritici, is one of the most destructive diseases of wheat in China. Conjunction area of Gansu, Sichuan and Shaanxi, acting as over-summering and over-wintering regions for the pathogen, plays a unique and critical role in epidemics of this disease in China. Because of the complexity in terrains and environmental conditions within this conjunction area, studies on the population structure, gene flow between local subpopulations and maintenance of genetic structure over time within this area are important to understand the epidemiology of this disease in China, and have practical significance in management of this disease at national scale. In this study, 461 isolates of Pst were collected from the junction area of Gansu, Sichuan and Shaanxi from 2013 spring to 2014 spring, and genotyped with amplified fragment length polymorphism markers. Results revealed that genotypic and genetic diversity were consistently high in Gansu and Shaanxi, but low in Sichuan, and a closer genetic relationship was found between Gansu and Sichuan than between them and Shaanxi illustrated by φpt, shared genotypes, Bayesian and nonparametric clustering methods. Genetic differentiation existed among autumn subpopulations, and genetic barriers were detected, although spring subpopulations were less differentiated. Subpopulations in Gangu of Gansu and Longxian of Shaanxi remained stable over the seasons studied. Potential migration events occurred at the junction area between successive seasons. The estimated frequency of sexual reproduction was 0.970 (s) (i.e. 97% of individuals being sexually derived during the yearly sexual cycle), suggesting the existence of sexual reproduction in this region. The main conclusions of this study are that genetic barriers exist at the junction area, and subpopulations in Gansu and Shaanxi are stable, and population exchange occurs mainly between Gansu and Sichuan.

Background

Wheat stripe rust, caused by Puccinia striiformis f. sp. tritici (Pst), is a destructive disease of wheat throughout the world. It is also one of the most important diseases in China, where the epidemics of this disease in 1950, 1964, 1990, and 2002 caused yield losses up to 6.0, 3.2, 1.8 and 1.3 million metric tons, respectively (Li and Zeng 2002; Wan et al. 2004; Zeng and Luo 2006).

Occurrence of stripe rust in most wheat growing areas in China relies on long-distance dispersal of urediniospores by wind (Li and Zeng 2002; Zeng and Luo 2006) because the pathogen Pst, not forming any resting structure and obligatorily parasitizing on living tissues, cannot survive summer and winter (over-summering and over-wintering) in most places where the average temperature is higher than 23 °C in the summer (Li and Zeng 2002) or where temperature is low in the winter and all expanded wheat leaves are frozen to death (Li and Zeng 2002; Zeng and Luo 2006). Previous studies have identified that Gansu, Sichuan and Shaanxi Provinces acted as part of pathogen’s main over-summering and over-wintering regions in China (Li and Zeng 2002; Zeng and Luo 2006; Wan et al. 2007). In these provinces, south and east Gansu, northwest Sichuan, and west Shaanxi with high elevations and relatively low summer temperatures can serve as over-summering areas where Pst survives on late-maturing wheat and volunteer seedlings in the summer. However, Sichuan Basin with warm winters serves as the largest over-wintering region where urediniospores propagate continuously on wheat throughout winter. The urediniospores of Pst are dispersed down from over-summering areas to the over-wintering areas in autumn (Li and Zeng 2002; Wang et al. 2010; Liang et al. 2016), produce more urediniospores there and then migrate to the spring epidemic region (such as the Huang-Huai-Hai winter wheat region, the largest wheat growing region in China) in the spring when weather conditions there become conducive to infection of wheat by the pathogen (Li and Zeng 2002; Zeng and Luo 2006). The disease epidemic in the Huang-Huai-Hai winter wheat region can be greatly affected by the pathogen population structure (particularly emergence of novel races) in these three provinces (Wan et al. 2007).

Given the importance of these three provinces in epidemics of wheat stripe rust, efforts have been made to understand the pathogen population structure and diversity in this junction region (Wang et al. 2010; Liang et al. 2013, 2016; Liang 2014), which is related to the capacity of pathogens to adapt to new host varieties, and long-distance migration (Burdon et al. 1989; Carlier et al. 1996; Ali et al. 2016). In previous studies, one to several provinces or a large part of a province were considered as a population for analyzing genetic structure, and the local structure was ignored although heterogeneity existed in terrains/environmental conditions, since high mountain, plateau and basin coexist in this conjunction area. High mountains could serve as geographical barriers, which were reported to restrict gene flow between Pst subpopulations (Xie et al. 1993; Liang et al. 2013). Given that the Tsinling Mountains, Daba Mountains, and Liupan Mountains exist in this area, it is reasonable to hypothesize that there are geographical barriers for the migration of Pst at the junction area.

The maintenance of genetic structure over time for these local Pst populations is also very important as it is related to the possibility and rapidity of adaptation in the population (Goddard et al. 2005), and it is particularly important for determining how much diversity can be maintained in traits related to virulence, fungicide resistance or high temperature tolerance of Pst (Ali et al. 2014). The information of population exchanges in temporal scale among regions can be used to infer the direction of gene flow, which is critical for understanding the maintenance of a wide-ranging genetic variation (Zamudio et al. 2009). Furthermore, it is required for formulating better disease management strategies at national scale (Xhaard et al. 2012; Chen et al. 2013).

In order to detect local structures, this area need to be divided into multiple small areas with relatively uniform environmental conditions. The variance within sample units caused by geographical barriers and heterogeneity can be reduced when county rather than province is taken as the population unit. In addition to sampling strategies, multiple and robust population genetic analysis should be used for obtaining reliable results. In previous studies, genetic relationships among geographical or temporal subpopulations were determined by diversity, pair-wise genetic differentiation (Gst, Fst or φpt), phylogenetic analysis, and shared genotypes (Liu et al. 2011; Lu et al. 2011; Liang et al. 2013; Wang et al. 2013; Wan et al. 2015) etc. These methods provide descriptive results to a certain degree for inferring genetic relationships (Liang et al. 2013; Wan et al. 2015), and usually emphasize the variance among subpopulations (Estoup and Guillemaud 2010). However, the variance among individuals cannot be ignored since several studies confirmed that a large proportion of variance exists among individuals of Pst (Lu et al. 2011; Wan et al. 2015). Thus it can introduce errors from defining groups artificially in the analysis of genetic clusters. Cluster methods, including parameter (such as STRUCTURE) and non-parameter (such as discriminant analysis of principle components DAPC) analysis methods that emphasize statistics at the individual level are more robust in identification of population subdivision than descriptive statistics (Pritchard et al. 2000; Jombart et al. 2010).

Except for above methods, additional analysis can be carried out to gain insight into the spatial genetic barriers and temporal maintenance. For the spatial analysis, the correlation between geographical and genetic distance at population level is usually used to estimate the possibility of pathogen’s migration (Hu et al. 2017), but it may not be appropriate to estimate the short-distance dispersal of a long-distance dispersal pathogen (such as Pst) (Li and Zeng 2002). The R package of Geneland is suitable for analyzing the geographical distribution of diversity, and multilocus genotypes (MLGs)/genetic groups can help to test the geographical barriers at individual level (Guillot et al. 2005). For the temporal aspect, the effective population size (Nc) is an important parameter for determining the ability of adaption, and the genetic structure of population with large Nc had less change over time than that with small Nc (Ali et al. 2016; Walker et al. 2017). Previous methods infer Nc based on sequentially sampling populations across generations and do not take the reproduction mode into account (Waples 1989; Wang 2005), while the clonal propagation may introduce bias for Nc estimates in partially asexual populations (e.g. Pst) (Ali et al. 2016). Therefore, Ali et al. (2016) reported a method using R package CloNcaSe that could estimate both the effective population size (Nc) and the frequency of sexual reproduction (s) more accurately (Liang et al. 2013; Wan et al. 2015; Ali et al. 2016).

In this study, a representative set of isolates were collected from multiple counties at the junction of Gansu, Sichuan and Shaanxi Provinces in spring 2013, autumn 2013, and spring 2014, genotyped using amplified fragment length polymorphism (AFLP), and several conventional and contemporary methods in population genetic analysis were carried out: (1) to test whether different genetic groups and geographical barriers exist in the junction area; (2) to analyze the temporal dynamics of local population genetic structure and reproductive mode; and (3) to infer the possible migration events within the junction area.

Results

Effects of the number of loci on the numbers of distinct MLGs and shared genotypes between subpopulations

The AFLP with four pairs of selective amplification primers generated 98 reliable polymorphic bands. According to the presence and absence of these bands, 461 genotypes were detected among 461 isolates and the genotypic diversity was 1.00 for the overall population and each subpopulation. In order to compare genotypic diversity and the number of shared genotypes among subpopulations, only partial loci were chosen for the analysis of the divergence between subpopulations relating to genotypic analysis. When different numbers of loci were chosen randomly from the 98 polymorphic loci, the detected genotypic diversity increased with the number of loci, got close to 1 and varied little when the number of loci increased beyond 30 regardless which loci were selected (Additional file 1: Figure S1), suggesting that using more than 30 polymorphic loci was unnecessary.

Shared genotypes among the three provinces exhibited a consistent pattern regardless of the number of loci used in the analysis, i.e. the number of shared genotypes between Gansu and Sichuan populations was larger than that between Gansu and Shaanxi populations, and that between Sichuan and Shaanxi populations (P < 0.05, Fig. 1). The numbers of shared genotypes were larger using 11 and 13 loci than using 15 and 17 loci, while the standard errors (from sampling of loci) were much smaller using 13 loci than using 11 loci (Fig. 1). Therefore, using 13 loci balanced well between detecting adequate shared genotypes and not being affected dramatically by selection of loci (Fig. 1).

Fig. 1

Number of shared genotypes between paired province populations of Puccinia striiformis f. sp. tritici. Shared genotypes between combinations of Gansu-Sichuan, Gansu-Shaanxi, and Sichuan-Shaanxi were calculated based on the 11, 13, 15, and 17 loci randomly selected from the total of 98 AFLP loci using GenClone version 2.0. Bars represent standard errors from ten independent replicates

Genotypic diversity

The order of genotypic diversities of 3 province populations was Shaanxi > Gansu > Sichuan when 13 loci were used (Table 1). The standard errors were small compared with the means (less than 5% of the means in most cases) (Table 1). In the same seasons of spring and autumn of 2013, subpopulations in Sichuan were less diverse than those in the other two provinces, but their diversities in spring 2014 increased to levels similar to the diversities of subpopulations from Shaanxi and Gansu Provinces (Table 1). Among seasonal subpopulations at the same locations, diversities in subpopulations from Gangu County (GS-P1-S13, GS-P7-A13, and GS-P10-S14) and subpopulations in Longxian County (SX-P9-A13, SX-P15-S14) maintained at high and stable levels (Table 1).

Table 1 Genotypic and genetic diversities of subpopulations of Puccinia striiformis f. sp. tritici

Genetic diversity

When all subpopulations from a province were considered as a single province population, there was no significant difference in genetic diversity among provincial populations no matter whether H or I index was used as the indicator (Table 1). Although genetic diversities of individual subpopulations were different from their genotypic diversities, the overall trend is that genetic diversities of subpopulations from Sichuan in spring 2013 and autumn 2013 were much lower than those in spring 2014. Genetic diversities of subpopulations from Gansu and Shaanxi maintained at higher levels over the three seasons, which was consistent with the result for genotypic diversity.

Spatial population structure analyses

When subpopulations from the same seasons were compared, 12 and 14 shared genotypes were detected between Gansu and Sichuan in spring 2013, and spring 2014, respectively, while other subpopulation pairs, such as GS-S14/SX-S14 and SC-S14/SX-S14 shared only 2 or fewer genotypes (Fig. 2). This implied a close genetic relationship between Gansu and Sichuan populations in the springs of 2013 and 2014.

Fig. 2

Number of shared genotypes between paired seasonal subpopulations from Gansu, Sichuan, and Shaanxi Provinces. S13 = Spring of 2013, A13 = autumn of 2013, and S14 = spring of 2014. Gansu = GS, Sichuan = SC, Shaanxi = SX

The result of analysis of molecular variance (AMOVA) showed that pair-wise genetic differentiation was all significant at P = 0.05 for all pairs of subpopulations collected in the spring season except for the pairs of GS-P3-S13/GS-P4-S13, and SC-P12-S14/SC-P14-S14 in spring of 2013, and SX-P15-S14/SX-P16-S14 in spring of 2014 (Fig. 3). It implied that significant genetic divergence existed in most population pairs. According to Hartl and Clark (1998), the value of genetic differentiation φpt > 0.25 implies a significant genetic difference between two populations. All the φpt values (genetic differentiation) between paired subpopulations in autumn of 2013 were above 0.25, but quite a few pair-wise φpt values were low between subpopulations in the springs of 2013 and 2014, especially for pairs of subpopulations collected from the same province (Fig. 3), suggesting that the levels of genetic differentiation between autumn subpopulations were higher than those between spring subpopulations.

Fig. 3

Pair-wise φpt among 17 subpopulations of Puccinia striiformis f. sp. tritici in this study. The larger the circle the larger the φpt value was, and the black, dark gray and light gray colors of circle represent the φpt level range of 0–0.25, 0.25–0.5, and 0.5–0.75, respectively. A “X” indicated the differentiation between the paired subpopulations was insignificant at P = 0.05

The multivariate DAPC suggested that K = 3 was the optimal number of clusters according to the Bayesian Information Criterion (BIC) (Additional file 2: Figure S2 and Fig. 4a). Using the threshold of membership assignment possibility > = 80%, individuals collected from Gansu in spring of 2013 were mainly assigned to groups G1 and G3 (Fig. 4b), and samples from Sichuan Basin were mostly assigned to G3 (Fig. 4c). Similar to the results of genetic differentiation and shared genotypes, subpopulations from the three provinces differed remarkably on genetic composition in autumn of 2013 (Fig. 4e) in that individuals from Gansu (GS-P7-A13) were assigned to G1, individuals from northwest Sichuan (SC-P8-A13) and Shaanxi (SX-P9-A13) were mainly assigned to G3 and G2, respectively (Fig. 4c). However, a different situation was observed in the spring of 2014 when most samples, regardless of their origin locations, were classified into G2 with only a few isolates in G1 and admixture group (Fig. 4d), suggesting a high level of genetic identity among all subpopulations in the spring of 2014.

Fig. 4

Results from discriminant analysis of principle component (DAPC) and STRUCTURE analysis. a Scatter-plot of 461 isolates of Puccinia striiformis f. sp. tritici according to DAPC. The symbols of “△”, “●”, and “▲” represented isolates belonging to genetic groups 1, 2, and 3, respectively. b, c and d The frequency and geographical distribution of isolates assigned to genetic groups G1, G2, G3 and admixture (A) in subpopulations collected in spring 2013, autumn 2013, and spring 2014, respectively, according to DAPC. e Bayesian assignment of individual isolates into three genetic groups (black for G1, grey for G2, and white for G3) in STRUCTURE analysis. f, g and h The frequency and geographical distribution of isolates assigned to groups G1, G2, G3 and A in subpopulations collected in the spring of 2013, autumn of 2013 and spring of 2014, respectively, according to STRUCTURE analysis. A individual was assigned to a group if it had a probability no less than 80, and 60% belonging to the group in DAPC, and STRUCTURE analysis, respectively, otherwise it would be assigned to group A. Population code (P1-P17) was abbreviated from GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13, SC-P5-S13, SC-P6-S13, GS-P7-A13, SC-P8-A13, SX-P9-A13, GS-P10-S14, SC-P11-S14, SC-P12-S14, SC-P13-S14, SC-P14-S14, SX-P15-S14, SX-P16-S14, and SX-P17-S14, respectively

Clustering analyses based on Bayesian method with STRUCTURE software also identified a clear pattern of population subdivision. According to the method used by Pritchard et al. (2000), K = 3 was the optimal value for the number of clusters based on the largest rate of change between partitions (ΔK) (Additional file 3: Figure S3). Although more isolates were classified into admixture group than those in DAPC analysis (Fig. 4e-h), proportions of individuals belonging to the rest genetic groups were generally consistent with the results of DAPC.

Assuming uncorrelated allelic frequencies between sites, four clusters were identified in each spring population using Geneland program in R 3.4.2 (Additional file 4: Figure S4 and Fig. 5). In the spring of 2013, significant genetic boundaries existed between Sichuan and Gansu except between GS-P2-S13 and SC-P6-S13 (Fig. 5a). In the spring of 2014, individuals from a province were classified into the same genetic cluster, except that SX-P15-S14 and SX-P17-S14 from Shaanxi were divided into two genetic clusters 1 and 2 (Fig. 5f, g). These results suggested that genetic barriers existed at the junction of Gansu, Sichuan and Shaanxi Provinces, and that the provincial isolation of the three provinces were basically coincided with the genetic boundaries.

Fig. 5

Maps of the posterior probabilities of population membership inferred by Geneland. a-d Four clusters inferred by Geneland among individuals in spring 2013. e The map of study area in this study. f-i Four clusters in spring of 2014. Contour lines indicate the spatial position of the genetic discontinuities. Lighter shading indicates higher probabilities of population membership. Due to the limited space, population code (P1-P6, P10-P17) was abbreviated from GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13, SC-P5-S13, SC-P6-S13, GS-P10-S14, SC-P11-S14, SC-P12-S14, SC-P13-S14, SC-P14-S14, SX-P15-S14, SX-P16-S14, and SX-P17-S14, respectively

Temporal dynamics

When subpopulations collected from the same sampling sites (or sites very close to each other) in successive seasons were compared, large numbers of shared genotypes (Fig. 2) were detected between spring and autumn of 2013, between autumn of 2013 and spring of 2014 for samples from Gansu, and between autumn 2013 and spring 2014 for samples from Shaanxi. Pair-wise φpt between GS-P1-S13 and GS-P7-A13, between GS-P7-A13 and GS-P10-S14, between SX-P9-A13 and SX-P15-S14 were all smaller than 0.25 (Fig. 3), which was similar with the results of DAPC and STRUCTURE (Fig. 4), and suggested a certain degree of population stability in Gangu and Longxian Counties. In addition, the effective population size (Nc) in Gansu (contained GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13, and GS-P10-S14) and Sichuan Basin (contained SC-P5-S13, SC-P6-S13, SC-P11-S14, SC-P12-S14, and SC-P13-S14) over years were 479 (327–748) and 706 (406–1515) (Table 2), suggesting that both populations had a high survivability over years (Schmeller et al. 2002). The frequency of sexual reproduction (s) was 0.970 (0.858–1), revealing that high proportion of individuals undergo sexual reproduction in the junction area.

Table 2 Estimates of the frequency of sexual reproduction s and effective population size Nc

For the analysis of different locations in successive seasons, between the spring and the autumn of 2013, the population pairs of GS-S13/SC-A13, SC-S13/SC-A13, had 3 and 4 shared genotypes, respectively (Fig. 2). It showed a relatively low φpt between SC-P8-A13 and the other 6 subpopulations in spring 2013 (Fig. 3). Coupled with above results, results of DAPC and STRUCUTRE all showed high similarity of subdivision patterns between subpopulations from Gansu (GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13) and Sichuan Basin (SC-P5-S13, SC-P6-S13) in spring 2013 and SC-P8-A13, the subpopulation from Sichuan in autumn 2013 (Fig. 4b, c, f, g). These results suggested that Pst could migrate from Gansu and Sichuan Basin to Northwest Sichuan in spring. Between the autumn of 2013 and the spring of 2014, population pairs of GS-A13/SC-S14, GS-A13/SX-S14, SX-A13/GS-S14, and SX-A13/SC-S14 all had one genotype shared respectively (Fig. 2). Results of φpt, DAPC, and STRUCTURE showed similar results, revealing that population exchange occurred between subpopulations of Gansu, Shaanxi in autumn 2013 and subpopulations of Gansu, Sichuan and Shaanxi in the spring of 2014.

Discussion

The results from our genetic structure analysis demonstrated that there were barriers to migration of Pst at the junction area of Gansu, Sichuan and Shaanxi despite Pst had the capacity for long-distance migration. Previous study focusing on large geographic range also inferred the existence of genetic barriers, and attributed it to the long geographical distance (Liu et al. 2011; Hu et al. 2017). Apparently geographical distance was one of the possible reasons, but played a weak role in explaining the genetic differentiation in spring, as the geographical distance was not far among subpopulations at the junction area. Thus, the local environment might be another important factor leading to genetic divergence at the junction area (Crispo et al. 2006; Schluter 2009; Orsini et al. 2013). Genetic divergence was detected at different agricultural ecological region within Longnan, revealing that the environment could affect the population structure. At the junction area, Tsinling Mountains, Liupan Mountains, Wei River, and Sichuan Basin could make the terrain and climate complex, therefore the main environmental factors (e.g. elevation, humidity, habitat or host) influencing the migration of Pst should be confirmed through all year around monitoring in the future.

Our results revealed that although mountains in the junction area can serve as the barriers to migration, and high genetic differentiation exists between autumn subpopulations in this areas, the spring subpopulations at the conjunction area tended to be undifferentiated or unified, inferred based on low pair-wise φpt, similar genetic structure, and a large number of shared genotypes. The large numbers of shared genotypes detected between provinces in two spring seasons from this study suggested the existence of strong population exchange in the spring. This was also supported by the observations of large amount of urediniospores on wheat leaves in the warm regions like Sichuan Basin (Li and Zeng 2002) and the dominant southeast wind direction within this junction area during spring (Chen et al. 2014). On the other hand, because only a small number of urediniospores survive summer, genetic drift and founder effect may be the main reasons for the high genetic differentiation in the autumn subpopulations in this area.

Genotypic diversities were at similar levels for subpopulations from Gansu and Shaanxi Provinces and higher than those from Sichuan Province except for spring 2014. This was consistent with the facts that the sexual stage of Pst was reported on Berberis spp. in Gansu and Shaanxi Provinces (Mboup et al. 2009; Jin et al. 2010; Wang et al. 2016; Zhao et al. 2016, 2018) and that Southeast Gansu and Southwest Shaanxi served as the key over-summering area and the origin area of new races for Pst (Li and Zeng 2002; Chen et al. 2013). Genotypic diversity in GS-P2-S13, GS-P3-S13 and GS-P4-S13 was much lower than that in GS-P1-S13, which might be an explanation for the whole diversity in Gansu was lower than Shaanxi. Our study showed that the junction area had high values of the frequency of sexual reproduction (s), genotypic diversity and genetic diversity, all suggesting the existence of sexual reproduction. The genotypic diversity of SX-P9-A13 (0.973) was significantly greater than that of SX-P15-S14 (0.896), seemed not coincident with the fact that sexual reproduction occurred in early spring (Zhao et al. 2018). Although Pst in Longxian County could over-winter at low altitude, it could not be ruled out that some individuals died due to the low temperature in winter. That might be an explanation for reduced diversity in Longxian from autumn to the next spring season.

The results of shared genotypes among the three provincial populations revealed that Gansu population had a closer genetic relationship with Sichuan population than with Shaanxi population, which was consistent with Liang’s results (2014). High level of genetic similarity between Gansu and Sichuan populations might be explained by frequent genetic exchange between them (Li and Zeng 2002; Zeng and Luo 2006; Wan et al. 2007; Liang et al. 2016). Although urediniospores of Pst could be dispersed from the west to east through the upper air current in autumn (Li and Zeng 2002), the Tsinling and Liupan Mountains might limit the gene flow between Gansu and Shaanxi population to a certain degree. Liupan Mountains had been reported to play a barrier role on the gene flow between Gansu and Ningxia populations (Liang et al. 2013), which might also be an explanation for genetic differentiation between Gansu and Shaanxi populations.

Instability of genetic structure among seasons at the same location was revealed by significant pair-wise φpt (P < 0.05) between seasonal subpopulations, while low levels of genetic drift were also suggested by the pair-wise φpt values smaller than 0.25, indicating the possibility of local survival during over-summering or over-wintering. In addition, similar values of Nc in Gansu Province and the overall junctional population implied that Gansu population played a key role in maintaining population stability over time. Subpopulations of Gangu in three successive seasons (GS-P1-S13, GS-P7-A13 and GS-P10-S14) shared several genotypes, and exhibited low genetic differentiation (φpt < 0.25), indicating Pst might have survived locally during over-summering and over-wintering periods in Gangu County. Similarly, the autumn subpopulation and the subpopulation in the next spring had common genotypes, and only 3% variance (Additional file 5: Figure S5) existed between seasonal subpopulations in Shaanxi proved local overwintering of Pst in Shaanxi, which was consistent with previous studies (Li and Zeng 2002; Wan et al. 2007).

Previous studies had determined that close genetic relationships existed among Gansu, Sichuan and Shaanxi Provinces, however detailed migration events in different seasons were not clear. Common genotypes and low level of genetic differentiation of Pst at different locations in successive seasons might offer evidence for migration (Liang et al. 2013, 2016). Diseased leaves were collected from different seasons, thus the direction of migration could be inferred as from the early sampling area to the late sampling area. In our study, results of φpt and shared genotypes illustrated that possible migration events were from Gansu to Northwest Sichuan, and from Sichuan Basin to Northwest Sichuan in spring of 2013, and from Gansu to Sichuan Basin, from Gansu to Shaanxi, from Shaanxi to Gansu, from Shaanxi to Sichuan Basin in autumn of 2013. These results were consistent with previous studies (Li and Zeng 2002; Liang et al. 2016). The findings of no common genotypes and great genetic differentiation between Northwest Sichuan (SC-P8-A13 from Heishui) and Sichuan Basin (SC-P11-S14, SC-P12-S14, SC-P13-S14, and SC-P14-S14) implied a low possibility of migration from Northwest Sichuan to Sichuan Basin in autumn. It was different with previous inference that Pst on autumn seedlings in the northwest Sichuan could be blown to Sichuan Basin under the impact of westerly circulation (Li and Zeng 2002; Zeng and Luo 2006; Wan et al. 2007). The northwest plateau of Sichuan consisted of high mountains with average altitude of 3000–4000 m. Previous study found Pst in Heishui had a certain genetic distance with the other populations from the northwest Sichuan, illustrated a spatial population subdivision (Wang et al. 2010). Further studies are necessary to assess which exact (local) subpopulations of the northwest Sichuan have important effect on Sichuan Basin population.

Sampling size was a vital impact factor for analyzing population genetic structure, especially for population with a high level of variation (Goss 2015). Thus future studies should emphasize the number (greater than 30 for each population) and representativeness of samples (Lu 2009) to include more sites representing the population, varieties of host cultivars and altitudes etc. Reliable results can be generated through extensive polymorphic bands amplified by AFLP, which is not available for other methods such as simple sequence repeat (SSR) (Powell et al. 1996; Liang 2014), although AFLP cannot distinguish heterozygous genotypes from homozygous genotypes for the dominant allele. Nevertheless, for diploid organisms, codominant markers can generate more accurate clustering results than dominate markers, such as AFLP, when limited number of loci were used (Guillot and Santos 2010). Therefore, codominant markers such as single nucleotide polymorphism (SNP) and simple sequence repeat (SSR) should be used in the future research as much as possible. In addition, how specific environmental forces (e.g. terrain and climate) shape the genetic structure and interregional migration of Pst in this study remain uncertain and will require additional research to resolve.

Conclusions

Pst in Gansu and Sichuan had a closer genetic relationship compared with Shaanxi Province, and pathogens in Gansu could migrate to Sichuan Basin in autumn, thus the disease incidence in Gansu could serve as a risk estimator for Sichuan Basin. A cooperative scheme for management of wheat stripe rust can be formulated among the three provinces. The time for disease control (fungicide application) in Gansu and Sichuan should be set in autumn and early spring, respectively.

Methods

Sample collection

In the spring and autumn of 2013, and the spring of 2014, diseased wheat leaves showing sporulation of Pst were collected from the junction of Gansu, Shaanxi and Sichuan Provinces (Fig. 6), including 4 counties of Gansu Province, 7 counties of Sichuan Province and 3 counties of Shaanxi Province (Additional file 6: Table S1). Samples collected from the same county within one season were considered as one subpopulation. Seventeen subpopulations were collected: 6 from Gangu, Chengxian, Huixian, Qinzhou, Jiange, and Zitong Counties in the spring of 2013, named as GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13, SC-P5-S13, and SC-P6-S13, respectively; 3 from Gangu, Heishui, and Longxian Counties in the autumn of 2013, named as GS-P7-A13, SC-P8-A13, and SX-P9-A13, respectively; 8 from Gangu, Shehong, Daying, Nanbu, Jialing, Longxian, Qishan and Meixian Counties in the spring of 2014 named as GS-P10-S14, SC-P11-S14, SC-P12-S14, SC-P12-S14, SC-P13-S14, SC-P14-S14, SX-P15-S14, SX-P16-S14, and SX-P17-S14, respectively (Additional file 6: Table S1). Subpopulations were named as a two-letter brief of province name (SC for Sichun, SX for Shaanxi and GS for Gansu)-population number (P1 to P17)- season (A for autumn and S for spring) and year (13 and 14 for 2013 and 2014, Additional file 6: Table S1). Among subpopulations in Sichuan Province, SC-P5-S13, SC-P6-S13, SC-P11-S14, SC-P12-S14, SC-P13-S14 and SC-P14-S14 were sampled in Sichuan Basin, while SC-P8-A13 was sampled in the northwestern mountain area. A total of 461 isolates of Pst were obtained from diseased wheat leaf samples collected at the junction of Gansu, Sichuan and Shaanxi Provinces, of which 186 isolates were obtained in the spring of 2013, 51 isolates in the autumn of 2013, and 224 isolates in the spring of 2014 (Additional file 6: Table S1).

Fig. 6

Geographical locations where isolates of Puccinia striiformis f. sp. tritici were collected in this study. a The overall view of the study area, which was marked with a square frame, on the map of China; as well as the enlarged views of the study area showing the locations of samples collected in: b the spring of 2013, c the autumn of 2013, and d the spring of 2014. Due to the limited space, population code (P1-P17) was abbreviated from GS-P1-S13, GS-P2-S13, GS-P3-S13, GS-P4-S13, SC-P5-S13, SC-P6-S13, GS-P7-A13, SC-P8-A13, SX-P9-A13, GS-P10-S14, SC-P11-S14, SC-P12-S14, SC-P13-S14, SC-P14-S14, SX-P15-S14, SX-P16-S14, and SX-P17-S14, respectively

Preparation of genomic DNA from samples

Urediniospores of the 186 isolates collected in spring, 2013 were propagated on seedlings of susceptible wheat cultivar Mingxian 169. About 8 seedlings were planted in a 1 L pot. After moisturized at 10 °C -13 °C for 12 h, urediniospores were picked up from a single uredium pustule on a leaf sample and rubbed in a drop of water on the leaf surface of 10-day-old Mingxian 169 seedlings for inoculation. The inoculated seedlings were incubated in a moist chamber at 8 °C to 10 °C for 24 h in darkness to promote infection, and transferred to an incubator at 17 °C and 14 °C (day and night, respectively) with 14 h of light per day. After about 13 days of incubation, urediniospores were harvested from the lesions and then reinoculated in the same way as described above to further increase spores. At the end, spores were harvested from inoculated leaves, transferred into a 0.5 mL Eppendorf tube, dried in desiccators at 4 °C for 3 to 4 days, and stored at − 20 °C for DNA extraction. Genomic DNA of urediniospores obtained from propagation artificially was extracted using a modified CTAB method (Justesen et al. 2002). The DNA concentrations and qualities of these samples were determined via UV absorption at wave-length 260 and 280 nm using a Nanodrop 2000 Spectrophotometer (Gene Company Limited, China). The DNA samples were then diluted to a concentration of 100 ng/μL and stored at − 20 °C for AFLP.

For the 51 samples collected in the autumn of 2013 and the 224 samples in the spring of 2014, genomic DNA were amplified directly from samples using a modified multiple displacement amplification (MDA) method, which was demonstrated to amplify the genome DNA of Pst well enough for further AFLP analysis without introducing significant errors in our study published previously (Zhang et al. 2015). Urediniospores from a single pustule on the sampled leaves were picked to make a suspension at a concentration of 20–30 spores/μL in phosphate-buffered saline. Then 1 μL spore suspension was transferred into a 200 μL centrifuge tube and used as MDA template. The MDA process was conducted using REPLI-g Mini Kit (Qiagen Company, Beijing, China) following the method used by Zhang et al. (2015). All amplification products obtained via MDA method were also diluted to the same concentration of 100 ng/μL and stored at − 20 °C for AFLP analysis.

AFLP procedure

All the 461 DNA samples obtained via CTAB and MDA methods were subjected to AFLP analysis following the procedure described by Justesen et al. (2002). DNA was digested with restriction endonuclease PstI and MseI, and ligated with T4-DNA ligase in a 12.5 μL reaction system containing 350 ng DNA, 1.5 U MseI, 2 U PstI, 8 U T4-DNA ligase, 5 pmol PstI-adapter, 50 pmol MseI-adapter (Additional file 7: Table S2) and 2.5 μg NEB-BSA buffer. The reaction was performed at 37 °C for 10 h and then inactivated by heating the sample at 65 °C for 20 min. Then pre-amplification was performed using the primer pair of PstI0 and MseI0 (Additional file 7: Table S2) in a 20 μL reaction system containing 2 μL ligated DNA, 4 ng primers PstI0 and MseI0, 1 U Taq polymerase, 1.6 μL dNTP (2.5 mM) and 2 μL 10 × PCR buffer (Mg²⁺ included). The pre-amplification was performed under the following conditions: 94 °C for 5 min, followed by 30 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 60 s. The products from pre-amplification were diluted for 30 times with sterilized double distilled H₂O and then used for selective amplification. Four pairs of selective primers (Additional file 7: Table S2) were used and the reaction systems were set up the same as the pre-amplification reaction described above. The selective amplification was performed with the following steps: 94 °C for 5 min, followed by 13 cycles of 94 °C for 30 s, 65 °C for 30 s with a 0.7 °C decrement per cycle and 72 °C for 60 s, followed by 25 cycles of 94 °C for 30 s, 56 °C for 30 s and 72 °C for 60 s. Final PCR products were separated and visualized using an ABI3730XL automatic DNA sequencer (Applied Biosystems, Carlsbad, CA, USA) by Beijing Tsingke Biotech Co., Ltd. AFLP fragments within the range of 100 to 500 bp were scored for each sample visually using GeneMarker (version 2.4.0, SoftGenetics, LLC. 100 Oakwood Ave, Suite 350 State College, PA 16803 USA). Scores were assigned as 1 if the band was present and 0 if the band was absent.

Data analysis

First, the effects of the number of AFLP loci (from 1 to 40) on the number of distinct MLG and on the number of shared genotypes detected between subpopulations were examined using GenClone (version 2.0, Arnaud-Haond and Belkhir 2007). To further determine the proper number of loci in analyses of shared genotypes, 11, 13, 15 and 17 loci (which could differentiate 58, 72, 81, and 87% of total genotypes) were sampled randomly from all polymorphic loci 10 times, and the means and standard errors for shared genotypes between province subpopulations in the same or successive seasons were calculated. The significance of variance in shared genotypes among different population pairs was estimated in Statistical Package for Social Sciences (IBM Corp. Released 2017. IBM SPSS Statistics for Windows Version 25.0. Armonk, NY: IBM Corp.), and a P < 0.05 was used for statistical significance of differences. Based on the results of this analysis, 13 loci were used for subsequent analysis to balance between detecting a large number of shared genotypes among subpopulations and having small errors from selection of loci. Second, genotypic diversity within a population/subpopulation was calculated in GenClone (Arnaud-Haond and Belkhir 2007) for each selection of 13 loci. The mean was calculated from the values in the 10 random sets of loci, and the standard error calculated as standard deviation divided by the square root of 10.

The genetic diversity within each sampled subpopulation, provincial population and total population was assessed by using two indices, Nei’s gene diversity (H) (Nei 1972) and Shannon’s information index (I) on each of the all 98 loci in POPGENE (version 1.3, Yeh et al. 1997). The mean values and standard deviations of these diversity indexes were then calculated from the results of the 98 loci, and the standard errors were calculated as the standard deviations divided by the square root of 98.

To gain insight into the spatial population genetic structure, multiple approaches were taken. First, the number of shared genotypes between interprovincial subpopulations in the same season and successive seasons were calculated using 13 loci. Distribution of each genotype was displayed in GenClone (Arnaud-Haond and Belkhir 2007), and the shared genotype between seasonal interprovincial subpopulations was performed by manual statistics. Second, pair-wise genetic differentiation was examined between the subpopulations of the same season and successive seasons by using the AMOVA program in GENALEX (version 6.5, Peakall and Smouse 2006, 2012). The parameter φpt, an analogue of Fst, was estimated for each pair-wise comparison with 999 permutations based on Nei’s genetic distance matrices. If the φpt value differed significantly from zero, the null hypothesis of no genetic differentiation could be rejected, and if φpt was greater than 0.25 and statistically significant (P < 0.05), then the two subpopulations could be considered as two populations (Hartl and Clark 1998). φpt and P matrix were displayed in R 3.4.2 using corrplot package to make the result more visual. Third, multivariate analyses were carried out using DAPC, implemented in the ADEGENET package in R environment (Jombart 2008). The number of clusters was identified based on BIC as suggested by Jombart et al. (2010) Individuals were assigned to a population if they had a probability of > = 80% belonging to that population. Fourth, a model-based Bayesian clustering method was used for inferring the number of genetic groups (K) and evaluating the extent of admixture among them implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000). The Monte Carlo Markov Chain (MCMC) sampling scheme was run for 400,000 iterations with a 100,000 burn-in period, with K ranging from 1 to 10 and 15 independent replications for each K. The best K value were obtained in STRUCTURE HARVESTER (http://taylor0.biology.ucla.edu/ structureHarvester) valuated by ΔK as proposed by Evanno et al. (2005). Individual files and population files of 15 runs corresponding to the best K were processed with CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007), then bar plot was generated in Distruct 1.1 (Rosenberg 2004) using the output of CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007). An individual was categorized into a given population if they had > = 60% probability of assignment to that population. Finally, for further detailed analysis, the Bayesian clustering method implemented in the R package Geneland version 4.0.4 (Guillot et al. 2005) was used to detect geographical discontinuities in genetic clustering in the springs of 2013 and 2014, but not for subpopulations in autumn 2013 because there were not adequate sampling sites in the autumn of 2013. Ten independent runs were performed with 1,000,000 MCMC iterations, of which every 1000th one was saved. Initially, the tested number of genetic clusters (K) was set to vary between 1 and 6 (in spring 2013) or between 1 and 8 (in spring 2014). The run with the highest average posterior probability was chosen for subsequent analysis. The relationship of spatial domain and population membership was calculated using a burn-in period length of 100 iterations.

Temporal dynamics were assessed by measuring shared genotypes, pair-wise φpt, and genetic components clustered by DAPC and STRUCTURE between the 3 seasons. The local temporal maintenance was estimated by the contemporary effective population size (Nc), a method based on clone-mate resampling probability within temporally spaced samples and displayed using R package CloNcaSe (Ali et al. 2016). The effective population size was estimated for subpopulations at each location and for the overall population, in which Pst could complete annual cycle (spring of 2013 to spring of 2014), thus subpopulations in Shaanxi were not involved in calculation. The same 13 loci as in the analysis of shared genotypes were used in CloNcaSe (if 98 loci were used, the output of CloNcaSe was out of threshold). In addition, the frequency of sexual reproduction (s) was also output at the same time using CloNcaSe (Ali et al. 2016).