A total of 21 P. juncea accessions were selected as test materials, 19 accessions of which were provided by the U.S. National Plant Germplasm Resources Conservation System (NPGS), and two accessions by the China National Medium-term Gene Bank for Forage Germplasm. The variety 549118 ‘Mengnong No. 4’ was a variety selected from the cultivar ‘Bozoisky’ which was introduced from USDA-ARS-FRRL in 1984, and ‘Mengnong No. 4’ was registered by Inner Mongolia Agricultural University in 2009. The seeds of test materials were germinated in a greenhouse in January 2018 and transplanted to the field to establish the test plot in July 2018. The information of all accessions including eco-geographical conditions are shown in Supplementary Table 1. The test plot was located at the forage research station of Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, China (111.41°E, 40.48°N). At least 30 plants in each accession were planted individually at spacings of 50 cm.

A total of 480 individual plants were selected from 21 P. juncea accessions, ranged from 14 to 30 of each. The leaves were wiped with alcohol and deionized water and 1–2 fresh leaves per plant were cut into a 5 mL cryopreservation tube, and brought to the laboratory and stored at −80°C. DNA was extracted using a plant genomic DNA extraction kit (Tiangen Biotech Co., Ltd., Beijing, China). The quality and quantity of DNA was evaluated on 1% agarose gel with a NanoDropTM2000 Spectrophotometer (Thermo Fisher, Waltham, MA, United States). After all the measurements, the DNA stock solution was diluted to 50 ng/μL for SSR marker analysis and stored at −20°C.

In order to establish EST library under the maximum genetic background, 17 accessions with significant difference in tiller number trait between individuals were selected from total 21 accessions. A group of thirty plants with dense tillers and another group of thirty plants with sparse tillers were selected from same 17 accessions (named sample 1 and sample 2, Supplementary Table 2). The normality tests of tiller related traits of sample one and sample two were listed in Supplementary Table 7 and simple statistics analysis were conducted using Origin 2019b (MicroCal, Northampton, MA, United States). 100 mg fresh tissues of the shoot meristem of each plant were briefly washed in deionized water and frozen in liquid nitrogen. Total RNA was extracted from tissue using the rapid plant RNA extraction kit and was treated with RNase-free DNase I at 37°C for 30 min to remove residual DNA (Takara, Otsu, Shiga, Japan). The RNA quality was verified using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) and confirmed using RNase-free agarose gel electrophoresis. The concentration of RNA was quantified with a Nano-drop 2000 (Thermo Fisher Scientific, Wilmington, DE, United States). Equal amounts of RNA of the thirty plants of each group were pooled to establish sample 1 and sample 2, and each sample had three replicates.

The cDNA library was constructed and sequenced by the Biomarker Biotechnology Corporation (Beijing, China). Briefly, the poly (A) mRNA was enriched via magnetic oligo (dT) beads and then broken into short fragments using an RNA Fragmentation Kit (Beckman Coulter, Brea, CA, United States). These cleaved mRNA fragments were used as templates for first strand cDNA synthesis using random hexamer primers. Second strand cDNA was synthesized with RNase H and DNA polymerase I. The library fragments were purified with an AMPure XP system (Beckman Coulter, Beverly, MA, United States). These short fragments were ligated to sequencing adapters, and the desired fragments were separated using AMPure XP beads. Next, the purified cDNA fragments were enriched via PCR. Finally, the cDNA libraries were sequenced using Illumina HiSeq™ 2000 and paired-end reads were generated. After obtaining high-quality sequencing data, the data was assembled to obtain unigenes.

According to the minimum length standard (unit size/minimum repetition time 1/12, 2/6, 3/5, 4/5, 5/5, and 6/5), unigenes over 1 Kb were analyzed using MISA¹ to identify Simple Sequence Repeats. The primer pair for each SSR was designed according to differentially expressed genes and using Primer 3² with the following criteria: primer lengths of 16–23 bases (optimum 20 bp), GC content of 40–70%, annealing temperature of 50–60°C, and PCR product size of 100–500 bp. The primers were synthesized by Sangon Biotech (Shanghai, China).

