Quantitative Trait Loci for Heat Stress Tolerance in Brassica rapa L. Are Distributed across the Genome and Occur in Diverse Genetic Groups, Flowering Phenologies and Morphotypes

Heat stress events during flowering in Brassica crops reduce grain yield and are expected to increase in frequency due to global climate change. We evaluated heat stress tolerance and molecular genetic diversity in a global collection of Brassica rapa accessions, including leafy, rooty and oilseed morphotypes with spring, winter and semi-winter flowering phenology. Tolerance to transient daily heat stress during the early reproductive stage was assessed on 142 lines in a controlled environment. Well-watered plants of each genotype were exposed to the control (25/15 °C day/night temperatures) or heat stress (35/25 °C) treatments for 7 d from the first open flower on the main stem. Bud and leaf temperature depression, leaf conductance and chlorophyll content index were recorded during the temperature treatments. A large genetic variation for heat tolerance and sensitivity was found for above-ground biomass, whole plant seed yield and harvest index and seed yield of five pods on the main stem at maturity. Genetic diversity was assessed on 212 lines with 1602 polymorphic SNP markers with a known location in the B. rapa physical map. Phylogenetic analyses confirmed two major genetic populations: one from East and South Asia and one from Europe. Heat stress-tolerant lines were distributed across diverse geographic origins, morphotypes (leafy, rooty and oilseed) and flowering phenologies (spring, winter and semi-winter types). A genome-wide association analysis of heat stress-related yield traits revealed 57 SNPs distributed across all 10 B. rapa chromosomes, some of which were associated with potential candidate genes for heat stress tolerance.


Introduction
Brassica napus (oilseed rape, canola), a globally important oilseed crop, is vulnerable to heat and drought stress, especially during the early reproductive stage. B. napus has a relatively narrow gene pool due to its recent evolution from its ancestors B. rapa (field mustard, turnip, sarson, Asian cabbage) and B. oleracea (Mediterranean cabbage) [1], and as a result of 'bottlenecks' during original polyploidisation events, subsequent agricultural selection in spring and winter pools, selection for canola quality and breeding in isolated regional environments [2]. This situation was accentuated in Australia, where selection in a closed spring canola breeding population significantly reduced genetic diversity over 30 years (1970-2000) due to genetic drift [3]. Consequently, research has focused on the ancestors of B. napus (such as B. rapa) as potential sources of heat tolerance [4].
B. rapa is a traditional oilseed and vegetable crop that has been cultivated globally for more than 6000 years [5]; it has a wide geo-distribution with some types flourishing in heat and drought-affected regions [6]. The species includes a diverse range of agricultural morphotypes (leafy, rooty and oilseed types) and phenologies (winter, spring and semiwinter types) which arose after initial domestication from wild types in Central Asia [7] followed by independent selection in Europe and Asia [7][8][9][10][11][12].
B. napus, B. rapa and B. juncea were very sensitive to heat stress during flowering, which decreased subsequent pod and seed formation [13,14]. Temperatures greater than 29.5 • C during flowering in the field decreased seed yield in B. napus [15]. Young et al. [16] showed that high-temperature stress in B. napus (35 • C for 4 h each day for 1 or 2 weeks after the initiation of flowering) reduced fruit and seed development, pollen germination and in vivo pollen tube growth. Using similar transient temperature stress for 7 days after the first open flower, Annisa et al. [4] found genetic variation in pod and seed number following heat stress in six spring-type B. rapa accessions. Chen et al. [17] found that transient heat stress negatively affected male and female reproductive organs of B. napus, with female organs more sensitive than male organs to heat stress.
Based on our previous experience that heat stress has its greatest impact during the reproductive stage in B. rapa and B. napus [4,13,14,17], we evaluated a genetically diverse global collection of B. rapa germplasm for heat stress tolerance following simulated transient daily heat stress after the first open flower on the main stem of well-watered plants. We included a range of morphotypes (leafy, rooty and oilseed types) and flowering phenologies (spring, winter and semi-winter types). We assessed the genotypic diversity of the population using SNP markers. We used a genome-wide association analysis (GWAS) to investigate significant quantitative trait loci (QTL) for heat tolerance and potential relationships between heat tolerance and geographic origin, morphotype and/or flowering phenology. We searched for candidate genes associated with significant SNPs which may contribute to heat tolerance.

