A diversity panel of 200 accessions of Brassica and related species including B. carinata (83), B. rapa (90), one accession each of B. barrelieri, B. deflexa, B. juncea, B. maurorum, B. oxyrrhina, B. ruvo, B. tournefortii, E. sativa, M. longipetala, S. alba, and S. erysimoides, two accessions each of A. thaliana, B. nigra, B. napus, and S. arvensis and eight accessions of B. oleracea were obtained from the Australian Grains Genebank, Horsham (Raman R. et al., 2014). In addition, the F₂ population comprising 300 individuals was developed from a single F₁ cross between BC73526 (shatter resistant with high RE) and BC73524 (shattering prone with low RE) to identify the QTL associated with pod shatter resistance. Both parental lines were selected on the basis of their contrasting rupture energy values among 83 accessions of B. carinata. Each F₂ line was selfed to generate F_2:3 population to confirm phenotypes.

The diversity panel comprising 200 accessions was grown in white plastic pots (10 inch diameter, Garden Plastic city, Australia) in 2012 and 2013 at the Wagga Wagga Agricultural Institute, New South Wales, Australia. Both trials consisted of a 4 range by 100 row array with two replications. Five plants were grown per pot. Passport data on days to first flowering (first open flower on at least two plants in a pot) were recorded. At maturity, five pods from each plant were collected to evaluate for shatter resistance using the pendulum test as described previously (Raman H. et al., 2014). Pod length from each test sample was measured with a scale excluding the length of ‘beak’ to adjust the position of the pod when pendulum strikes. In the present study, we only focused on 83 B. carinata accessions for pod shatter resistance (Table 1).

The two parental lines and their F₂ population of 300 plants were grown in 2015 in white plastic pots (10 inch diameter, Garden Plastic city, Australia) under birdcage conditions at the Wagga Wagga Agricultural Institute, New South Wales, Australia. Plants were watered daily, fertilized weekly using in-line liquid fertilizers, and protected from aphids. A total of 71 F₂ plants showed abnormal phenotypes with flower sterility; these individuals were discarded from genetic analysis. Five pods from 229 F₂ plants (normal phenotype) were collected in the 50 mL tubes containing a silica sachet for further testing of pod rupture energy. Days to flowering was recorded daily for each F₂ plant. All 229 F₂ plants were enclosed with pollination bags to get pure F₃ progenies, while leaving the primary stem out for the natural pod development for shatter testing. Ten F₃ plants from 229 F₂ families were grown in 2016 in a 20 row × 12 column array design including nine controls and two parents at Wagga Wagga. Five pods were collected per F₃ plant. For validation, 58 F_2:3 families (29 high RE and 29 low RE) and parents were tested with pendulum test as described earlier (Raman H. et al., 2014).

Anatomical features of pod DZ were observed in 30 random F₂ plants and five F_2:3 progenies from 20 F₂ plants selected on the basis of their rupture energy (10 each with low RE and high RE). Pods were collected at 35–40 days after anthesis. Hand sections were prepared from one cm from the pedicel end of the pod. Fresh sections were observed for autofluorescence using a fluorescence microscope at the Charles Sturt University, Wagga Wagga. Photographs were taken using a Zeiss Axiphot microscope fitted with a Sony Cyber-shot digital camera.

The rupture energy data of an F₂ population and of a set of 83 diversity lines were square-root transformed to normalize and further analyzed using ASREML in R. Genotype was considered as a fixed effect and environment as random effects. The estimated means for each genotype were used for further genome-wide association analysis. The correlation between rupture energy in 2012 and 2013 was calculated using Pearson’s correlation coefficient. RE of five pods of each F₂ plant was averaged and used for QTL analysis.

Young leaf tissue of the field grown plants was collected for DNA isolation. Tissue were ground in liquid nitrogen and extracted using a method described in Raman et al. (2005). The diversity panel of 83 B. carinata accessions and the F₂ population comprising 300 lines were genotyped with the genotyping-by-sequencing based DArTseq marker approach (Raman H. et al., 2014) at the DArT P/L, University of Canberra, Australia.