According to the relative expression level between the two samples, the differentially expressed genes were screened. A total of 2,020 significant differential gene expressions (DEGs) were identified, and 400 significant DEGs with SSR loci were screened. One SSR locus could synthesize three pairs of primers, and one pair of primers was selected for each locus. After removing 75 unigenes that were not suitable for designing primers, 325 pairs of primers were synthesized by Sangon Biotech (Shanghai, China). Eight individual plants from different populations were selected randomly, and used to test the amplification effect and polymorphism of 325 pairs of primers. A total of 103 primer pairs with polymorphism were selected and used in population genetic diversity and structure analysis (Supplementary Table 3).

A total of 103 selected SSR primer pairs were used in 480 individual plant DNA PCR reactions. The PCR reaction system used a 15 μL system for amplification, and reactions containing of 7.5 μL of 2 × Taq PCR MasterMix (including 3.0 mM MgCl₂, 500 μM dNTPs, 0.1 U Taq DNA polymerase, 100 mM KCl, 20 mM Tris–HCl, Tiangen Biotech, Co., Ltd, Beijing, China), 1 μL of each forward and reverse primer, 50 ng of genomic DNA and 4.5 μL of ddH₂O. PCR amplification was performed using an Applied Biosystems (Life Technologies Holdings Pte. Ltd, Singapore). The PCR conditions were as follows: pre-denaturation at 94°C for 4 min followed by 33 cycles of 94°C for 30 s, 51.5–58.5°C (depending on the Tm of the primer set used) for 30 s, and 72°C for 45 s, and finally extension at 72°C for 10 min. The fragment size of PCR amplified product was detected using a QIAxcel Advanced capillary electrophoresis instrument. The 100 bp to 2.5 kb DNA Size Marker and 15 bp to 3 kb Alignment Marker were selected as standard size markers (QIAGEN, Dusseldorf, Germany).

The BioCalculator software was used to accurately analyze a single data of amplified product by calculating the characteristics of various peaks, such as peak number, peak height, peak width and peak area. The DNA fragments amplified by the selected 103 pairs of primers were arranged in order from small to large, with high statistical resolution and clear bands, with “1” indicating presence of fragment or “0” indicating absence of fragment.

Taking band frequency less than 95%, the number of polymorphic bands (NPB) and the percentage of polymorphic bands (PPB) were calculated using Excel 2013. The GeneAlEx 6.503 (Peakall and Smouse, 2012) software was used to convert various file formats for different analysis and to calculate genetic diversity parameters, including the number of alleles (Na), the effective number of alleles (Ne), Shannon’s information index (I), Nei’s gene diversity (h), the number of private alleles, gene flow (Nm), and the F-Statistics (Fst). Genetic distance analysis, principal coordinate analysis (PCoA) and analysis of molecular variance (AMOVA) were conducted using GeneAlEx 6.503 software. AMOVA was carried out to partition the genetic variances into three levels: among populations, among individuals within populations, and within populations. The unweighted pair-group method with arithmetic means (UPGMA) cluster analysis was performed based on Nei’s unbiased genetic distance matrix with MEGA (Kumar et al., 2018). The Mantel test of the association between genetic distance and geographic distance were calculated by Geographic Distance Matrix Generator_v1.2.3 software (Ersts, 2010; Zhou et al., 2020). The geographical distance between different sites was calculated according to the latitude and longitude using the decimal degrees formula³ to detect whether the IBD (isolation by distance) patterns exist (Ersts, 2010).

The population genetic structure was analyzed using STRUCTURE 2.3.4 (Pritchard et al., 2000) software, using a model-based clustering algorithm that implements a Bayesian framework and the Markov Chain Monte Carlo (MCMC) algorithm. To confirm the optimum number of subpopulations (K), five independent runs for each value of K ranging from 2 to 10 were conducted (Wang M. L. et al., 2009). Each run consisted of a burn-in period of 10,000 steps followed by 100,000 MCMC iterations. The ΔK parameter, based on the rate of change in the log probability of data between successive K values, was estimated to determine the best K, based on the model developed by Evanno et al. (2005).