Plant Materials
A total of 217 B. rapa lines were used in this study (Table S1). These accessions were chosen based on their self-compatibility reported in Annisa et al. [10] and Guo et al. [12] and represent a wide range of genetic diversity groups, morphotypes and flowering phenologies. Of these, 75 lines with a limited number of seeds and with self-incompatibility were excluded for phenotyping experiment, so 142 lines were phenotyped for several physiological and yield-related traits.

Temperature Treatments
Two sets of plants were sown in pots, with a single plant in each pot. Set 1 plants were grown for the control (no heat stress) treatment and Set 2 plants were grown for the heat stress treatment at first flower on the main stem. There were three replications (pots) of each line in each set. Plants were grown in 8.1 L pots 230 mm in depth (Garden City Plastics, Perth, Australia) in a greenhouse at The University of Western Australia, Crawley, Western Australia (31 • 57 S, 115 • 47 E) with an average relative humidity of 65%, 25 • C/15 • C (day/night) temperature and average light intensity 635 mmol m −2 s −1 PAR under optimal moisture conditions. Light readings in the glasshouse were taken at mid-day on a sunny day in mid-spring in an unshaded area of the glasshouse. Each pot contained 4.5 kg of canola potting mix, comprising 50% fine composted pine bark, 20% coco peat and 30% brown river-sand plus 1.0 g of gypsum per kg with its final pH of 6.0. After sowing, each pot was fertilised every two weeks with 0.25 g of soluble nutrient powder (Thrive, Yates Australia, NSW, Australia; containing 25% N, 5.0% P, 8.8% K, 4.6% S, 0.5% Mg, 0.18% Fe, 0.01% Mn, 0.005% Cu, 0.004% Zn and 0.001% Mo). When the first open flower was seen on the main stem, plants were moved to a controlled environment room (CER) for seven days of temperature treatment, and the first five floral buds on the main stem (flowers opened on the first day) were tagged. Set 1 plants were transferred to CER1 (control treatment) and Set 2 plants to CER2 (heat stress treatment). After 7 d of temperature treatment, the plants were returned to the greenhouse and maintained with adequate nutrition and moisture until seed harvest at maturity. Both CER1 and CER2 were set to 16 h day and 8 h night with 425 mmol m −2 s −1 light and 65% relative humidity. CER1 (control) was set at 25 • C day/15 • C night. CER2 (heat stress) was set at 25 • C night, and during the day the temperature was gradually increased in the first 6 h to 35 • C, maintained at 35 • C for 4 h and then gradually decreased over the next 6 h to 25 • C. Plants in both treatments were watered regularly through low-flow (35 mL h −1 ) drippers to maintain the soil water content at not less than 90% field capacity. Both high temperature and control temperature CERs had the same light intensity, same day-length, same humidity, same nutrition and adequate water to avoid any confounding of heat stress treatment with other effects that might have differed between the treatments.
Both sets of plants, destined for high and control temperature treatments, received the same growth conditions before, during and after the temperature treatments to avoid any confounding of heat stress with other growth inputs such as water, nutrients and light.