In order to determine molecular diversity in B. carinata, we genotyped 83 accessions with high-quality DArTseq markers having call rate of ≥90%, ≤5% of missing data and minor allele frequency (MAF) of >0.05. Hierarchical cluster analysis based on Euclidian distance was conducted using the software, PRIMER 6 (Clarke and Gorley, 2006). Principal coordinate analysis was performed to understand the global diversity among accessions. Bayesian clustering was performed to infer the number of sub-populations among 83 accessions using the software package STRUCTURE v 2.3.4 (Pritchard et al., 2000). The program was run using admixture model with correlated allele frequencies. The presumed sub-population number (k) was set from 1 to 5. Ten runs for each k were performed with 20,000 burn in period and 50,000 Markov Chain Monte Carlo iterations per run, with no prior information on the origin of individuals. The best k value was determined by using the (i) logarithm likelihood for each k [L(k)], (ii) an ad hoc quantity (Δk) according to Rosenberg et al. (2001) and Δk method described by Evanno et al. (2005), respectively. Genotypes were classified into subpopulations based on their membership coefficients estimated in STRUCTURE.

The linkage map of F₂ population was constructed using DArT P/L’s OCD MAPPING program (Petroli et al., 2012). Markers were clustered into linkage groups according to the method described by Wu et al. (2008). Markers with identical genotypes are placed in redundant bins, and the resulting markers/bins within each linkage group were ordered using the traveling salesman path solver program Concorde (Applegate et al., 2006). The linkage map was constructed for each parent by combining the relevant in silico DArT and SNP markers. A linkage map was chosen to be seed map and then a consensus map was constructed using the markers in common for the complete F₂ population.

Two QTL mapping strategies implemented in software packages, GAPIT in the R (Lipka et al., 2012) and SVS (Golden Helix, Bozeman, MT, United States) were used to identify loci associated with pod shatter tolerance. For GAPIT analysis, we did not correct population structure using principal components in the F₂ mapping population. Linear marker regression analysis was performed to determine trait-marker associations in the SVS package. The same approach was also followed to reveal the genome-wide association between DArTseq markers and rupture energy among 83 accessions. For GWAS, we selected a set of 54,034 high quality markers which were genotyped across all accessions. To control spurious trait-marker associations, the first 10 eigenvectors (principal components) were calculated in the SVS package. Cryptic relatedness due to ancestry by descent was controlled with the Identity-by-Decent matrix (K matrix). The Mixed Linear Model (Price, 2006; Yu et al., 2006) adjusted with K-matrix and population structure matrix with PC1 – PC10 was used to test the trait-marker associations in the SVS package. The p-values were adjusted to control the false discovery rate (FDR) of 5%. The significance threshold was determined by applying Bonferroni correction [p = 0.05/6464 (total of markers mapped): 7.73515E_-06]. Trait-markers with significance ≤ log(₁₀)p of 5.11153 were ‘declared’ as true associations for pod shatter resistance in an F₂ population. Manhattan plots were generated in the SVS package.

The physical map positions of significant markers associated with pod shatter resistance were determined using the reference B. nigra, B. oleracea, B. juncea, and B. napus genomes by BlastN (Altschul et al., 1990) searches, as detailed in Raman H. et al. (2014). The physical positions of pod shatter resistance genes in A. thaliana (accessed from TAIR¹) were also determined by searching sequence identities with the reference genomes. The top blast significant hits (≥E^-10) were considered to infer the putative physical positions of markers/candidate genes on the reference genomes, while blast hits to multiple loci with the same top E value were considered to be unmapped onto the reference genome.

There were significant differences (p < 0.001) within the 83 B. carinata accessions tested with respect to pod rupture energy that ranged from 1.52 to 4.5 mJ in 2012, and 1.6 and 4.6 mJ in 2013 (Figures 1A,B). A positive strong correlation (r = 0.69) among accessions evaluated across both the 2012 and 2013 growing environments was observed, indicating that RE is genetically controlled (Figure 1C). Three B. carinata accessions, ATC94126, ATC94457, and ATC94458 had 9.14 to 9.63 times higher RE compared to the B. napus control genotype, BLN2762 (Table 1).