The cDNA library was sequenced using the Illumina Hiseq high-throughput sequencing platform, and the transcriptome sequencing of six samples was completed. After removal of low-quality reads, a total of 38.81 GB of high-quality clean data was obtained. The total base number of clean data was 6,468,655,828, and the average Q30 and GC percentages were 92.37 and 54.02%, respectively. These clean reads were de novo assembled into 100,560 unigenes by the assembling program Trinity, and the total length of the unigenes was 91,562,881 bp, with an average length of 910 bp and the N50 of unigenes was 1,739 bp. Among all the unigenes, those with length ranging from 200 to 1,000 bp accounted for 70.28%, and 29.72% of the unigenes were more than 1,000 bp (Table 1). The length distributions of the unigenes are shown in Table 2 and Figure 2.

A total of 9,300 (31.12%) unigenes above 1 kb were identified as containing SSRs. EST-SSRs were of perfect repeat types, ranging from mononucleotide to hexanucleotide and among them, mononucleotide and trinucleotide repeats accounted for the largest proportion, at 40.30 and 36.11%, respectively and 400 EST-SSRs with compound repeats accounted for 4.3% of all sequences (Table 3).

In the PCR reaction, a total of 344 loci were amplified using 103 SSR primer pairs among 480 individuals, of which 308 (89.53%) exhibited polymorphic patterns. The number of polymorphic loci for each primer combination varied from 1 to 6, with an average of 2.990 loci (Table 4). The amplified fragments of the 103 primer pairs were different in size and varied between 100 and 500 bp. All the primer pairs had a high Nei’s gene diversity value (h) and identified a high level of polymorphism. The percentage of polymorphic loci revealed different levels of polymorphisms ranging from 50 to 100%. The h value varied from 0.092 to 0.449 with an average of 0.274. The average effective Na per locus was 1.898. The primers also showed high Shannon’s information index (I), and the Shannon’s index of the 480 individual plant materials ranged from 0.173 to 0.640, with an average value of 0.420.

Using the amplified markers of 480 individuals as the original data, genetic similarity coefficients were obtained (Supplementary Table 4). The genetic similarity coefficient varied from 0.5008 to 0.9111, with an average of 0.6618. The genetic similarity coefficient of accession 565052 from Russia and accession 598610 from Kazakhstan was the largest (0.9050), while that of accession 502576 from Russia and accession 531828 from the United States was the smallest (0.5008).

Analysis of the variation of the similarity coefficient of each plant in P. juncea material indicated that the genetic similarity coefficient between the accession 598610 individual plants from Kazakhstan had the greatest variation, which ranged from 0.5145 to 0.9111, indicating that this accession had great genetic differences and rich genetic diversity, while those plants of accession 531828 from the United States had the smallest variation, indicating that the genetic variation among individual plants of this accession had relatively small genetic differences and narrow genetic background (Table 5).

The 480 materials of P. juncea were divided into eight populations according to regions to compare the genetic diversity (Table 6). The Shannon’s information index and the change trend of genetic diversity followed the order: China > United States > Mongolia > Russia > Kazakhstan > Former Soviet Union > Estonia > Canada, the percentage of accessions with polymorphic locus from China was the largest (86.63%), and the smallest was from Canada (56.10%).

Analysis of molecular variance (AMOVA) was implemented to evaluate variance components among groups and individuals (Table 7), which showed highly significant differences (P < 0.05). Of the total accessions, 93% of the variance was due to differences within the accessions in the groups and the remaining 7% was due to differences among the groups. Of the individuals, 83% of the total genetic variance was due to differences within populations, while 5% was due to differences among populations and only 12% was ascribed to differences among individuals within the populations.

To understand the genetic structure among the eight populations of P. juncea, PCoA and cluster analysis between the populations were carried out, based on the unweighted pair-group method with arithmetic means (UPGMA) on the regional populations (Figures 3A,C). In addition, an unrooted tree was drawn (Figure 3B). The results showed that the “Estonia,” “Kazakhstan,” and “Former Soviet Union” populations were clustered together; the “United States,” “Mongolia,” “Russia,” and “China” populations were clustered together; and “Canada” was one separate group.

