Genetic Variability and Aggressiveness of Tilletia indica Isolates Causing Karnal Bunt in Wheat

Karnal bunt caused by Tilletia indica is a quarantine disease of wheat causing huge economic losses due to the ban on the import of bunted grains. This study was designed to characterize pathogenicity, aggressiveness and genetic diversity of 68 Tilletia indica isolates collected from different geographic regions of Pakistan. Forty-six isolates were tested for their pathogenicity on eight wheat varieties, out of which three were non-aggressive. The coefficient of infection (CI) ranged from 15.73% (PB-25) to 10% (PB-68, PB-60, and PB-43). The isolates collected from central Punjab showed higher infestation compared to other isolates. Among the wheat varieties used for the aggressiveness study, WL-711 showed susceptible reaction with 10.88% CI, while NIFA-Barsat, HD-29, Janbaz, Bakhtawar-92, Tatara, and AARI 2011 showed resistance to the highly resistant response. These isolates were amplified using 31 random amplified polymorphic DNA (RAPD) markers and 32 inter-simple sequence repeat (ISSR) markers for diversity analysis. The principal component analysis (PCA) and analysis of molecular variance (AMOVA) showed greater divergence among isolates collected from Punjab and Khyber Pakhtunkhwa (KPK), with a moderate level of admixture. The isolates from Faisalabad (Punjab) were more aggressive compared to isolates from KPK and were clearly separated based on PCA, indicating the significant genetic distance in the populations. Our findings will assist breeders and pathologists in better understanding the pathogenic variability in Tilletia indica and in subsequent disease management.


Introduction
The pathogen, Tilletia indica (T. indica), causes Karnal bunt (KB) in wheat, durum wheat, and triticale, and is also named as a partial bunt. It is an internationally quarantined disease [1] affecting quality and production of wheat, and there is zero tolerance in various countries against KB for importing bunted wheat grains [2]. The first report of occurrence was from the Indian city Karnal [3] and since then it has been reported from other countries, namely Afghanistan, Iraq, Nepal, Pakistan, Iran, Mexico, Brazil, the United States (New Mexico, Arizona, Texas, and California) and, finally, from South Africa (Northern Cape Province) in the present century [4].
Management of KB has become a serious issue due to various factors such as the pathogen's dispersal mode, lack of tolerant wheat varieties, and long-term survival of T. indica spores in the soil. Therefore, cultural and fungicide control of KB is quite challenging. A compatible coupling of secondary sporidia increases the chance of variability because of heterozygosity. The higher crossover frequency in the pathogen suggested that intraspecific and interspecific hybridization may promote recombination and greater genetic diversity [5]. Variability studies between fungal isolates and virulence are important for disease management and resistance improvement programs. Genetic variation of a pathogen can be understood through the use of genetic markers. Random amplified polymorphic DNA (RAPD), ISSRs, RFLP and AFLP markers and have been confirmed to be useful in finding out the genetic and phylogenetic relationships of organisms [6][7][8]. However, ISSRs and RAPD produce more bands and represent more loci than RFLP and AFLP; and previously 34 RAPD primers and 28 inter-simple sequence repeat (ISSR) markers on 10 isolates collected from different regions of India [8,9]. The mystery of genetic variation in pathogens is widely understood through the use of genetic markers. Diversity of T. indica isolates was analyzed using pathological and molecular markers [10,11]. Ref. [12] used ISSR-18, ISSR-1, ISSR-17, ISSR-1, ISSR-19, and RAPD 16 primers on 8 isolates of T. indica. Effect of host determinant on variability was determined. They reported both marker systems have indicated the genetic variation between treated and untreated samples of monosporidial strains.
Despite its economic significance, very little work has been done on the molecular component and pathogenesis of the pathogen. In this study, attempts were made to analyze the genetic variations of T. indica isolates collected from different districts of Punjab and Khyber Pakhtunkhwa (KPK). More importance was given categorically for precise detection of variation among T. indica isolates by using RAPD-ISSR markers. This study, to our knowledge and review the first time a large number of Pakistani isolates has been studied by using molecular techniques of ISSR and RAPD, will explore the links to find out the basis of pathogenicity in T. indica.