A set of 54,037 high quality DArTseq markers with call rates of >90% and a reproducibility of >95% were selected for genetic diversity and population structure analyses to determine whether shatter resistant sources are genetically diverse (Table 1). Hierarchical cluster analysis based on the Euclidean distance revealed five distinct groups at 60% similarity (Figure 2A). The cluster I was the largest with 75 accessions, followed by three accessions in cluster II (ATC94120, ATC93973, and ATC94192) and cluster IV (ATC90258, ATC94011, and ATC93888). Both cluster III (ATC94409), and cluster V (ATC94109) contained only one accession (Figure 2A). The overall genetic diversity among accessions was assessed with PCO analysis (Figure 2B), which revealed similar clustering. There were four clear groups with the majority of the accessions in cluster I. The first three coordinates (PC1 = 15.9%, PC2 = 5.3%, and PC3 = 4.3%) accounted a total of 25.39% of the genetic variation (Figure 2B), suggesting a weak population structure. The Bayesian – based clustering analysis using the maximum likelihood distribution LnP(D) of 83 accessions identified two sub-populations as shown in Figure 2C. The Wilcoxon test also revealed the presence of two subpopulations. Seventy nine accessions were in sub-population I and four accessions were in sub-population II. The STRUCTURE analysis supported the results of cluster analysis; all 83 accessions were grouped in two clusters at 90% similarity (Figure 2A).

Based on the pod shatter resistance (RE) scores, two single plant selections were made from accessions BC73526 (high RE) and BC73524 (low RE) to generate an F₂ population, representing cluster I (Figure 2A). Both parental lines of the F₂ mapping population from the cross, BC73524/BC73526 differed significantly from each other with respect to pod shatter resistance; the shatter prone, maternal parent (BC73524) had the lower RE of 2.2 mJ^(1/2) (4.8 mJ) and the resistant, paternal parent (BC73526) had the higher RE of 4.4 mJ^(1/2) (19.8 mJ; Table 1). The F₂ population showed a continuous distribution of RE scores, ranging from 2.2 to 4.7 mJ^(1/2) with the mean score of 2.71 mJ^(1/2) (Figure 3). This was typical for quantitative traits such as pod shattering resistance. In order to validate these F₂ RE scores, we evaluated the F_2:3 progenies (Supplementary Table S1). Our results showed that there was a strong positive correlation (r = 0.83) between RE scores of F₂ plants and their F₂:₃ progenies, suggesting that phenotypic scores in F₂ were accurate.

A total of 6,464 markers that showed polymorphism between the parents, and segregated in the complete set of F₂ population (300 lines) were selected for the genetic linkage map construction (Supplementary Table S1). All of the mapped markers were assigned to the 17 linkage groups, equivalent of haploid genome of B. carinata. Of them, 4,981 marker loci were located on the 8 linkage groups of B_c subgenome and 1,483 loci were on the 9 linkage groups of C_c subgenome, covering a total genetic distance of 1622.82 cM (Table 2). The marker density ranged from 1.07 (C4) to 7.35 (B01) with an average density of 3.98 cM. Chromosome C5 had the least number of markers (78) as compared to B4 (881). This genetic linkage map was further used for the QTL identification.

We identified five significant QTL (LOD = 3) associated with pod shatter resistance, Qpsr.wwai-B1a, Qpsr.wwai-B1b, Qpsr.wwai-B3, Qpsr.wwai-B8, and Qpsr.wwai-C5 in the BC73524/BC73526 population (Table 3, Figure 4A and Supplementary Table S1). Two QTL; Qpsr.wwai-B1a tagged with the in silico DArT marker 5863583, and Qpsr.wwai-B1b tagged with DArTseq-SNP marker 5858104|F| 0-14:A > T, were located 7.4 cM apart on chromosome B1. Other three QTL, Qpsr.wwai-B3, Qpsr.wwai-B8, and Qpsr.wwai-C5 were identified on chromosomes B3, B8, and C5, respectively (Table 3). Of these, Qpsr.wwai-B1b accounted for the maximum (5.27%) of the phenotypic variation and the Qpsr.wwai-C5 accounted for the least (3.71%). All five QTL explained a total of 23.73% of the phenotypic variation for RE. The shatter resistant parent, BC73526 contributed the favorable allele as envisaged by pod strength, and thus reduced pod shattering in progenies. DArTseq markers were assigned the physical positions on B. carinata genome, by comparing their sequence identities with the reference genomes of B. nigra, B. juncea, B. oleracea, and B. napus. Our results showed that the Qpsr.wwai-B1a, qPSR.wwai-B1b, Qpsr.wwai-B3, Qpsr.wwai-B8, and Qpsr.wwai-C5 were located to the pseudomolecules of B1, B3, B8, and C5, respectively (Supplementary Figure S2 and Table S2).