The N_m and F_st among the eight regions were analyzed (Table 8). It was observed that the N_m between “Canada” and the others was relatively small and the F_st was large, which could explain why the “Canada” population was separately classified into one group.

Analysis using the Mantel test showed that there was a very significant relationship between genetic distance and geographical distance (R² = 0.0269, P = 0.025). Therefore, the data showed that there was evidence of significant isolation by distance (Figure 4).

The relationship among the accessions from the different groups and subgroups based on genetic distance was further explored by UPGMA cluster analysis, PCoA analysis, and genetic structure analysis. Markers of 103 SSR primer pairs scanned from 480 P. juncea individuals were analyzed using Structure 2.3.4 software. The line graphs were established using the abscissa as the K value (K = 2–10) and the ordinate as ln P(K). The ln P(K) increased constantly with increasing K and there was no maximum value, so the best number of subpopulations could not be determined (Figure 5A). The method of Evanno et al. (2005) was used to determine the best classification number. When the delta K value changed with the K value, an obvious peak can be identified as the best classification number. Structure cluster analysis showed that the ΔK value reached a peak when K = 5 (Figure 5B), which indicated that the 480 P. juncea individuals could be divided into five different groups.

Psathyrostachys juncea is completely cross-pollinated, so great differences in genetic variation exist within natural populations. In order to reveal the genetic components between individual plants as much as possible, the materials with Q value greater than 0.5 were classified into corresponding groups, which indicated that the genetic structure was relatively unitary. The materials with Q value less than 0.5 were considered complex (Supplementary Table 5). In group 1, 96 individual plants clustered into one category, and the accessions came from China and Mongolia in Asia. In group 2, 107 individuals were from Russia. In Group 3, 54 individuals were from Canada and Mongolia. In group 4, 42 individuals originated from the United States, and in group 5, 181 individuals were mainly concentrated in many countries in Asia and Europe, and were of multi-regional origin. Meanwhile, the genetic differentiation coefficient F_st was obtained to describe the heterozygosity level of the group alleles which is often used to measure the degree of population differentiation. The average F_st of the five groups followed the order group 4 (0.5487) > group 3 (0.4587) > group 2 (0.4442) > group 1 (0.4225) > group 5 (0.3370) (Supplementary Table 6), which showed that group 4 individuals have a high level of genetic differentiation and group 5 have the lowest level of genetic differentiation (Figure 5C).

The UPGMA dendrogram also showed that the accessions could be classified into five clusters (Figure 6A). The clustering of 480 P. juncea individual materials was basically in line with the results of population genetic structure analysis. Based on the UPGMA dendrogram, some individual plants from Russia were clustered with individual plants from Mongolia and Canada, while in the population structure analysis, plant materials from Russia were classified separately into group 2. The materials from Kazakhstan, Estonia, and the Former Soviet Union were clustered into one group in the UPGMA dendrogram, indicating that the genetic distance between materials from these three regions was small. In the cluster diagram, the materials from Russia were clustered into one group, and the materials from China and Mongolia were clustered into the same group. The classification of materials from Canada and the United States was extremely complex and they were clustered into different branches.

A total of 480 individual plants of P. juncea were analyzed based on the first and second principal components (Figure 6B). The distance of their positions in the figure represents the distance of genetic relationship. Psathyrostachys juncea individual plants can be classified into five groups, and the results of the PCoA were basically consistent with those of the UPGMA dendrogram and population structure analysis. The material distribution of individual plants in the five groups was relatively centralized, indicating that they were closely related to each other. In PCoA analysis, the materials from Mongolia were classified into multiple groups. The material sources of group 2 and group 4 had a relatively independent distribution, while groups 1, 3, and 5 had multi-distributed sources.