Phenotyping
Several physiological traits were measured between 11.00 and 13.00 h on the seventh day of temperature treatment ( Table 1). The temperature of a recently opened flower and a nearby leaf was measured using an infrared thermometer (Impac ® Model IN 15 plus, LumaSense Technologies GmbH, Santa Clare, CA, USA) with a minimum 2.2 mm diameter measurement area. A separate digital thermometer with a 1 s response time measured the ambient temperature. The leaf and bud temperatures and the ambient temperature were recorded simultaneously with four repeat measurements per leaf and bud. The temperature difference between bud and ambient environment (T1, • C), leaf and ambient environment (T2, • C) and bud and leaf (T3, • C) were calculated. Stomatal conductance of the youngest fully expanded leaf was measured using an SC1 leaf porometer (Decagon Devices, Pullman, WA, USA) on adaxial and abaxial leaf surfaces. Leaf conductance (LC, mmol m −2 s −1 ) was the sum of the adaxial and abaxial conductance. The chlorophyll content index (CI, SPAD unit) was measured on intact leaves with a portable chlorophyll meter SPAD-502Plus (Konica Minolta Sensing Americas, Ramsey, NJ, USA). CI is the ratio of the leaf transmittance in red light at 650 nm (at which chlorophyll absorbs) and in near-infrared light at 940 nm (for the correction of leaf thickness).
Several yield-related traits were also measured at maturity. Each plant was cut at the soil level and dried at 32 • C for 14 d before measuring the above-ground biomass (BM, g). The seed pods from each plant were harvested at maturity, threshed manually and cleaned using a vacuum separator (Kimseed, Wangara, WA, Australia). Seed yield from the first five floral buds on the main stem (Y5P, g) and seed yield on each whole plant (YWP, g) were measured. The harvest index (HI, %) was calculated based on the ratio of YWP and BM for each plant.
Line performance in the control treatment was identified by the suffix '_C' (e.g., BM_C), and line performance in the heat stress treatment was identified by the suffix '_H' (e.g., BM_H). Line performance was also compared across C and H treatments. Line performance in control (Y c ) and heat stress (Y h ) treatments, and the mean performance of all lines in the control treatment (Y c ), were used to calculate the stress tolerance index (STI) [18] (Equation (1)) and percentage change from the control to heat stress treatment (%C) (Equation (2)) as follows:

Statistical Analysis of Phenotypic Traits
Statistical analysis of phenotypic traits was based on linear mixed models with residual maximum likelihood estimation in ASReml-R version 4 [19]. For all the measured traits, the temperature treatment ('Treatment') was considered a fixed effect and Line, Line × Treatment and Error as random effects. The significance of fixed effects was evaluated by Wald statistic and random effects by Z statistic of variance components. Principal component analysis (PCA) biplot and clustering were performed using the 'factoextra' and 'FactoMineR' packages in R Version 4.0.2.

SNP Genotyping, Population Structure and Genetic Diversity Analysis
Leaf samples were combined from three plants at the four-leaf stage for each B. rapa line to extract genomic DNA using the Qiagen DNeasy plant kit (Düsseldorf, Germany). SNP genotyping was assessed on a B. napus Illumina Infinium 6K array described in Dalton-Morgan et al. [20]. Chips were scanned using an Illumina HiScan SQ machine. Samples with their call rate below 70% were excluded from further data analysis. SNP locations were established by BLASTing to the B. rapa genome sequence (version 1.5) [21]. All SNPs located on the B. rapa A genome with amplification in >90% of experimental lines and a minor allele frequency >10% across the diversity set were selected for further analysis.
The software package STRUCTURE version 2.3.4 [22] was used to detect population structure [23]. Kinship (K) was calculated with SPAGeDi 1.3 [24,25]. A hierarchical cluster analysis was performed using the unweighted pair group method with arithmetic averages (UPGMA) as proposed by Sneath and Sokal (1973), and the ordination analysis was performed using PRIMER 6 software [26].