a.
Sample collection Infected wheat grains with KB were randomly collected from Punjab and Khyber Pakhtunkhwa (KP). Surveys were carried out to find the current status of KB during 2013-2014, and 48 isolates were collected [13] (Table 1). Twenty isolates collected during 2010-2011, including Pb-25 (aggressive isolate) [14], were also used for the molecular analysis. The geographic coordinates of the sample collection sites and host cultivars are also given in Table 1. The map positions of geographic coordinates are given in Figure 1.

b. Isolation and multiplication of T. indica isolates and fungal inoculation
The infected seeds of 68 isolates (Table 1) containing teliospores were suspended in 10 mL of sterilized distilled water and vortexed for 2 min. The spore suspension was passed through a 60 µm sieve and centrifuged at 12,000 rpm for 2 min. Pellet was suspened in 1% sodium hypochlorite solution and centrifuged for 30 s at 12,000 rpm. Washing of spores with sterilized distilled water was done thrice and for this the spores were first centrifuged at 12,000 rpm for 2 min and then the pellet was re-suspended in sterilized distilled water and plated on the water agar with 1000 µL micro-pipette (Gilson). The plates were incubated at 20 ± 2 • C for 10 days. After 10 days of incubation, the plates were examined under stereomicroscope for the confirmation and morphological studies of the pathogen as described previously. Colonies from single teliospores of T. indica were cultured on PDA medium. After isolation of the pure cultures of each isolate of T. indica, their mass multiplication for pathogenicity study was done. A small disk of freshly cultured T. indica isolate on water agar was cut and placed on the upper edge of flasks/plate containing PDA for further multiplication and inoculum preparation. These flasks were labeled and incubated at 20 ± 2 • C for 15 days. The flasks were stored in the refrigerator to use for further experimentation [15].

c.
Pathogen (T. indica) material and DNA extraction As described above in Section b, the mycelia of 68 isolates were collected for DNA extraction, from 15-20-day-old cultures that were incubated at 20 • C. Each mycelial layer of T. indica was filtered through autoclaved Whatman filter Paper No. 1 and washed thrice with sterile distilled water. The hyphae mat was dried between sheets of paper at room temperature. The dried mycelia were stored at −20 • C for further use. All isolates were cultured separately, multiplied and maintained on Potato Dextrose Agar. The sporidial suspension of each isolate was prepared by adding sterilized distilled water, the fungus was scratched from PDA [16] and the concentration of inoculum was adjusted to 50,000 spores mL −1 with a hemocytometer [17].

Pathogenicity Test of KB Isolates
Eight wheat varieties were selected as these were mega commercially grown varieties in the country, to assess the pathogenicity of isolates, including six commercial wheat varieties, namely, AARI-2011, Bakhtawar-92, Fakhre-Sarhad, Janbaz, NIFA-Barsat, and Tatara, along with 2HD-29 and WL-711 as resistant and susceptible check, respectively. Seeds of wheat varieties were obtained from National Coordinated Wheat Program, National Agriculture Research Centre Islamabad. The wheat varieties were sown in pots, during November 2014. After germination, three plants were maintained per pot. Five heads per pot were inoculated and three inoculated spikes were selected to assess the disease incidence and severity. The experiment was conducted in three replications in a completely randomized design in a two factorials setup. a.