Tillering is a major determinant of the forage yield of bunch grass, which is regulated by plant hormones, growth and development modes and environment. In recent years, new tillering genes and regulatory mechanisms of monocotyledons represented by Gramineae have been reported, which has deepened understanding of plant tillering. The regulation mechanism of plant tillering has been revealed from analysis of hormones, genes and transcription, which has provided a basis for the study of forage yield formation, environmental adaptation and the survival competitiveness of gramineous plants. Xu et al. (2020) reported that approximately 13,609 DEGs were identified from a high-tillering genotype and a low-tillering genotype of Orchard grass. In this research, P. juncea with dense tillers and sparse tillers were selected for transcriptome analysis. Among 100,560 unigenes, 9,300 genes with SSRs were identified. The significantly differentially expressed genes were screened according to the relative expression level between the two samples. A total of 2,020 significantly DEGs were identified, of which 400 were DEGs with SSR loci. In this study, we find abundant of EST-SSR exist in tiller trait related DEGs. These EST-SSR markers provide theoretical tools for the correlation analysis between tiller trait and molecular markers, and could play an important role in future breeding program to increase tiller number and forage yield. Furthermore, the P. juncea genome is more than 6.7 G base pair and limited genome sequence information was published (Li et al., 2020). SSR markers is easy to detect and analyze in a huge population with sufficient genetic variations. For genetic mapping and tiller trait QTL analysis, it is necessary to clarify the genetic distance of populations and individuals of different germplasms, especially to identify genetic differentiation between populations and individuals based on locus controlling tillering trait. It’s also necessary to analyze the genetic structure of populations to evaluate differentiation level of tiller trait related gene locus in P. juncea. The details of tillering trait transcriptome analysis of P. juncea are expected to be discussed in another manuscript.

EST-SSR is the identification of SSR by electronic screening using existing EST sequences, followed by PCR detection (Zheng et al., 2020). Using EST to develop SSR avoids the cloning and sequencing steps in the process of developing SSR primers, makes full use of existing data, and reduces the development costs (Fan M. et al., 2020). EST-SSR is well conserved and has good versatility among different species, and can distinguish materials with relatively close genetic relationships (Huang and Lu, 2010). Liu et al. (2013) pointed out in the development and genetic diversity analysis of EST-SSR primers on panicle traits of wheat that the primers had high amplification efficiency and high polymorphism information content.

The polymorphism rate among the genotypes under exploratory conditions is considered a key factor in measurement of the diversity analysis efficiency of DNA markers (Athipat et al., 2019). A large number of studies have pointed out that the polymorphism of markers affects the level of genetic diversity of plants. Generally, the genetic parameters of tested materials were more reliable using primers with high polymorphism than those of primers with poor polymorphism. Zhang et al. (2019) analyzed the population and individuals of P. juncea germplasm using SSR primers distributed on seven linkage groups of barley, and reached the conclusion that SSR markers can provide a theoretical basis for genetic diversity research on P. juncea. However, in this study, for the first time EST-SSR markers were used to analyze the population structure of P. juncea, and 103 primer pairs were selected from 400 EST-SSR primers to mark the P. juncea population. The result showed that the polymorphism of EST-SSR markers developed from P. juncea was higher than SSRs from other related species, which indicates the population structure analysis was more accurate. The genetic diversity index in this research was slightly higher than that in Zhang et al. (2019). The high genetic diversity detected in this research may be due to the use of EST-SSR molecular markers developed for the genome of P. juncea, which can better distinguish gene loci. The Shannon information index (I) ranged from 0.173 to 0.640, with an average of 0.420 in our study, which was higher than that estimated by ISSRs for the related P. huashanica (I = 0.391). This finding indicates that the selected primers can objectively reveal the genetic diversity of P. juncea germplasm resources. Overall, the 103 EST-SSR primer pairs can provide an adequate analysis of the genetic differences among P. juncea materials.

Through simple statistics analysis of seven agronomic traits related to tillering of P. juncea in samples 1 and 2 group, the variation coefficient of nutritional tiller number was the most obvious in both samples (Supplementary Table 7), which suggested a relatively obvious genetic variation of tiller number within P. juncea population. It is consistent with the results of previous germplasm evaluation studies. Fan Y. K. et al.’s (2020) find tiller number trait has the most significant variation among all traits after evaluated 30 germplasm accession of P. juncea in field trail.