Genome-Wide Association Mapping
A general linear model (GLM) and a mixed linear model (MLM) were compared in TASSEL [27] to identify QTL significantly associated with yield and heat tolerance traits and to assess the effect of the kinship matrix for family relatedness estimates (K), principal component analysis (P) and population structure (Q). We compared three different models in GLM: (1) the naive model (no control for kinship or population structure), (2) the Q model (controls for population structure) and (3) the P model (controls for principal components) with the top 10 principal components included as fixed effects [28,29]. We also compared three different models in MLM: (1) a model based on kinship (K), (2) a model that unified population structure and kinship (Q + K) and (3) a model based on principal component and kinship (P + K). The p-values obtained from all models were converted into −log10 (p). The variation of observed −log10 (p) for each SNP from marker-trait associations was compared with the expected p-values in quantile-quantile (QQ) plots. We evaluated these six models for false positives and false negatives based on the QQ plots. A sharp deviation from the expected p-value distribution in the tail area would indicate that a model appropriately controlled both false positives and false negatives [30]. The significance of the association between markers and traits was based on the threshold −log10 (p) = 3.0.

Candidate Gene Discovery
The flanking SNPs of each QTL were used to blast the reference genome sequence of B. rapa (version 1.5) [21] and the A genome of B. napus using the Darmor-bzh genome reference [1] to locate the position of SNPs. Genome annotation files were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/, accessed on 8 January 2021). All genes within 100 bp of the QTL were selected for further analysis. The gene function and other information were evaluated, including the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology, gene ontology (GO) component, GO function and GO process.

Phenotypic Variation and Heat Stress Tolerance
Of the 142 B. rapa lines observed for five physiological traits measured on day 7 of the temperature treatment and four yield-related traits measured at maturity, eight lines showed self-incompatibility and were excluded from the analysis.
The main effect of temperature treatment ('Treatment') was highly significant for all nine traits measured ( Table 2), indicating that heat stress significantly increased or decreased the mean value of the trait across all lines relative to the control treatment. For example, the percent change (%C) under heat stress relative to the control in YWP, HI and Y5P, averaged across all lines, was -39.09%, -50.03% and -98.08%, respectively, while %C in BM was +8.64% under heat stress ( Table 2). The %C in T3, CI and LC were +85.50%, +5.28% and +73.47%, respectively; that is, the high-temperature treatment enhanced the physiological activities of the plants in the absence of moisture and nutritional stress and increased the accumulation of dry matter (BM). However, heat stress inhibited the reproductive process and %C for seed yield (both YWP and Y5P), and HI was negative ( Table 2).
There were significant or highly significant interaction effects of Line × Treatment for all nine traits measured ( Table 2), suggesting that the performance rank of these 134 B. rapa lines varied in the control and heat stress treatment for all nine traits.
PCA was based on %C and STI for Y5P, YWP, BM and HI in 134 lines ( Figure 1). The first principal component (PC1) explained 30.61% of the variance and was strongly associated with %C and STI of YWP and HI. The second component (PC2) explained 19.67% of the variance and was negatively associated with Y5P_STI and Y5P_%C. PC2 separated six lines, particularly R101 and R090, with relatively higher YWP and Y5P under heat stress.
Cluster analysis, also based on %C and STI for Y5P, YWP, BM and HI in 134 lines, revealed two major clusters, heat-tolerant (Tol) and heat-sensitive (S), each of which had two sub-clusters (Tol1, Tol2 and S1, S2) ( Figure 2). Tol1 and Tol2 had an average %C in YWP of -8.38% and -22.81%, respectively, while S1 and S2 had average %C in YWP reduction of -53.67% and -42.48%, respectively (Table 3). Tolerant groups showed a greater increase in BM under heat stress than sensitive groups and experienced less reduction in HI. Tolerant groups had higher STI (higher overall performance) for HI than sensitive groups. Tol2 suffered the least impact of heat stress on Y5P (yield of five pods on the main stem) ( Table 3). One line, R101 (a yellow sarson from India), had relatively high YWP and Y5P under heat stress ( Figure 1, Table S2). Table 2. Variance components and predicted values of four seed yield-related traits measured at maturity and five physiological traits measured during heat stress at the early flowering stage of 134 B. rapa lines. The Line and Line × Treatment interaction effects were considered random, and the Z test was applied to test for significance. Treatment was considered a fixed effect and the Wald test was used to test for significance. The predicted value ± standard error of each trait under control (C) and heat stress (H) conditions and the percentage change (%C) under H compared to C are also provided.    Five physiological traits measured on the 7th day of temperature treatment also showed a similar trend ( Table 3). The tolerant groups had cooler buds (T1) and leaves (T2) than the sensitive groups, and Tol2 had the coolest buds and leaves. The temperature difference between bud and leaf (T3) in heat-tolerant lines also showed higher %C than in the heat-sensitive lines (Table 3), consistent with our previous reports [13,14]. Tolerant groups had better photosynthesis and transpiration performance than heat-sensitive groups, as reflected by the relatively higher chlorophyll content index and leaf conductance (Table 3).