Boot inoculation
The plants were inoculated individually with a spore suspension of each isolate using the boot inoculation method at booting stage [18,19]. Varieties were sown in pots and each variety was replicated thrice in different pots in a glasshouse. Three spikes per variety of each pot were inoculated just before the emergence of awns. Suspension of T. indica spores (1mL) from 10-15-day-old cultures was injected into the boot with the aid of a hypodermic syringe. After inoculation, the spikes were labeled and covered with glycine bags to maintain the maximum moisture content for fungus proliferation [20]. The inoculated pots were kept under high humidity (approx. 80%) in the glasshouse, an automated misting system was equipped, and the sprinklers were sprayed five times per day for 20 min each time, under a temperature regime of 18-20 • C [21]. b.
Harvesting and scoring The inoculated spikes were harvested at maturity (58-60 days) after inoculation and threshed manually. The infected seeds were divided into different degrees/grades of infection depending on the extent of the damage on the scoring scale 0-5 [22] (Table 2). The coefficient of infection (CI) was calculated [23] as follows: Coefficient of infection = Gross Total/Total numbers of seeds × 100 All isolates of the fungus T. indica were selected based on the CI. The isolates were classified as the most aggressive (20. Disease severity of the isolates and varietal response were analyzed based on their CI using Statistix 8.1. Average, percent seed infestation and LSD of both CI and PI were calculated using Statistix software. Design of the experiment was a completely randomized design in two factorials. Analysis of variance as appropriate for the completely randomized design was applied to assess the effect of isolates and varieties on the coefficient of infection and percent seed infestation. The boxplots were prepared in the R-Statistical software.

Genetic Diversity Analysis of KB Isolates
a.
DNA extraction The dried mycelium of T. indica (200 mg) of each isolate was taken for DNA extraction using the CTAB method [26] and was allowed to run on 1% agarose gel, along with standard marker for its quantification. The concentration of amplified DNA was assessed visually by checking the band intensity in comparison with Lambda (λ) DNA of known concentration under the UV transilluminator using a gel documentation system. Sharp and separate bands were estimated as good, whereas improper dissolution and thin bands were estimated as sheared [27].

b. Primers used for molecular analysis of T. indica
A set of 31 RAPD (Table S1) and 32 ISSR (Table S2) were used for this study. Primers were diluted up to 10 pmol/µL. PCR amplification of T. indica isolates was performed in a volume of 20 µL reaction on automated Applied Biosystems Thermal Cycler (Verity 96 well), with composition of 2 µL of 10X buffer, 2 µM of each primer, 0.4 uM dNTPs mixture, 0.2-unit Taq polymerase, 2.4 mM MgCl2, 1ul of template DNA and 12 µL of ddH2O.
The ISSR amplification regimes were performed at 94 • C for 3 min after that, 35 cycles of 94 • C for 40 s, later annealing was done according to primer for 40 s, followed by extension at 72 • C for 1 min, and a final extension at 72 • C for 10 min. For RAPD amplification, PCR cycles consisted of DNA denaturation at 94 • C for 3 min followed by 37 cycles of 92 • C for 1 min, then annealing was done according to primer for 1 min then extension at 72 • C for 2 min and followed by a final extension at 72 • C for 15 min.
The amplified DNA or products were run on 1.8% agarose gels at 100 V in 1 × TAE buffer for 1 h along with 1 kb DNA ladder. The DNA band pattern was visualized with UV light and recorded by the imaging system of gel documentation. PCR amplification was repeated thrice to reproduce the DNA profiles with all of the selected primers. A negative control (without DNA template) was added in all of the PCR reactions. c.
Data analysis Data were computed to evaluate the various parameters for the comparison of RAPD and ISSR primers. The gel of all sets of primers was analyzed. Presence of bands was scored as 1 and absence as 0. A dendrogram was constructed using UPGMA (an unweighted pair method with arithmetic mean) to group individuals into discrete groups [28]. Polymorphic information content (PIC) was calculated as follows: PIC = 1 − (p/q) 2 where p is total alleles detected at a given marker locus, q is a collection of genotypes studied. The multiplex ratio (MR) for each primer was assessed by dividing the total number of amplified bands (mono and polymorphic) amplified by the total number of primers (the combination of primers used -n). By using a marker, the average DNA fragments amplified by genotype are known to be a multiplex ratio (n). The number of polymorphic loci within the set of interest of the germplasm analyzed by experiments is called effective multiplex ratio (E) and estimated as E = nβ, where β is a fraction of polymorphic markers. It was calculated by considering the polymorphic loci (np) and a non-polymorphic locus (nnp), as ß = np/ (np + nnp). The usefulness of a given marker system is the balance between the level of polymorphism identified and the extent to which the marker identifies multiple polymorphisms. An informational product measured by PIC and an effective multiplex ratio called the Marker Index (MI) provided a suitable evaluation of the utility of the marker. MI = PIC × E or MI = n × β × PIC.
AMOVA was performed using software R to analyze the genetic diversity of T. indica isolates. Escoffier et al. [29] estimated the variance components of the RAPD and ISSR profiles and the assessment of intra-and inter-genetic diversity parameters such as the diversity of the gene Nei (Hexp), Standard error for the rarefaction analysis (SE), Shannon-Wiener Diversity index(H), Stoddard and Taylor's Index (G), Simpson's index (Lambda), (E.5)Evenness, Association Index (Ia), Standardized Index of Association for each population factor (rbarD), Observed number of alleles (obs), Standard deviation of data (Std. obs). Genetic structure was studied using the diversity statistics of Nei gene, diversity within population (HS), total genetic diversity (HT), population diversity (D-het) and gene differentiation coefficient (GST) calculated based on the Nei method by using software R [30]. Gene flow estimated from GST was calculated as follows: Nm = 0.5 (1-GST)/GST. Principal component analysis (PCA) was performed to obtain a 3D image graphically representing the genetic diversity between the isolates. It was generated using R software [31].

Variability in Pathogenicity for Isolates of T. indica
Highly significant variability was observed among varieties and isolates for both CI of and PI and their interaction. A total of 49 isolates of T. indica (Table 3), isolated from infected wheat seed samples collected from different agro-ecological regions of wheat in 2013-2014, showed variable virulence on different varieties. Forty-six isolates were found to be pathogenic to all varieties of wheat, while three isolates were found non-virulent ( Table 3). The disease reaction of forty-nine isolates of T. indica on these wheat varieties is shown in Table S3. The highest CI (15.73%) was recorded in isolate PB-25 followed by PB-55 (11.37%) from Faisalabad, while the lowest severity was observed in PB-68, PB-60, and PB-43(10%) collected from Jhelum, Bahawalnagar, and Bahawalpur ( Table 3). The mean value of CI showed that isolate PB-25 was aggressive on WL-711 and Fakhre-Sarhad compared to other varieties (Table S3).
The aggressive reaction of forty-nine isolates of T. indica on these wheat varieties was subjected to analysis of variance (ANOVA). The results revealed that there was a highly significant difference of isolates among themselves and their reaction on different varieties tested at p ≤ 0.01. To clarify the aggressiveness behavior of forty-nine isolates of T. indica based on their reaction on eight wheat varieties, the results on these varieties were subjected to cluster analysis using Minitab16. The dendrogram ( Figure 2) is based on the combined means of the coefficient of infection of T. indica on varieties. Based on aggressiveness behavior, the isolates were clustered into different groups and subgroups, namely Group I (subgroup = 0), Group II (subgroup = IIA, IIB), Group III (subgroup IIIA and IIIB). Group I comprised one isolate (PB-25) and belonged to a moderately aggressive category. Group II consisted of slightly aggressive isolates having a single isolate in each subgroup IIA, IIB. Group IIIA comprised three non-aggressive isolates, while Group III B comprised forty-three isolates and belonged to weakly aggressive category. Although the two groups (Group II, subgroups II A and II B, and Group III subgroup III B) were based on CI, they belonged to a weakly aggressive category, but based on cluster analysis and similarity these comprised different groups (Figure 2).

Reactions of Wheat Varieties against KB Isolates
Eight varieties, including one resistant and susceptible check, were tested for their responses against 49 isolates of T. indica. WL-711 exhibited moderately susceptible response and showed a maximum CI value (10.88%) ( Table 4). The varieties NIFA-Barasat, HD-29, Janbaz, Bakhtawar-92, Tatara, and AARI 2011 showed a resistant to highly resistant response to all the isolates. Fakhre-Sarhad showed a moderately susceptible to highly resistant response (Table S3). Data for mean percent seed infected with KB disease for each treatment was calculated. The data were subjected to ANOVA. The results revealed that there was a highly significant effect of isolates among themselves and their reaction on different varieties tested at p ≤ 0.01. The isolate PB-25 (32.93%) caused the highest PI on all the wheat varieties tested followed by PB-55 (17.99%), whereas the isolate PB-64 (10.29%) exhibited the least infection. The three isolates PB-43, 68, and 60 were found to be nonpathogenic on these varieties ( Table 3). Comparison of the overall distribution of CI values and PI revealed a strong correlation between the two parameters, with an R 2 value of 0.764.

Variability in T. indica Percent Infection as a Function of Sampling Location
Analyzing the response of T. indica isolates, while considering their location of sampling, revealed substantial variability in their pathogenicity as assessed by the coefficient of infection and percent infection (Figure 3). The result revealed that based on the sampling from different geological areas the variability significantly exists in the isolates. Isolates sampled from Faisalabad had the maximum variability for percent infection of seeds, followed by samples collected from Rawalpindi, Charsada, Chakwal and Gujranwala. Across the studied major wheat growing regions of Pakistan, the samples collected from southern Punjab had the maximum percent infection (with a mean value of 14.13) followed by Khyber Pakhtunkhwa (Table 4). A similar pattern was observed for CI, the maximum observed at district Faisalabad of Punjab (13.7), while the minimum (10.0) was observed for isolates sampled from Northern Punjab district Jehlum and Rawalpindi, and southern Punjab district of Bahawalnagar (Figure 3).

Banding Patterns of RAPD and ISSR Primers
Thirty-one RAPD primers (Table S1) were selected for this study based on their reproducibility, the number of polymorphic fragments per primer and the level of polymorphism detected in a particular population. Eight primers were selected. A total of 99 different bands were amplified and found to be (100%) polymorphic. Bands amplified per primer ranged from 9 (OPA-9) to 15 (OPA-20, RAPD-16) with an average of 12.37. The size of the product ranged between 0.25 kb and 2.5 kb. Four unique bands were found: OPA 3 (on isolate 17), OPA 18 (isolate 61), and OPA 20 (isolates 15 and 30). The average PIC values, is a description of allele diversity and frequency, found in this study among the T. indica species were 0.93, ranging from 0.8 (RAPD 16) to 0.98 (OPA 20) (Table 5).  Among thirty-two ISSR primers (Table S2) tested, 20 were useful for characterizing samples, while 11 were excluded due to the absence of amplification. The 20 selected ISSR primers produced 215 bands, 100% of which were polymorphic with an average of 10.75 bands per primer. Bands produced by each of the twenty ISSR ranged from 6 (ISSR-1) to 15 (876 and 856-2). The size of the product ranged from 0.25-3.0 kb. Nine unique bands were produced from the primers 812 (15,58), ISSR-1 (41), 845 (38,44), 847 (46), 817 (15), 856 (1), and 856-2 (17). Values in parentheses denote the code number assigned to the isolates as sequenced in Table 6. The average PIC values of the T. indica species used in this study 0.902 ranged from 0.69 (ISSR 17) to 0.973 (ISSR-18) ( Table 6). Values in parentheses of unique bands denote number of codes assigned to the isolates given in Table S3 AMOVA analysis among and within isolates of T. indica under study was performed to test the population differentiation utilizing molecular markers. For the analysis of molecular variance based on RAPD, 7.76% diversity was computed among the population, whereas variations within the isolates were 90.44%. ISSR analysis showed 7.408% variability among the population and 91.09% within the isolates. Genetic diversity between the isolates and among the population of T. indica was statistically significant (Table 7).

Principal Component Analysis (PCA) of RAPD and ISSR
The principal component analysis was done to verify the grouping of isolates based on RAPD and ISSR. The first three Eigenvectors showed the minimum range of the data and accounted for 30% of the total variance in RAPD, 23.23% in ISSR, whereas the first component accounted 14.89% and 11.77% variance, respectively (Table S4). Two-dimensional (2D) distributions of the T. indica isolates based on PCA values in RAPD and ISSR analysis gave a perfect discrepancy of the sixty-eight T. indica isolates used in the study ( Figure 4A). PCA plot of RAPD and ISSR showed that the groups (Punjab and KPK) are intersecting but are different ( Figure 4B).

Combined RAPD and ISSR Cluster Analysis
Euclidean similarity coefficient was calculated using the combined 8 RAPD and 20 ISSR primers data. Similarity coefficient (Sm) values ranged from 0.4 to 0.9 in combined RAPD and ISSR analysis. UPGMA cluster analysis grouped the isolates of T. indica into three major groups (I, II and III) that were further subdivided into two subgroups each. Group, I A comprised one isolate, while I B comprised 9 isolates. Group II A comprised two isolates, whereas Group II B comprised 22 isolates. Group III A comprised thirteen isolates and Group III B comprised 21 isolates ( Figure 5).

Genetic Diversity in T. indica Isolates
Genetic diversity of T. indica isolates based on RAPD and ISSR was calculated. In the analysis of RAPD and ISSR, isolates from two provinces (Punjab, KPK) of Pakistan were evaluated. Among 68 isolates, 55 isolates from Punjab and 13 from KPK were evaluated in this study and all were multilocus (MLG). Expected MLG based on rarefaction was 13 for both RAPD and ISSR analysis ( Table 8). The minimum spanning tree for the frequency of isolates for RAPD and ISSR markers using Bruvo's distance was formed ( Figure 6A,B). The Shannon-wiener index (H) and MLG diversity (G) in RAPD were 4.01 and 55 for Punjab, whereas in the case of ISSR it was 2.56 and 13 in KPK population, respectively. Simpson's index (lambda) was estimated (0.98) and (0.923) in Punjab and KPK, respectively, and found to be same for all sets of markers. Evenness was 1 and showed that genotypes were fairly evenly distributed. Evenness of Punjab and KPK isolates (1) was similar in RAPD, ISSR. The Nei's gene diversity ranged from 0.2 to 0.209 in Punjab and 0.204 to 0.242 in KPK, and was detected by the RAPD, ISSR, (0.2, 0.209) respectively. Test for multilocus association (IA) was significantly different from zero (p < 0.001), indicating a linkage disequilibrium in all populations (Table 8). Ia values for RAPD (5.32) and ISSR (5.59) in the isolates collected from Punjab and association of KPK samples were 5.2 for RAPD and 3.83 for ISSR. Standardized index of association ranged from 0.0179 to 0.0683 and showed the maximum value in RAPD (0.0374) as compared to ISSR (0.0179 in Punjab, and showed the maximum value in RAPD (0.0452) as compared to ISSR (0.0149) in KPK (Table 8).

Genetic Structure and Gene Flow in T. indica Isolates
Total population diversity (HT), within population diversity (Hs), diversity between the population (D-het) and Gene flow (Nm) were calculated by using R software. HT (total diversity for among the population) was higher in ISSR (0.23), while the lowest diversity was found in RAPD (0.21). The highest value of HS was found in ISSR 0.22 and (0.19) in RAPD. The mean coefficient of gene differentiation (GST-est and G prime-st) was higher in RAPD (0.05, 0.13) than ISSR (0.05, 0.12), respectively. Gene flow was relatively higher in ISSR (9.69). The diversity between the populations (D-het) was higher in ISSR (0.029) (Table S5).

Comparison of RAPD and ISSR in Evaluating Genetic Diversity of T. indica
The marker systems (RAPD and ISSR) were compared according to different criteria ( Table 9). Detection of polymorphism, unique band, means PIC value, MR, EMR, and MI. Analysis by RAPD and ISSR showed a polymorphism level of 100%. In the present study, the average number of DNA fragments amplified by RAPD and ISSR was 12.37 and 10.75, respectively. This high variation produced in the number of fragments by these arbitrary primers can be ascribed to differences in the binding sites along the genome of the T. indica isolates used. The PIC value of the RAPD primers was higher (0.92). The multiple ratio (MR), effective multiplex ratio (EMR) and marker index (MI) in RAPD (Table 9) were higher than ISSR. The RAPD and ISSR data were combined for UPGMA cluster analysis. The combined cluster of RAPD and ISSR was 80-85%, similar to each of the separate (RAPD or/ISSR) cluster analyses ( Figure 6).

Discussion
Compared to other sumt fungi, the control of KB is challenging due to its heterothallic nature and sporadic occurrence [22]. Unfortunately, no work was previously conducted to understand and explore the pathogenicity and genetic variability of T. indica in Pakistan. However, limited studies on the variability of T. indica in terms of size of primary and secondary sporida, germination percent, morphology [13,15,32] host reaction, and aggressiveness among isolates in pathogenicity tests were reported earlier from India and Pakistan [14,33,34].
Our results revealed significantly variable pathogenicity in T. indica isolates collected from different locations of Pakistan in terms of CI and PI (%). The isolate PB-25 reported earlier was also found significantly pathogenic as compared to other isolates assessed in the current study ( Table 2). Among other isolates, PB-55 has an aggressive response and was collected from the same district as the reference isolate [14].
Based on pathogenic behavior, the isolates were arranged into three major groups. Group I contain moderately aggressive isolate (PB-25), Group II contained weakly aggressive isolates, and Group III comprised weakly and non-aggressive isolates ( Figure 2). Grouping of isolates into various clusters was based on their response to host [16,35]. The clustering was independent of the location of sampling for the given isolate. This would imply that variable pathogenicity exists among the isolates prevalent at a given location in major wheat-growing regions of Pakistan. Previously, similar results of lack of location-specific virulence variability were reported for Puccinia striiformis populations from Pakistan [36]. The prevalence of aggressive isolates has been suggested to be present in southern Punjab, Pakistan [37]. However, our findings revealed the prevalence of the aggressive isolate from Northern Punjab (Table 3), which is an indication of shifting of the pathogen from one area to another area. This has proven to be an alarming situation for the country as the area belongs to the major wheat-growing zones of the country.
Considering the airborne and seed-borne dispersal of the disease [38], the hot spot of KB could be shifting from one place to another over the seasons. This will require a more collaborative effort while considering the deployment of varieties and disease management strategies. The situation could further be complicated by the assortment of virulence genes during sexual reproduction of the pathogen, which leads to the breakdown of host races [10,39] since the life cycle of T. Indica differs from those of other bunt diseases of wheat in several ways as they germinate and produce haploid primary sporidia (basidiospores) that conjugate immediately on the promycelium and then produce dikaryotic infectious hyphae or dikaryotic secondary sporidia that systemically infect seedling plants. In contrast, teliospores of T. indica germinate and produce monokaryotic, haploid, and primary sporidia that do not conjugate, but rather germinate and produce numerous monokaryotic, haploid, and secondary sporidia. The secondary sporidia are airborne to plant surfaces where they may germinate and produce additional generations of secondary sporidia [2,8,33]. Further, because of the heterothallism and the complex phenomenon of race existence, the varieties needed to be tested against multiple isolates for their pathogenicity [25,40,41].
The reaction of the eight varieties against 49 isolates of T. indica revealed WL-711 exhibited a moderately susceptible response. The varieties NIFA-Barasat and HD-29, Janbaz, Bakhtawar-92, Tatara, and AARI -2011 showed a resistant to highly resistant response against all the isolates. The variety Fakhre-Sarhad showed a moderately susceptible to highly resistant response (Table S3). These findings show that these varieties may possess multiple disease-resistance genes [14,42] and reported a resistant to a high-resistant reaction in varieties such as NIFA-Barasat, Punjab-2011, BARA-09 and Seher against a set of 21 T. indica isolates. The significant effect of isolates and varieties suggested the presence of pathogenic variation in T. indica isolates as well as the resistance genes in these cultivars [34], which could result from the gene-to-gene interaction dependent on a strain-host relationship [37].
There are several field experimental studies to assess the response of varieties to natural infection of KB; however, limited studies have been done to assess the variability in a large number of pathogen isolates in terms of aggressiveness. Although such studies are regularly done for rust pathogens, fewer report on KB. The pathogenic variability observed in many isolates is helpful in disease-resistance development and deployment [43,44]. The pathogen adaptation potential is directly linked with the pathogen isolate diversity, which enables the acquisition of virulence against the deployed host resistance [45]. The study must be conducted across multiple years to track the pathogen evolution in time. The isolates must also be compared with other regional partners in India, Afghanistan, and Iran to understand the potential of the pathogen in the context of invasions, as revealed in another wheat pathogen such as rusts. Such studies on pathogenic variability must thus enable better development and deployment of KB tolerance and resistance at the farmer field.
T. indica isolates were found to be highly diverse even in single infected grain [46]. Previously, the variability of T. indica was analyzed by using RAPD and ISSR on a very limited number of isolates. RAPD and ISSR techniques generated numerous polymorphic bands, and the presence of a high percentage of polymorphism indicated high variability in this pathogen targeted by the primers [9]; however, in our study RAPD and ISSR primers generated 100 percent polymorphic bands, and no monomorphic band was found in each set of primers. Similar results were observed by other scientists who reported that DNA markers such as RAPD and ISSR [9,47,48] are useful for detecting, differentiating, and determining phylogenetic relationships between isolates of morphological species of different pathogens [7,8,49].
In this study, some unique bands of RAPD (4) and ISSR (9) were obtained. These unique bands were special to a particular isolate, making them distinguishable from others. These unique bands can be used as markers to identify the important isolates. Parveen et al. [9] also observed a similar pattern in banding profiles obtained from RAPD and ISSR markers. In our study results of AMOVA based on ISSR and RAPD marker analyses showed the genetic variation of T. indica isolates under study. Significant variability was found in gene diversity among the population and within the isolates at p-value 0.05. Although the difference between the isolates was not significant, this was due to the small geographical distance of the sample collection sites for the isolated pathogen. ISSR-based AMOVA reveals the highest differences of T. indica isolates when comparing variations within species. However, in another study conducted by Medhi et al. [50] on zanthoxylum spp., they reported high variation in RAPD markers in AMOVA analysis.
Standardized index of association of Punjab and KPK showed the maximum value in RAPD (0.034 and 0.045), respectively, in the whole population. The tests for multilocus association (IA) were all significantly different from zero (p < 0.001), revealing that there is linkage disequilibrium in all populations. Genetic diversity within a species is often associated with geographic extent, reproductive patterns, coupling systems, and seed dispersion and reproduction [51]. Genetic structure and gene flow of the T. indica isolates were found by the mean of the parameters such as total and within population diversity, coefficient of gene differentiation, and gene flow. Total (HT) and within population (HS) diversity were found to be the highest in ISSR with the high gene flow compared to RAPD. Coefficient of gene differentiation was higher in RAPD. In population genetics, the value of gene flow (Nm) ≤ 1.0 (less than one migrator per population in the population) or equivalently, the gene differentiation value (GST) > 0.25 is generally considered a threshold [52]. Genetic flow values above 1 are "strong enough" to avoid significant differences due to genetic drift. The extension of genetic variation in species always favors selection (drift, gene flow) and nonselectivity (natural selection) in population segmentation. Taking into account these criteria, the T. indica population have exceeded the threshold level, and ISSR is excellent in gene flow and gene differentiation assessment (ISSR (Nm = 9.6, GST = 0.05) in comparison to RAPD (Nm = 8.84, GST = 0.05). This information generated on such a large number of isolates of T.indica from diversified ecological zones of Pakistan to assess genetic diversity is a first attempt that could be helpful for not only the breeders but the pathologists of South East Asia in combating this quarantine significance disease by exploring appropriate links for resistance genes.

Conclusions
Knowledge of the aggressiveness, pathogenicity, and diversity of the T. indica is essential for effective management of the KB disease. The pathogenicity of a countrywide collection of KB isolate was resolved, which indicated more aggressive isolates from Punjab compared to the KPK part of Pakistan. Both the molecular marker systems, RAPD and ISSR, further supported the pathogenicity by similar patterns of genetic diversity. The molecular markers identified the genetic diversity and differentiated the isolates based on geographic prevalence and host cultivar. The isolates from Punjab were more genetically diverse compared to the isolates from KPK.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jof8030219/s1, Table S1, List of RAPID Primers used for the molecular characterization of T. indica; Table S2, List of ISSR Primers used for the molecular characterization of T. indica; Table S3, Disease reactions of eight wheat varieties to 49 isolates of T. indica; Table  S4, Eigen value, explained variance and cumulative variance in the PCA using characters to classify 68 isolates by RAPD and ISSR; Table S5, Genetic structure of Tilletia indica isolates and estimates of gene flow within the genus.