The genetic similarity coefficient (GS) is an index of the degree of similarity between individuals. When the genetic similarity coefficient is greater, the genetic relationship between materials is closer. The GS of 480 individual plants of P. juncea was analyzed, and the average value of GS was 0.6618, which showed that the genetic similarity level among individual plants of P. juncea was high. The materials with GS values ranging from 0.6 to 0.7 accounted for 60% of the total, indicating that the distribution of genetic diversity among materials was relatively centralized, the range of genetic variation was small, and the genetic relationship was close due to the limitation and centralization of the natural territory of P. juncea. Liu et al. (2009) analyzed the GS of 15 P. juncea materials using the ISSR molecular marker method and found that those with a GS range from 0.6 to 0.7 accounted for 37.4%. By comparing the average GS between wild materials and cultivated varieties, it was found that the GS value of cultivated varieties (0.667) was greater than that of wild materials (0.658), indicating that the genetic relationship between cultivated materials was closer and that of wild materials was more distant. The same detection was reported in the analysis of genetic similarity coefficient of wild materials and cultivated varieties of alfalfa (Chen and Shi, 2015). The genetic similarity coefficient of 88 P. juncea materials analyzed by RAPD molecular marker showed that the smallest genetic similarity coefficient was among Kazakhstan material and the greatest was among Russian material (Wei et al., 1997). That finding differs from the result in this study, indicating that in the process of genetic evolution, plant DNA has changed and produced a variety of rich germplasm banks with change in the surrounding environment. These differences also can be attributed to differences in molecular markers and plant materials used in each study.

In this research, based on the geographical origin of P. juncea accessions, accessions were clustered into eight subgroups for genetic structure analysis. The Shannon information index of each subpopulation ranged from 0.2615 to 0.4142, indicating a high genetic diversity level. The high genetic diversity in the P. juncea population may be related to the anemophily pollination of the species. Genetic variation within the populations (83–93%) was higher than that among populations (5–7%) in the present study. This result is consistent with the findings of previous studies of P. juncea based on AFLP analysis (77% within populations) (Xiong et al., 2020), Black cottonwood (Populus deltoides) based on SSR (84.88% within individuals) (Chen et al., 2020), white clover based on AFLP analysis (84% within cultivars) (Kölliker et al., 2001) and white clover based on SSR (86.5% within cultivars) (George et al., 2006). This is also consistent with the results of studies on other cross pollinating species, such as perennial ryegrass (Bolaric et al., 2005; Van Treuren et al., 2005). Our research showed that in the process of tillering, the genetic differences of P. juncea were mainly caused by individual variation.

Classically, a high degree of gene flow, which could neutralize interspecific differentiation and intraspecific genetic drift, is extremely common in cross-pollinated plants, thereby contributing to low genetic diversity among individuals or populations (Sui et al., 2009). In our study, the P. juncea gene flow (N_m) between Canada and other countries was small, and the N_m between other countries was high. Ellstrand (1993) suggested that N_m < 1 indicates a high level of genetic differentiation among plant populations, while N_m > 1 indicates a low level of genetic differentiation among plant populations, and George et al. (2006) also found this result in cross pollinated white clover. Due to frequent gene exchange among subpopulations, the relationship between populations can be maintained, so the genetic differentiation between populations was not significant. In the process of evolution, P. juncea has undergone huge genetic variation and frequent exchange of genetic material with the materials of other populations. Therefore, germplasm resources with rich polymorphism and high genetic variation have been formed at the level of genetics and breeding. Germplasm collections and the introduction history of different countries may also affect gene flow and shape the new genetic structure of subpopulations.

The genetic structure of species is affected by the interaction of multiple factors, such as the transmission model of seeds and pollen, population demographic history, geological events, geographical or ecological barriers and divergent selection for environmental factors (Yang et al., 2017; Smith et al., 2020). Based on Nei’s GD, the eight populations were clustered into three groups. This clustering showed strong geographic regionality. China and Mongolia are also adjacent to Russia in geographical location and similar in climate environment, but the Canadian materials located in North America were clustered separately, which may be caused by the long-term accumulation of genetic variation in the process of adapting to the regional environmental conditions after the introduction of P. juncea germplasm from Eurasia to North America.