SNP Genetic Diversity and Genetic Population Structure versus Morphotype, Flowering Phenology and Heat Tolerance and Sensitivity
All 217 B. rapa lines were subjected to SNP genotyping, and the average call rate of all lines was 87.04%, and five lines with their call rate below 70.0% were excluded from further data analysis. A total of 1602 SNP alleles were polymorphic across 212 lines, and these SNPs were used in a population structure analysis. The optimum solution in STRUCTURE was four genetic sub-populations (Figures 3 and 4), and this was supported by two major populations each with two sub-populations in a phylogenetic cluster analysis based on PRIMER 6 ( Figure 5).
These four genetic sub-populations reflected the geographical origin of lines: A1 from East Asia (China, South Korea, North Korea and Japan), A2 from South Asia (India, Pakistan, Afghanistan, Bangladesh and Nepal) and B1 and B2 from mainly European countries. Various morphotypes (leafy, rooty and oilseed types) were distributed across each of the four sub-populations, as were winter, semi-winter and spring types ( Figure 5 and Table S1). The heat-tolerant lines in Tol1 and Tol2 and sensitive lines in S1 and S2 were also randomly distributed across the four genetic sub-populations ( Figure 5).   SNPs. Of these, 134 lines were classified into two heat-tolerant groups (Tol1 and Tol 2) and two heat-sensitive groups (S1 and S2) based on cluster analysis and principal component analysis of four yield-related agronomic traits measured at maturity (see Figures 1 and 2). The flowering phenologies based on cold treatment requirements and days to flowering in the screen-house are spring (S), winter (W) and semi-winter (SW) types. The morphotypes are leafy (L), rooty (R) and oilseed (O) types.