To establish whether pod length relates to pod shattering, we mapped QTL associated with pod length in the F₂ population. Our results showed that one significant marker, 5859132|F| 0–7:T > G at QTL (LOD = 3.55) was associated with pod length in the F₂ population (Figures 3, 4B, and Supplementary Tables S2, S3). This QTL was identified on chromosome B8 and mapped 1 cM apart from the pod shatter resistance QTL, Qpsr.wwai-B8. Two other markers, 5832583 and 5863583 on chromosome B1 also showed association with pod length (LOD scores of 2.7 to 2.9, Supplementary Table S3).

GWAS analysis using 54,034 markers based polymorphisms was performed to verify the alleles for pod shatter resistance in diverse B. carinata accessions Although, we used a small number of accessions for this analysis, we found 19 statistically significant SNP associations between markers and pod shatter resistance (RE scores) based on the Bonferroni corrected threshold –log 10(p) = 9.25292E^-07 (Supplementary Table S3). By controlling type 1 error using kinship coefficients (IBS) and first 10 principal components at least 16 consistent significant associations were identified across both 2012 and 2013 trials with LOD score of ≤5.35 (Supplementary Tables S3, S4).

Of 6464 DArTseq markers mapped, the chromosomal positions of 5,080 markers could be linked with the pseudomolecule positions to the published genome sequences of B. oleracea, B. napus, B. juncea, and B. nigra (Supplementary Table S1). We also anchored several scaffolds which have been unmapped yet to the pseudomolecules of B. juncea genome assembly in an F₂ population. Furthermore, marker sequences targeting QTL were aligned with the sequenced reference B, C and AC genomes and physical intervals harboring candidate genes for pod shatter resistance. Of the seven pod shatter resistance genes of A. thaliana searched, FUL – a MADS box gene negatively regulated by APETALA1 (TAIR ID: AT5G60910.1), was located 63.1 kb away from the significant SNP marker, 5863583 on chromosome 1B (Table 3). Other candidate genes were located 0.4 to 8.3-Mbp apart from corresponding QTL regions in the F₂ population (Supplementary Tables S4, S5). We identified 40 GWAS SNP associations (LOD ≤ 3) in the proximity (5.1 kb to 16 Mbp) of genes controlling pod shattering in A. thaliana. Of them, three orthologs of FUL were located on chromosome B1 (5.1 kb), B6 (97.4 kb) and on the LFLV01001230.1scaffold_28.1 of the reference genome of B. nigra (34.69 kb), while two orthologs of IND were located on chromosome B1 (53.69 kb) and B2 (96.77 kb). One ortholog of SHP2 was also identified within 77 kb region of chromosome B5 corresponding to SNP association with 100067358|F| 0-31:T > C marker (Supplementary Table S5).

Pod structure was observed (40 days after anthesis) under a fluorescence microscope to determine any link between the pod DZ differentiations and shatter resistance in B. carinata. The anatomical feature of parents displayed a distinctive difference in the valve margin formation (Figure 5). The shatter prone parent, BC73524 had the well-developed DZ comprising thin walled parenchymatous cells (Figure 5a) compared to the shatter resistant parent, BC73526 (Figure 5b). Thirty randomly selected F₂ plants exhibited a varied level of DZ development pattern (Figures 5d–i). For example, there was either clear DZ along the whole valve margin (Figure 5e), similar to shatter prone parent (Figure 5a); loss of DZ proximal to the main vascular bundle (mv) as well as near the outer part of the replum (Figure 5d), similar to the shatter tolerant parent (Figure 5c); and DZ proximal to the main vascular bundle (mv) but did not extend near the outer part of the replum (Figure 5h). In B. napus, a well-developed DZ was clearly evident (Figure 5c) similar to the shatter prone B. carinata parent, BC73524.