In this study, according to the results of population genetic structure analysis, the materials with Q greater than 0.5 were assembled to five groups and the individual plants from the same region were also be clustered into different groups. In the UPGMA cluster diagram, some individual plants from Russia were clustered with individual plants from Mongolia and Canada, but in the population structure analysis Russian materials were clustered independently into group 2. The clustering may be the result of breeding and domestication, which has a great impact on the diversity structure. Selection and breeding tend to keep plants with economically valuable traits such as tillering ability. In addition, different environments can also cause genetic changes and thus affect the division of population structure. It may also be that the DNA of P. juncea materials contains hereditary substances of different individual plants, and the genetic information of all individuals was integrated by the population structure. UPGMA clustering was carried out based on the GD of P. juncea materials to cluster the closely related materials into one group, which may lead to the aggregation of materials with different genetic structure. Shen et al. (2010) analyzed the population genetic structure of 64 oat germplasms and also found that germplasm from the same region were clustered into different groups.

Through the analysis of population genetic structure of P. juncea individual materials from different countries and regions, it was found that the population structure of 480 P. juncea individual plants had a certain correlation with geographical distribution. The Mantel test of genetic distance and geographical distance among materials showed that there was a significant correlation (R² = 0.0269, P = 0.025). This showed an IBD model over all sampling locations that was comparable with the results of Xiong et al.’s (2020) research, which analyzed the genetic diversity of eleven P. juncea wild germplasms using AFLP markers, and identified that the key factors inducing moderate genetic differentiation were isolation by distance (IBD). Generally, more opportunity for allelic exchange could be obtained by individuals of neighboring populations because gene flow will theoretically be obstructed by a longer geographical distance (Zhang et al., 2018). Our research showed that geographical isolation hindered gene exchange among individuals from different populations and promoted genetic differentiation. However, it was interesting that although Canada and the United States are geographically close, there were great differences in the genetic structure of these two subpopulations, and they were not clustered into the same groups in the group division in the UPGMA and PCoA analysis (Figure 6). The clusters show substantial overlap of different populations, except for those from Russia and the United States. Moreover, high values of the genetic mixture were also confirmed by structure analysis. This is mainly attributed to the cross pollination and self-incompatibility of plant species (Khan et al., 2009), human seed transplantation (Wang Z. Y. et al., 2009), different biological dispersal patterns and evolutionary forces (Chapman et al., 2010) and random dispersal in a region (Golkar and Mokhtari, 2018). In the long-term breeding process, gene exchange among germplasm resources is becoming more and more frequent, resulting in the decline of regional correlations and blurring of original geographical sources among materials, which reduces the genetic differentiation of P. juncea plants. Similar results have also been obtained for perennial ryegrass (Yu et al., 2011). In addition, compared to food crop species, P. juncea is generally used as a forage grass species, and has a relative short domestication and breeding history, and shows a certain level of consistency between genetic and geographical differentiation.

The P. juncea accessions in the present study were clustered into five groups based on PCoA, UPGMA, and STRUCTURE analyses. The classification of 480 P. juncea individual plant materials basically accorded with the analysis of population genetic structure, but there were some differences. This could be attributed to the different statistical principles applied in different methods (Lv et al., 2020; Wu et al., 2021). PCoA can provide a more valid classification based on the dissimilarity matrix of the original data, which is not strictly in line with the Hardy-Weinberg equilibrium assumption. Population structure analysis can better understand genetic diversity and estimate the variation of germplasm resources, which is conducive to their effective utilization. Structure assigns the accessions to subgroups probabilistically by a Bayesian clustering approach, and it is used for the subdivision of natural out-crossing populations. The clustering of accessions using UPGMA analysis was implemented based on genetic distance, which showed more detailed relationships among the accessions. Overall, these three methods could work together to provide a comprehensive understanding of the P. juncea population genetic structure. Moreover, the differences in P. juncea revealed can be assessed for their high breeding and hybridization advantages when selecting parent combinations, and provides one basis for the selection of P. juncea to make cross combinations with strong genetic differences.