Genome-Wide Association Analysis and Candidate Genes Related to Heat Stress Tolerance
In order to remove or reduce skewness and normalise the distribution of the data, some phenotypic traits were transformed before the marker-trait association analysis. Log transformation was applied to Y5P_C, HI_C, Y5P_H, YWP_H and HI_H, and cube root transformation was applied to BM_C and YWP_C.
After comparing six models, an MLM based on population structure and kinship Q+K was selected to identify QTL significantly associated with yield and heat tolerance traits. The MLM(Q+K) model showed a relative sharper deviation from the expected P-value distribution in the tail area for all yield-related traits ( Figure S1).
A total of 57 SNP markers were significantly associated with at least one of the yieldrelated traits with a logarithm of odds (LOD) score greater than 3.0 (Table 4). Five SNPs were associated with YWP_%C and distributed on four chromosomes (A03, A04, A05 and A07). A single SNP, UQnapus5640, was associated with YWP_STI at LOD 4.37 and explained 19.06% of the phenotypic variance (Table 4).
Among the 57 significant SNPs, the range of phenotypic variance explained ranged from 6.52% for UQnapus1680 in relation to BM_H to 23.43% for UQnapus3633 in relation to HI_C. Most SNPs explained <10% of the total phenotypic variance, indicating that heat tolerance is a complex quantitative trait, with several QTLs for each trait.
Under heat stress conditions, 24 SNPs were associated with seed yield-related traits and distributed on all 10 B. rapa chromosomes (Table 4). Under control conditions, 33 SNPs were associated with seed yield-related traits and distributed on all 10 chromosomes (Table 4).
One SNP marker, UQnapus2563 located on Chromosome A05, was associated with two traits. This marker explained 20.35% of the phenotypic variance for HI_C with an LOD of 6.31 and 15.9% of the phenotypic variance of HI_H with an LOD of 4.58 (Table 4). A search of the A genome of B. rapa and B. napus in a 230 kb genomic region flanking UQnapus2563 revealed 49 candidate genes present in both species (Table S3).
Another SNP marker, UQnapus2473 located on Chromosome A07, was associated with Y5P_H at an LOD of 4.69 and explained 15.67% of the phenotypic variance. UQnapus2473 was also associated with YWP_C at an LOD of 3.71, explaining 10.81% of the phenotypic variance (Table 4). A search of the genomic region of 225 kb flanking UQnapus2473 found nine candidate genes in the B. napus A genome, eight of which were also in the B. rapa genome (Table 5). These nine candidate genes were further searched for orthologs in Arabidopsis thaliana. Function annotation revealed that several candidate genes, including the protein kinase gene AT2G18890.1 and several C2H2 zinc finger protein genes, were related to the abiotic stress response. Candidate gene AT2G18510.1 is involved in embryo development ending in seed dormancy. Candidate gene AT2G18500.1 encodes ovate family protein 7 (OFP7) located in the plasma membrane and expressed in embryos, flowers and seeds (Table 5). Table 4. Fifty-seven SNP markers were associated with at least one yield-related trait with the LOD ≥ 3.0 under control (C) and/or heat stress (H) conditions. The four yield-related traits are above-ground biomass (BM), seed yield whole plant (YWP), harvest index (HI) and seed yield of five pods (Y5P). The stress tolerance index (STI) and percentage change (%C) based on YWP are considered for heat tolerance. The LOD value calculated by -log 10 (p) of each associated SNP marker and the phenotypic variance explained are also listed. Markers which were associated with two or more traits are in bold and italic format.

Discussion
Heat tolerance is an important trait for Brassica grain crops exposed to high temperatures during flowering to avoid grain yield losses due to heat stress [13][14][15][16]. Heat stress events during crop flowering are expected to increase in frequency in the future due to global warming [31]. This study simulated heat stress during flowering in a globally diverse collection of B. rapa lines and compared their responses to transient daily heat stress relative to a control temperature treatment for 7 d from the first open flower. The plants were provided with adequate moisture and nutrients to avoid other stresses during the temperature treatment. As a result, several promising heat-tolerant B. rapa lines were identified, which maintained grain yield, biomass and harvest index after 7 d of heat stress at the first flower, while heat-sensitive lines suffered a significant reduction in yield-related traits in the heat stress treatment relative to the control.
B. rapa is a traditional oilseed and vegetable crop with broad genetic diversity and wide geographic distribution, with some types flourishing in heat-affected regions. In previous research, a leafy vegetable morphotype of B. rapa from Indonesia tolerated heat stress during flowering and early seed fill [4]. Hence, we expanded our survey of heat stress tolerance in B. rapa to include leafy, rooty and oilseed morphotypes with spring, winter and semi-winter flowering phenologies from diverse geographic origins. Forty-one out of 134 lines were heat tolerant and distributed across four SNP sub-populations identified from a cluster analysis based on 1602 polymorphic SNP markers with geographic origins in East Asia (A1), South Asia (A2) and Europe (B1 and B2). Heat-tolerant lines were found among all morphotypes (leafy, rooty and oilseed) and flowering phenologies (spring, winter and semi-winter types) ( Figure 5).
Phylogenetic analyses of the 212 B. rapa lines in this study follow previous results on this population using microsatellite markers [4,12]. The various morphotypes (leafy, rooty and oilseed type) were distributed across each of the four SNP sub-populations, as were winter, semi-winter and spring flowering phenologies ( Figure 5). The four sub-populations based on SNP genetic diversity reported in this study reflect the geographical origins of lines but not their morphotype or flowering phenology. This is consistent with reports [11,12] in which similar morphotypes from different regions were often not genetically related, but contrasts with Bird et al. [32] who reported five SNP sub-populations in B. rapa, which tended to be divided by morphotype and geographic origin.
Genome-wide association studies of heat stress-related yield traits revealed 57 SNPs, with some associated with potential candidate genes for heat stress tolerance. Identifying SNP markers and candidate genes associated with heat stress-affected traits helps select stress-tolerant accessions in breeding programs. The 57 SNPs were distributed across all 10 B. rapa chromosomes, with the markers explaining between 6.52% (UQnapus1680 for BM_H) and 23.43% (UQnapus3633 for HI_C) of the phenotypic variance. On average, each SNP explained a small percentage of phenotypic variance. Heat stress tolerance is a complex quantitative trait, with no major genes found in this study.
Complex quantitative genetic control of heat stress and other abiotic stress tolerance has been reported in Brassica [33] and other crop species [34]. A study of heat stress tolerance in spring-type B. napus under controlled conditions reported that heat stress increased pollen sterility and sterile/aborted pods and reduced pod numbers on the main stem [35] but did not investigate seed yield, biomass or harvest index.
Nine candidate genes flanking UQnapus2473 were found in the B. rapa genome and the B. napus A genome (Table 5). Function annotation of orthologs in A. thaliana revealed that several candidate genes, such as protein kinase gene AT2G18890.1 and several C2H2 zinc finger protein genes, were related to abiotic stress response. AT2G18510.1 is involved in embryo development, and the embryo development stage is very sensitive to heat stress in B. napus [14]. AT2G18500.1 encodes ovate family protein 7 located in the plasma membrane and expressed in embryos, flowers and seeds. Ovate family proteins are a family of plant growth regulators that play diverse roles in many aspects of physiological processes. Ovate family protein 1 regulates the drought stress response in Populus trichocarpa [36].
However, further functional confirmation research is needed to validate these candidate genes in response to heat stress in B. rapa. Once confirmed, breeders can potentially transfer valuable heat stress tolerance QTL from B. rapa to B. napus.
Future research should use a higher density linkage map to minimise the QTL region on the genome limiting the number of candidate genes and undertake functional characterisation of candidate genes via qRT-PCR. Eventually, it should be possible to use functional markers to assist the screening of heat-tolerant germplasm and marker-assisted introgression of multiple QTL from B. rapa to B. napus for heat tolerance.

Conclusions
B. rapa is a traditional oilseed and vegetable crop with broad genetic diversity and wide geographic distribution. Large genetic variation for heat tolerance and sensitivity was found for above-ground biomass, whole plant seed yield and harvest index and seed yield of five pods on the main stem at maturity. Phylogenetic analyses confirmed two major genetic populations from East/South Asia and Europe. Heat stress-tolerant lines were distributed across each population and included diverse geographic origins, morphotypes (leafy, rooty and oilseed) and flowering phenologies (spring, winter and semi-winter types). A genome-wide association analysis of heat stress-related yield traits revealed 57 SNPs distributed across all 10 B. rapa chromosomes, some of which were associated with potential candidate genes for heat stress tolerance. In summary, QTL for heat stress tolerance in B. rapa are distributed across the genome and occur in diverse genetic groups, flowering phenologies and morphotypes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13020296/s1. Figure S1: Comparison of six models for variation of observed p values (x-axes) against the expected p values (y-axes) using a general linear model (GLM) and mixed linear model (MLM); Table S1: List of 217 lines tested in this study; Table S2: Predicted value of 134 lines for nine traits under control (C), heat stress (H), and their percentage change (%C) and stress tolerance index (STI); Table S3