A novel motif in the 5'-UTR of an orphan gene 'Big Root Biomass' modulates root biomass in sesame.

Developing crops with improved root system is crucial in current global warming scenario. Underexploited crops are valuable reservoirs of unique genes that can be harnessed for the improvement of major crops. In this study, we performed genome-wide association studies on seven root traits in sesame (Sesamum indicum L.) and uncovered 409 significant signals, 19 quantitative trait loci containing 32 candidate genes. A peak SNP significantly associated with root number and root dry weight traits was located in the promoter of the gene named 'Big Root Biomass' (BRB), which was subsequently validated in a bi-parental population. BRB has no functional annotation and is restricted to the Lamiales order. We detected the presence of a novel motif 'AACACACAC' located in the 5'-UTR of BRB in single and duplicated copy in accessions with high and small root biomass, respectively. A strong expression level of BRB was negatively correlated with high root biomass, and this was attributed to the gene SiMYB181 which represses the activity of BRB by binding specifically to the single motif but not to the duplicated one. Curiously, the allele that enhanced BRB expression has been intensively selected by modern breeding. Overexpression of BRB in Arabidopsis modulates auxin pathway leading to reduced root biomass, improved yield parameters under normal growth conditions and increased drought stress sensitivity. Overall, BRB represents a solid gene model for improving the performance of sesame and other crops.

Introduction

Climate change‐related risks coupled with the rapid population growth and the rarity of natural resources are posing serious threats on food and nutritional security. This has catalysed global efforts to increase crop productivity and resilience through various crop improvement strategies. Among these, molecular plant breeding plays a central role in the identification of the genetic basis of important agronomic traits and their deployment for the development of improved cultivars (Moose and Mumm, 2008; Varshney et al., 2006). However, modern crop breeding has mainly focused on the above‐ground plant components, largely neglecting the hidden half of the plant, that is the root system (Voss‐Fels et al., 2018). As the root system plays a crucial role in anchoring the plant, exploring and exploiting soil resources such as nutrients and water to sustain growth and productivity under various environmental conditions (Lynch, 2013; Wang et al., 2015; Svacina et al., 2014; Xie et al., 2017; Liu et al., 2018; Su et al., 2019), more studies and efforts are needed to deliver crops with improved below‐ground traits.

The complex nature of root traits and the opaque nature of soil are the main factors that confound large‐scale and in‐depth studies on plant root system. Nonetheless, by adapting artificial growth conditions, phenotyping platforms and novel data processing tools, the genetic architecture of root traits is progressively uncovered in many plant species (Hochholdinger and Tuberosa, 2009; Mai et al., 2014; Downie et al., 2014; Atkinson et al., 2019). For example, using hydroponics system and root imaging, Kitomi et al. (2018) identified two quantitative trait loci (QTL) associated with root length in rice. The gel imaging platform was employed to detect several QTLs controlling root traits in maize (Zurek et al., 2015) and rice (Topp et al., 2013). Beyer et al. (2019) recently employed a paper roll‐supported hydroponic system to perform genome‐wide association studies (GWAS) in wheat that led to the identification of 68 marker–trait associations associated with various root traits. Transparent Plexiglas nailboard sandwiches filled with glass bead, and a nutrient solution was adapted to grow a rice diversity panel and perform GWAS for root traits (Courtois et al., 2013). Several studies were also conducted in Arabidopsis thaliana to dissect the genetic network controlling root traits by using petri dish or hydroponics (Gifford et al., 2013; Meijón et al., 2014; Kobayashi et al., 2015; Jia et al., 2019; Ogura et al., 2019). In wheat, the filter paper/polycarbonate screening plates and the ‘clear pot’ methods were used to discover major QTLs for root system architecture traits (Alahmad et al., 2019; Maccaferri et al., 2016).

Discoveries made so far on the genetic architecture of root traits have been essentially focused on major crops and model plants while minor crops have been overlooked. Besides their role in the diversification of food production, minor crops represent untapped gene pool, which could be harnessed for improvement of major crops (Mayes et al., 2012). Sesame (Sesamum indicum L.) is a typical example of an underexploited crop with special attributes such as high‐quality vegetable oil (Anilakumar et al., 2010; Miyahara et al., 2001; Noguchi et al., 2001; Sankar et al., 2005), good adaptation to marginal lands (Langham, 2007) and resistance to abiotic stresses such as drought and salt stresses (Dossa et al., 2019a, 2019b; Islam et al., 2016; Lakhanpaul et al., 2012). In particular, its extensive root system is thought to play a foremost role in its adaptation to different environments (Dossa et al., 2017b, 2017a; Gloaguen et al., 2019; Lakhanpaul et al., 2012; Langham, 2007). In addition, Su et al. (2019) recently reported a positive contribution of high root biomass to increased seed yield in sesame. However, the genetic basis of sesame root system is still unknown.

In a previous study, we designed and successfully adapted a hydroponic growth platform to demonstrate the high variability of root traits in sesame (Su et al., 2019). The availability of the genome sequence (Wang et al., 2014) and high‐quality single‐nucleotide polymorphism genotyping data of 705 diverse accessions (Wei et al., 2015) offer the opportunity to explore the genetic variants controlling the natural variation of sesame root traits. Herein, we performed GWAS based on a large‐scale phenotyping of root traits in a sesame diversity panel. We identified key genomic regions and several candidate genes associated with root traits in sesame. In particular, we identified and functionally characterized a novel gene named ‘Big Root Biomass’ (BRB), which modulates root biomass in sesame.

Result

Natural variation of root traits in sesame

An association panel composed of 705 diverse sesame (Sesamum indicum L.) accessions was previously assembled in our group and used for the detection of genetic variants governing several agronomic traits (Wei et al., 2015; Li et al., 2018; Zhou et al., 2018; Dossa et al., 2019a). To investigate variation in root traits in this study, we extracted 327 landraces and modern cultivars with large morphological diversity (Wei et al., 2015) and originating from 28 different countries in Africa, Asia, Europe and America (Table S1). A total of seven root traits, including root number (RN), root volume (RV), root length (RL), root surface area (SA), root dry weight (RDW), shoot dry weight (SDW) and root–shoot ratio (RSR), were evaluated on 1‐month‐old seedlings grown hydroponically under controlled conditions.

Large variation and normal distribution of all the assayed traits were observed in the association panel (P < 0.001) (Table S2; Figure S1). The broad‐sense heritability (H) values of the traits were moderate to high, ranging from 0.58 to 0.83 (Table S2). Moreover, strong correlations were observed between the traits (Figure S2), suggesting the existence of pleiotropic loci.

GWAS for sesame root traits

To uncover the genetic variants underlying the root traits in sesame, we performed a genome‐wide association study (GWAS) using one million high‐quality single‐nucleotide polymorphisms (SNPs) with a minor allele frequency ≥ 0.05 (Wei et al., 2015). The mixed model taking into account the population structure and parental relatedness was employed (Dossa et al., 2019a). We identified 409 SNPs significantly associated with the different traits at P < 1 × 10−5 (Table S3; Figure S3). Most of these associated loci were distributed in clusters across ten linkage groups (LG). The significant SNPs were resolved into 19 quantitative trait loci (QTL) based on the estimated linkage disequilibrium (LD) decay (88 kb) in the sesame genome (Wei et al., 2015; Dossa et al., 2019a) (Figure 1a; Table S3). We identified seven pleiotropic QTLs associated with up to five different traits (Table S3). The top five peak SNPs for each trait contributed from ~ 16.7 (RSR) to 67% (RN) of the phenotypic variation in the entire population (Table S3; Figure S4). Collectively, the QTLs harboured 318 genes (Table S4), 30 of which contained in total 90‐genic and 59‐promoter significant SNPs (Table S5). Candidate genes in the different QTLs were inferred by prioritizing the genes containing significant associated SNPs or closer to the peak SNPs (Dossa et al., 2019a). In total, 32 candidate genes controlling root traits and involved in various biological functions related to hormone signalling, transcription regulation, growth and development, cell organization, were uncovered in this study (Figure 1b; Table S6).

Figure 1

Genome‐wide association study (GWAS) for root traits in sesame. (a) Circos plot summarizing the GWAS results; (b) Gene Ontology enrichment analysis of the candidate genes identified for root traits in sesame; (c) Manhattan plot for root number and root dry weight traits. A strong association was co‐detected for the two traits on the linkage group 15, and the regional plot displays the top SNP (SNP5024573) located in the gene SIN_1025576.

BRB gene is significantly associated with root number and root dry weight

The most significant SNPs associated with the root traits were located on the LG15 within the co‐localized QTLs qRN15.2 and qRDW15.1 (Figure 1c; Table S3). The peak locus SNP5024573 (P < 6.6 × 10−10) fell within the promoter region of the gene SIN_1025576, dubbed BRB (‘Big Root Biomass’). BRB contained one synonymous SNP SNP5022360 (A/T) located at 393 bp in the coding region, one SNP SNP5022842 (C/G) located at −90 bp in the 5’‐UTR and 11 SNPs topped by the SNP SNP5024573 (G/A) located at −1821 bp in the promoter (Figure 2a; Table S7). The two loci, SNP5022842 and SNP5024573, were found in complete LD, resulting into two haplotypes (Hap1 and Hap2). Accessions with Hap1 representing the smallest group (n = 63) had significantly higher RN and RDW (P < 0.001) than those with Hap2 (n = 253) (Figure 2b,c).

Figure 2

Natural variations in the gene BRB (SIN_1025576) are associated with root traits. (a) Schematic representation of the significant SNPs located in coding and promoter regions of BRB. The SNP5022360 (A/T) was synonymous SNP so was excluded from further analyses. The two other SNPs, SNP5022842 (C/G) and SNP5024573 (G/A), were in perfect linkage disequilibrium leading to two haplotypes (Hap1 and Hap2). Sesame genotypes with Hap1 display dense root biomass while genotypes with Hap2 exhibited small root biomass; (b) root dry weight was significantly higher in Hap1 genotypes than those with Hap2; (c) root number was significantly higher in Hap1 genotypes than those with Hap2; (d–e) association analysis of root dry weight and root number traits in 197 F2:3 sesame population. The significant associated InDels and SNPs were marked in red triangle and dot, respectively. (f–g) sesame lines with a nucleotide deletion at the peak InDels (InDel5022798 and InDel5022804) displayed significantly high root dry weight and root number than those with a nucleotide insertion. Error bars indicate the SD (*** significant difference at P < 0.001).

To further confirm that variants in BRB modulate root biomass in sesame, we developed 197 F2:3 population from the accessions G340 (RDW = 0.121 ± 0.04 g; RN = 1137±109) and G441 (RDW = 0.040 ± 0.019 g; RN = 566±72) having the Hap1 and Hap2, respectively. We tried unsuccessfully to sequence 2 kb promoter region and the coding region (417 bp) of BRB in the F2:3 population. However, after decreasing the promoter size, we sequenced and analysed 292 bp upstream of the translation start codon (ATG), the full coding sequence and 180 bp downstream of the stop codon (TGA). A total of 146 SNPs and 53 insertions and deletions (InDel) were identified. Four InDels (InDel5022798, InDel5022796, InDel5022797 and InDel5022799) topped by IndDel5022798 located at −46 bp in the 5’‐UTR of BRB contributed significantly to the variation of RDW in the F2:3 population grown in pots containing soil as substrate (P < 1.67 × 10−4) (Figure 2d; Table S8). All these InDels were found in strong LD (r 2 > 0.9). Similarly, a bunch of six InDels (InDel5022804, InDel5022801, InDel5022800, InDel5022798, InDel5022799 and InDel5022796) and the GWAS‐based identified significant locus SNP5022842 were found significantly associated with RN, peaked by InDel5022804 located at −52 bp upstream of the ATG (P < 2.3 × 10−4) (Figure 2e; Table S8). All these loci were also found in strong LD (r 2 > 0.9). Sesame lines with a nucleotide deletion at the peak InDels (InDel5022798 and InDel5022804) displayed significantly high RDW and RN (P < 0.001) than those with a nucleotide insertion (Figure 2f,g). Taken together, our results demonstrated that variants in the 5’‐UTR of BRB, particularly the InDels, contribute to the natural variation of RN and RDW in sesame, independently of the growing media.

The significant InDels correspond to a novel motif in the 5’‐UTR altering the expression of BRB

Since the detected InDels contributing to the variation of RDW and RN in sesame were closely located in the 5’‐UTR of BRB based on the F2:3 population analysis, we sequenced and analysed the 5’‐UTR in 50 accessions from the association panel, including 25 accessions with Hap1 and 25 accessions with Hap2. We found that all the clustered significant InDels for RDW and RN correspond to the presence/absence of a 9 bp sequence duplicated motif ‘AACACACAC’ (Figure 3a). Blast search in various databases (http://www.dna.affrc.go.jp/PLACE/; http://arabidopsis.med.ohio‐state.edu) for the identification of the motif did not return any hit, so we presumed that it is an unreported or an uncharacterized motif. All sesame accessions with Hap1 displayed the single AACACACAC motif while accessions with Hap2 had a duplicated motif (18 bp), indicating that the duplication pattern of the motif is associated with variations of RDW and RN in sesame (Figure S5).

Figure 3

A motif in the 5’‐UTR alters the expression level of BRB gene. (a) Alignment of the reverse complemented sequences of genotypes from Hap1 and Hap2. The red box shows the location of the conserved motif (AACACACAC) in Hap2 and Hap1. The yellow box shows the TATA box. The genotype names underlined are those used later for RNA‐sequencing; (b) RT‐PCR analysis of transcript levels of the BRB gene in various sesame tissues. The band for Leaf in BRB row was rearranged from the initial gel images and this has been highlighted by a white space around the band; (c) β‐glucuronidase (gus) gene driven by the promoter of BRB was expressed in Arabidopsis main root and root hair but not in the root tip. The composite picture was generated by combining multiple pictures from different root sections; (d) expression levels of BRB gene in two contrasting genotypes for root trait based on RNA‐seq. G546 (Hap1) has dense root biomass while G259 (Hap2) has small root biomass; (e) qRT‐PCR analysis of the transcript levels of BRB in 25 genotypes (Hap1) and 25 genotypes (Hap2). Error bars indicate the SD (*** significant difference at P < 0.001); (f) Yeast one‐hybrid assay for the interactions between the transcription factor SIN_1023179 (SiMYB181) and the 5’‐UTR sequences of BRB (Hap1, Hap2 and synthetic sequences). The co‐transformed yeasts strain Y187 containing the bait and prey were cultivated on the synthetic dextrose minimal medium (SD) lacking leucine, tryptophan and histidine (SD/‐Leu/‐Trp/‐His) supplemented with 0 or 30 mm 3‐amino‐1,2,4‐triazole (3‐AT). N = negative control and P = positive control.

BRB is preferentially expressed in sesame root as compared to other organs such seed, leaf, flower and capsule (Figure 3b). Similarly, β‐glucuronidase (gus) gene driven by the promoter of BRB was expressed in Arabidopsis main root and root hair but was not expressed in the root tip (Figure 3c). Because of the specific position of the repeated motif in the 5’‐UTR of BRB, we hypothesized that it may affect the gene expression level among the accessions with Hap1 and Hap2. To test our hypothesis, we analysed the expression level of BRB using our previously released root transcriptome data from G546 (Hap1) and G259 (Hap2) (Su et al., 2019). We found that the expression level of BRB was ~ 45‐fold lower in G546 than in G259 (Figure 3d), suggesting that BRB is a negative regulator of root biomass in sesame. To further confirm this finding, we performed a quantitative reverse transcription PCR (qRT‐PCR) of root RNA samples from 25 Hap1 and 25 Hap2 accessions. The qRT‐PCR results validated that accessions with the single motif (Hap1) exhibited significantly lower level of BRB transcript as compared to those with the duplicated motif (Hap2) (P < 0.001) (Figure 3e) with similar observations in leaf samples (Data not shown).

The single motif AACACACAC facilitates the binding of an MYB transcription factor in the 5’‐UTR of BRB

Gene expression levels are known to be modulated by transcription factors (TF), binding to specific motifs in the promoter. We speculated that the duplication pattern of the motif AACACACAC in the 5’‐UTR of BRB may affect the binding of its regulator TF, leading to the differential expression level between Hap1 and Hap2. To examine this possibility, we first identified the TFs co‐expressed with BRB by performing a weighted gene co‐expression network analysis using 60 RNA‐seq datasets previously released by our group from root samples of two sesame genotypes at 10 time points of root growth (National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) accession number: SRP181800). BRB gene was co‐expressed with two WRKY and two MYB TFs, which represent the potential regulators (Table S9). Next, we performed yeast one‐hybrid assay using the cloned 5’‐UTR sequences from the Hap1 and Hap2. In addition, we also designed two different synthetic 5’‐UTR sequences by switching the alleles at the locus SNP5022842 (C/G) to assess the effect of this polymorphism on TF binding ability. The MYB gene SIN_1023179 (SiMYB181; Mmadi et al., 2017) was able to bind to the 5’‐UTR of BRB but specifically to the Hap1 sequences (Figure 3f). Also, the base change at the locus SNP5022842 did not affect the SiMYB181 binding. Altogether, our results suggest that the deletion of the motif AACACACAC in Hap1 accessions facilitates the binding of SiMYB181 to the 5’‐UTR of BRB, impeding the gene normal transcription. Therefore, SiMYB181 is a negative regulator of BRB. Nonetheless, it is still unclear how the duplicated motif affects the SiMYB181 binding to the 5’‐UTR of BRB in Hap2 genotypes.

BRB is a Lamiales‐specific gene, and the duplicated motif in the 5’‐UTR is under positive selection in sesame

The full‐length genomic DNA sequence of SIN_1025576 (BRB) gene consists of 417 bp long with a single exon. The ORF encodes a protein of 138 amino acids with a predicted molecular mass of 15.96 kDa which has no functional annotation. Searching for recognizable protein domains or signatures of conserved motifs in the protein sequence in publicly available databases did not result in any significant hits. We searched for BRB homologs within the plant kingdom through blastn of the NCBI non‐redundant nucleotide sequence collection database (E‐value < 10−10 and Identity > 70%). Two homologs (LOC105953437 and LOC105975223) were only detected in Erythranthe guttatus (Phrymaceae) and one (LOC111397704) in Olea europeae (Oleaceae), none being functionally annotated. Erythranthe guttatus and Olea europea are close relatives of sesame, and all belong to the Lamiales order. This suggests that BRB is an orphan gene restricted to the Lamiales order. We further confirmed this finding by successfully cloning the homologous of BRB using genomic DNA of three wild sesame species (Pedaliaceae): Ceratotheca sesamoides (Genbank number: MN336257), Sesamum radiatum (Genbank number: MN336258) and Sesamum alatum (Genbank number: MN336259) (Figure 4a). The motif AACACACAC was present as single or duplicated forms only in the 5’‐UTR of LOC105953437 from Erythranthe guttatus and the homologous BRB from all wild sesame species. The homologous BRB genes have the same structure, similar size and conserved nucleotide sequences (Figure S6; Figure S7).

Figure 4

Evolutionary analysis of BRB. (a) BRB was only detected in few plant species of the order Lamiales. Analysis of the 5’‐UTR of the BRB genes from different species showed the presence of no motif, single motif (Hap1) and duplicated motif (Hap2); (b) frequency of the presence of single motif (Hap1) and duplicated motif (Hap2) in sesame landraces and modern cultivars.

We further compared the frequency of the single/duplicated motif between landraces and modern cultivars to assess the impact of modern breeding at this locus. Intriguingly, we noticed that none of the modern cultivars harboured the single motif (Hap1) associated with the high root biomass trait while 20% of the landraces harboured the single motif (Figure 4b). This finding indicates that the duplicated motif (Hap2) was intensively selected by modern breeding; therefore, we hypothesized that Hap2 favours another desirable trait in the above‐ground part of sesame.

Overexpression of BRB in Arabidopsis modulates root and shoot traits and drought stress response

In order to confirm the function of BRB in another plant system, we created BRB‐overexpressing Arabidopsis mutants. There was not obvious difference in the main root length of wild‐type (WT) and BRB‐overexpressing lines (Figure 5a). However, the lateral root number, total lateral root length and root fresh weight were significantly reduced in the BRB‐overexpressing lines as compared to WT plants under control growth conditions (Figure 5a,b,c,d). This result demonstrates that BRB modulates root biomass and has potential applications in other plant systems. Since root traits are crucial for plant response to drought stress, we further tested whether overexpression of BRB will affect the drought response in Arabidopsis mutants. We observed that all mutant lines died just after 7 days of drought stress while no death was recorded in WT plants until 15 days (Figure 5e), indicating that BRB overexpression induces drought stress sensitivity.

Figure 5

Functional analysis of BRB gene in Arabidopsis thaliana. (a) Phenotypes of 14 days old wild‐type (WT) plants and BRB‐overexpressing lines (BRB‐OE1 and BRB‐OE2). BRB‐OE1 and BRB‐OE2 display reduced root lateral root number. Bar = 1 cm; (b–d) lateral root number, total lateral root length and root fresh weight were significantly lower in BRB‐OE1 and BRB‐OE2 compared to WT; (e) phenotypes of 30 days old WT, BRB‐OE1 and BRB‐OE2 plants after 7 days drought stress. Error bars indicate the SD (*,**,*** significant difference at P < 0.05, 0.01, 0.001).

Because we suspected that high expression of BRB favours another desirable trait in the above‐ground part of sesame, we evaluated the phenotypes of Arabidopsis mutants at the reproductive stage. Interestingly, we observed that under normal growth conditions, mutant lines displayed clustered silique and short internodes (Figure 6a). In addition, they had significantly higher plant height, silique number, branch number and seed yield per plant as compared to WT plants (Figure 6b,c,d,e,f). These observations confirm that BRB expression level not only modulates root traits but also impacts on the plant above‐ground part.

Figure 6

Overexpression of BRB affects shoot and yield‐related traits. (a) Phenotypes of the main stem and shoot of the BRB‐overexpressing (BRB‐OE) lines and wild‐type (WT) plants. Bar = 1 cm; (b–e) comparison of plant height, silique number, branch number and seed yield per plant between WT, BRB‐OE1 and BRB‐OE2 at maturation stage. Error bars indicate the SD (*,**significant difference at P < 0.05, 0.01).

BRB is modulates the auxin pathway

BRB gene has no functional annotation. In order to understand the functional pathway involving BRB, we compared transcriptome of WT and BRB‐overexpressing lines based on RNA‐sequencing. A total of 4685 and 3579 differentially expressed genes (DEG) were identified in WT_vs_BRB‐OE1 and WT_vs_BRB‐OE2, respectively (Tables S10 and S11). Gene Ontology enrichment analysis of the DEGs revealed various enriched biological pathways such as response to wounding, response to cold, response to abscisic acid (ABA), response to water deprivation, response to chitin and response to sucrose (Figure 7a,b). These enriched pathways hint that BRB modulates plant hormone signal transduction. Therefore, we examined the DEGs related to plant hormone signal transduction. We detected 68 DEGs related to various hormones such as abscisic acid, auxin/Indole‐3‐acetic acid (AUX/IAA), brassinosteroid, cytokinin, ethylene, gibberellic acid, jasmonic acid and salicylic acid (Figure 7c; Table S12). Notably, AUX/IAA‐related genes were the most altered by BRB overexpression in Arabidopsis, which were further validated by qRT‐PCR analysis (Figure S8). Hence, we deduce that BRB modulates the auxin pathway with consequences on root and shoot traits.

Figure 7

Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEG) between wild‐type (WT) plants and BRB‐overexpressing (BRB‐OE) lines. (a) GO analysis of DEG between WT and BRB‐OE1; (b) GO analysis of DEG between WT and BRB‐OE2; (c) proportion of phytohormone‐related genes detected within the DEGs. ABA: abscisic acid, AUX/IAA: auxin/Indole‐3‐acetic acid, BR: brassinosteroid, CK: cytokinin, ET: ethylene, GA: gibberellic acid, JA: jasmonic acid, and SA: salicylic acid.

Discussion

Deciphering the genetic architecture of root‐related traits in crops is crucial for improving productivity and abiotic stress tolerance (Liu et al., 2018; Svacina et al., 2014; Wang et al., 2015). Genome‐wide association study (GWAS) has been applied to unlock the genetic loci associated with root traits in crops such as wheat, maize, barley, rice (Abdel‐Ghani et al., 2019; Li et al., 2017; Li et al., 2019; Wang et al., 2019) and extensively in the model plant species Arabidopsis thaliana (Niu et al., 2018; Ogura et al., 2019; Ristova et al., 2018). However, most genes discovered for root traits in crops so far have only minor contributions (Uga and Yano, 2013). Also, key findings on root architecture in Arabidopsis are not directly transferable to crops. Studies on minor crops and non‐model plants have been largely neglected although they represent great sources of untapped alleles for improving root traits in major crops (Mayes et al., 2012). Herein, we studied the non‐model plant sesame and observed high variability and heritability for root traits, which is advantageous for GWAS implementation. Using a previously developed root traits phenotyping platform (Su et al., 2019) combined with high‐quality genotyping data (Wei et al. 2015), we successfully dissected the genetic architectures of seven root traits using GWAS. With the strong and positive correlations between root traits, we identified several QTLs with concurrent effects which will be beneficial for multi‐traits breeding in sesame as reported in jatropha (Sun et al., 2012). For instance, the QTL cluster QTL11 located on LG7 encompassed five different QTLs and could be applied for simultaneously improving root dry weight, root length, root volume, root surface area and shoot dry weight traits in sesame. Similarly, we identified several major loci (R 2 > 15%), indicating that breeding efforts for high root biomass can be simplified by pyramiding favourable alleles of major genes (Ye et al., 2019). Moreover, we reported for the first time 32 candidate genes for root traits in sesame. Homologs of some of these genes have been shown to control root traits in Arabidopsis. In qRL5.1 detected for root length, the candidate gene identified was SIN_1016062 and Arabidopsis loss‐of‐function mutant of its homologous gene (AT2G25010, MAIL1) developed short primary roots (Ühlken et al., 2014). Also, we identified the gene SIN_1024379 in qRSR9.1 for root–shoot ratio and its homolog in Arabidopsis (AT2G30520, RPT2) is involved in root phototropism (Sakai et al., 2000). The gene AT1G02730 (ATCSLD5) was shown to be involved in normal root growth (Bernal et al., 2007), and Yoo et al. (2012) later demonstrated its function in root hair development. In this study, we identified in qRV2.1 for root volume the candidate gene SIN_1005391, which is the homolog of ATCSLD5. Some candidate genes, however, have no clear functions related to root growth and development and require further functional investigations to uncover their exact contribution to root traits in sesame. One of these genes is SIN_1025576 named as ‘Big Root Biomass’ (BRB) in this study. BRB is an orphan gene restricted to the order Lamiales, and two independent association studies using different statistical models, populations and growth media demonstrated that a natural variation in the 5’‐UTR of BRB is strongly associated (with a major contribution R 2 > 30%) to root number and root dry weight traits in sesame. Furthermore, BRB overexpression in Arabidopsis confirmed not only its function in root biomass modulation but also highlighted the potential of this gene in non‐Lamiales species. Orphan genes do not possess a known protein domains (Khalturin et al., 2009); therefore, their functions remain largely unknown (Arendsee et al., 2014). However, it has been proven that they are important players in key agronomic traits. For example, Yadeta et al. (2014) identified the Brassicaceae‐specific gene EWR1 conferring resistance to vascular wilt pathogens which was also functional in Nicotiana benthamiana, a member of the Solanaceae family. Also, the grass‐specific gene Ms2 conferring male sterility (Ni et al., 2017) and the Pooideae orphan gene TaFROG which modulates resistance to diseases (Perochon et al., 2015) have been disclosed.

More interestingly, we have discovered a novel motif ‘AACACACAC’ in the 5’‐UTR of BRB and a polymorphism of the motif was found strongly associated with root biomass in sesame. Insertion–deletion of DNA sequences in the regulatory regions of key genes has been linked to different phenotypes. For example, the insertion of a terminal‐repeat retrotransposon in miniature (TRIM) element in the promoter of Ms2 is responsible for the anther‐specific Ms2 activation that confers male sterility in wheat (Xia et al., 2017). Ye et al. (2019) found a 8‐bp InDel in the 5’‐UTR of SlbHLH59 which regulates ascorbate biosynthesis in tomato. Likewise, the presence of tandem repeats of a 23‐bp forming a minisatellite‐like structure in the upstream regulatory region of the gene MYB10 directly controls anthocyanin biosynthesis in apple (Espley et al., 2009). In this study, we found that variations in the 5’‐UTR of BRB regulate its expression level through interaction with the negative regulator SiMYB181. In particular, the presence of a single motif facilitates the binding of SiMYB181, repressing BRB expression while the duplicated motif prevents the binding of SiMYB181, leading to a normal transcription of BRB. It is also possible that the duplicated motif promotes the binding of an inhibitor of SiMYB181. Similar observations were reported by Espley et al. (2009) and Ye et al. (2019) who showed that variations in the repeat numbers or insertion/deletion of DNA sequences affect the binding of a transcription factor to the causative gene and consequently alter its expression level.

The phytohormone auxin regulates many aspects of root traits, including root elongation, gravitropism and lateral root development (Davies, 1995). It is well known that when inhibitors of auxin transport are applied to root, it impedes lateral root initiation (Casimiro et al., 2001). In this study, we observed that BRB mainly controls the root lateral number then impacting the total root weight. We deduce that BRB contributes to a normal polar auxin transport in sesame root. Several genes of the auxin/Indole‐3‐acetic acid (Aux/IAA) family have been characterized as controlling root traits. For example, IAA1, IAA19, IAA17, IAA3, IAA7, IAA28, IAA14 control lateral root formation, root hair development and root gravitropic response (Fukaki et al., 2002; Leyser et al., 1996; Nagpal et al., 2000; Rogg et al., 2001; Tatematsu et al., 2004; Tian and Reed, 1999). In addition, functional characterization of these genes revealed overlapping, different or opposite effects on lateral root formation, indicating a complex regulation of auxin pathway to obtain a specific root phenotype. BRB has no signature of IAA proteins but significant changes were noticed in the auxin pathway in BRB‐overexpressing Arabidopsis lines as compared to wild‐type plants. Based on this observation, we conclude that BRB modulates the auxin pathway similar to reports of Ogura et al. (2019) who discovered an exocytosis factor, EXO70A3, as a modulator of the auxin pathway, affecting root system architecture and depth.

BRB overexpression in Arabidopsis showed altered shoot phenotype and improved seed yield under normal growth condition. Our experimental set‐up unfortunately does not allow us to record seed yield traits in our GWAS population in order to compare the two haplotypes. Nonetheless, based on our results in Arabidopsis mutants, we infer that BRB not only affects root traits but also impacts shoot and yield traits. It is well documented that auxins solely or in combination with other phytohormones shape plant architecture and regulate seed development, size and yield (Cao et al., 2020; Gallavotti, 2013; Shirley et al., 2019). Therefore, it is probable that high expression levels of BRB in above‐ground tissues are beneficial to the plant productivity and may explain why the Hap1 (allele conferring high expression level of BRB) has been intensively selected by modern breeding in sesame. In addition, we observed that Arabidopsis lines overexpressing BRB were highly sensitive to drought stress. Lateral root development is important for plant responses to abiotic stresses such as drought (Comas et al., 2013).

Functional characterization of orphan genes is challenging (Arendsee et al., 2014), especially in non‐model plant species such as sesame with limited tools for functional genomics and biotechnology. Hence, various unanswered questions remain on BRB. For example, after several attempts we were unable to clarify the subcellular localization of the product of this gene. Also, why SiMYB181 cannot bind to the repeated motif located in the 5’‐UTR of BRB (Hap2) and how BRB modulates the auxin pathway will need additional in‐depth investigations.

In summary, we report the genetic loci associated with various root traits in sesame and propose several candidate genes. BRB, a gene containing significant SNPs associated with root number and root dry weight, has been functionally investigated and our results confirm that it controls root biomass (Figure 8). Besides, BRB appears to affect differentially shoot traits, seed yield and drought stress response, depending on its expression levels. Further understanding of the way to control the expression levels of BRB either through its regulator SiMYB181 or the novel motif discovered in the 5’‐UTR region will provide crucial functional tool for developing new sesame lines with combined enhanced attributes, including dense root system, high seed yield, improved shoot traits and tolerance to abiotic stress. Beyond sesame, it is envisioned that BRB might have a significant contribution to the improvement of major crops.

Figure 8

Proposed hypothetical model. The presence of the duplicated motif AACACACAC in the 5’‐UTR of BRB gene impedes the binding of the negative regulator SiMYB181, leading to a normal transcription of BRB. In contrast, the presence of the single motif AACACACAC in the 5’‐UTR of BRB gene facilitates the binding of SiMYB181, leading to a weak transcription of BRB. High and low expression levels of BRB differentially affect the auxin pathway to reduce and increase root biomass, respectively, through a yet undetermined mechanism.

Methods

Plant materials

We extracted 327 sesame (Sesamum indicum L.) accessions comprising landraces and modern cultivars (Table S1) from our association panel established at the Oil Crops Research Institute Chinese Academy of Agricultural Sciences (OCRI‐CAAS) (Wei et al., 2015). Based on the root phenotyping data, two contrasting genotypes for root biomass (G430 with a big root biomass and G441 with a small root biomass) were selected to develop 197 F2:3 population. In addition, three wild sesame including Sesamum alatum, Sesamum radiatum and Ceratotheca sesamoides were obtained from the Sesame Germplasm Resource Preserving Center of OCRI‐CAAS.

Growth conditions of the sesame association panel

The association panel was grown hydroponically in blue plastic basins (34 × 26 × 12 cm, length × width × height) containing modified half‐strength Hoagland’s solution (Hoagland and Arnon, 1950) following technical descriptions of Su et al. (2019). Twenty‐four seedlings of six accessions (four replicated seedlings for each accession) were grown in each basin. The experiment was led under controlled environmental conditions with the temperature and relative humidity kept at 35 °C and 60%, respectively, and under long‐day condition (16/8 h day/night). A completely randomized design was employed with three replications, resulting in a total of ~ 3600 seedlings for each experiment. Plants were grown for a week in ¼ strength Hoagland nutrient solution and subsequently transferred into ½ strength nutrient solution for three additional weeks. The nutrient solution was renewed once a week, and the whole experiment was replicated over three different time periods during 2016–2017. At the end of each experiment, the whole root was delicately separated from the shoot, and both samples were placed into paper bags for later use.

Growth conditions of the F2:3 population

Sterile seeds of the 197 F2:3 population and the two parents were grown in pots (25 cm diameter and 45 cm depth) filled with 6 kg of experimental soil composed of 1.5 kg vermiculite, 2 kg soil, 2 kg sand and 0.5 kg nutritive soil as detailed by Su et al. (2019). The experiment was led under controlled environmental conditions with the temperature and relative humidity kept at 35 °C and 60%, respectively, and under long‐day condition (16/8 h day/night) at a PPFD of 200 μmol/m2/s. Plants were watered regularly to keep normal growth. The experiment was laid out in a completely randomized design with four replications and single plant per pot. The root of 4‐week‐old plants was delicately separated from the shoot, thoroughly washed with tap water, and both samples were put into paper bags for later use.

Root image scanning and data acquisition

After harvest, root samples were quickly scanned using a desktop scanner (EPSON Perfection V800 Photo, EPSON, Carson, CA, USA). Intact root systems were spread to minimize overlaps in a tray (A4 size) with a transparent glass bottom which was placed on the scanner. Next, the root systems were scanned to produce digital images in uncompressed tagged image file (TIF) format (400 dpi). Root images were analysed using the WinRHIZO Pro software (Regent Instruments Inc., Quebec, QC, Canada). Four root data including the total root length (RL, cm), root surface area (SA, cm2), root volume (RV, cm3) and root number (RN) were extracted. Shoots and roots were oven dried at 80 °C for 72 h and weighed to determine shoot dry weight (SDW, g) and root dry weight (RDW, g). Root–shoot ratio (RSR) data were then estimated. Broad‐sense heritability (H 2) was estimated following descriptions of Dossa et al. (2019a).

GWAS implementation

To perform the genome‐wide association analysis (GWAS), the best linear unbiased prediction (BLUP) value of each phenotypic trait was computed using mixed linear model (MLM) in R 3.1 (R Core Team, 2015). Next, we extracted from the sesame HapMap project (www.ncgr.ac.cn/SesameHapMap/), one million high‐quality SNPs with a minor allele frequency ≥ 0.05 and a good coverage of the whole genome (Wei et al., 2015). GWAS was conducted based on the mixed model with the SNP data, the kinship matrix, population structure (first three PCA components) using the EMMAX package (Kang et al., 2010).

Suggestive P value threshold was set at P = 1 × 10−5 to identify significant SNPs because root traits are complex polygenic and only small effects of the individual underlying loci were expected. A similar suggestive P value was previously used for dissecting complex traits in sesame (Li et al., 2018; Dossa et al., 2019a). For clustered significant SNPs, the peak SNP was defined as the one having the lowest P value in the linkage disequilibrium (LD) region. The LD decay in sesame has been estimated at ~ 88 kb (Wei et al., 2015; Dossa et al., 2019a), and LD windows upstream and downstream of the peak SNPs were defined as quantitative trait loci (QTL). To estimate the phenotypic variance explained (PVE) of each peak SNP, the value R 2 derived from linear regressions were calculated in the R package (R Core Team, 2015).

BRB gene sequencing and candidate gene association analysis

The full coding sequence (417 bp), 292 bp upstream of the translation start codon and 180 bp downstream of the stop codon of the gene BRB (SIN_1025576) were sequenced in the 197 F2:3 population and the two parents using the primers A1 (Table S13). With the DNAMAN software (Lynnon Biosoft, Quebec, Canada), we assembled the sequences and the polymorphisms including SNPs and insertion–deletion (MAF ≥ 0.05) were detected using the DNASP software v6.11.01 (Rozas et al., 2017). Next, the association between the polymorphisms and two root traits (RDW and RN) was tested using TASSEL5.0 (Bradbury et al., 2007) based on the general linear model. Only significant variants with P ≤ 0.001 were retained.

Gene co‐expression analysis

A total of 60 root RNA‐seq datasets generated from two sesame genotypes at 10 time points (0, 4, 8, 12, 16, 20, 24, 40, 44, 48 h beginning 15 days after the initiation of flowering) of root growth under control condition released by our group (NCBI SRA accession number: SRP181800) (Dossa et al., 2019a; Wang et al., 2020) were employed for gene co‐expression analysis according to descriptions of Lv et al. (2019). The analysis was conducted with the WGCNA package version: 1.61 (Langfelder and Horvath, 2008) based on the normalized log2‐transformed FPKM matrix. The co‐expressed module containing the BRB gene (SIN_1025576) was extracted, and GO/KEGG enrichment analyses of the module were performed using clusterProfiler version 3.8.

Quantitative real‐time PCR

The qRT‐PCR analysis of target genes was performed according to descriptions by Dossa et al. (2019a) using the ChamQ SYBR qPCR Master Mix (Vazyme Biotec, Nanjing China) on a Light Cycler 480 II (Roche, Basel, Switzerland). The relative expression levels of target genes were normalized to the expression levels of two endogenous genes: SiActin 7 and SiHistone 3 for sesame and one endogenous gene AtActin 2 for Arabidopsis using the primers A2‐A32, (Table S13). Data are presented as relative transcript level based on the 2−ΔΔCt method (Livak and Schmittgen, 2001) with three technical and three biological replicates. For RT‐PCR, high‐quality RNA was extracted from various tissues using the EASYspin Plus kit (Aidlab, Beijing, China) according to the manufacturer’s instructions and reverse transcribed using the Superscript III reverse transcription kit (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions. The primer A2 (Table S13) was used to amplify BRB gene, and PCR products were visualized on 2% agarose gel.

Yeast one‐hybrid assay

The 5’‐UTR sequences from Hap1 and Hap2 containing a single and repeated AACACACAC motif, respectively, were isolated by PCR using the primer A33 (Table S13). Further, we synthesized two mutated synthetic 5’‐UTR sequences by switching the alleles at the locus SNP5022842 (C/G). The four sequences were cloned into the bait vector pHIS2 (Clontech, Palo Alto, California, USA) between EcoRI/SmaI and EcoRI/SacI sites. Next, the complete CDS of the transcription factors (TF) co‐expressed with BRB were amplified by the primers A34‐A37 (Table S13) and cloned into the prey vector pGADT7‐Rec2 (Clontech, Palo Alto, California, USA) using EcoRI/XhoI sites. Then, the co‐transformed yeasts strain Y187 containing the bait and prey were cultivated on the SD/‐Leu/‐Trp/‐His selective media supplemented with 0 or 3 mm 3‐amino‐1,2,4‐triazole (3‐AT) for 3 days according to the instructions for the MatchmakerTM Gold Yeast One Hybrid System (Clontech, Palo Alto, California, USA). Yeasts co‐transformed with pGADT7‐Rec2 −53 (the pGADT7‐Rec2‐53 plasmid containing a murine p53 and GAL4 AD domain fusion gene) and p53HIS2 (the pHIS2 harbouring three tandem copies of p53 recognized DNA consensus motifs) was used as positive control. The negative control was pGADT7‐Rec2‐TF and p53HIS2 co‐transformation. The interaction between prey and bait was observed according to the growth of the yeast transformants in a series of 10‐fold dilution.

Arabidopsis transgenics experiment

To functionally characterize the BRB gene in Arabidopsis thaliana, we isolated the coding sequence by PCR from G546 (Hap1) using the primer A38 (Table S13). Gene cloning and Arabidopsis transformation were performed as described by Dossa et al. (2019a) using the pCAMBIA 1301s vector. Positive T1 plants were screened on Murashige and Skoog (MS) medium containing 1% agar, and 1% sucrose and 50 μg/mL hygromicin and further confirmed by PCR, RT‐PCR and β‐glucuronidase (GUS) staining. Three independent T3 transgenic homozygous lines were used for the evaluation of root and yield traits and RNA‐seq analysis. We also generated BRB‐promoter:GUS expression cassette, which contains the GUS reporter gene under the control of the promoter‐containing DNA fragment from BRB (292 bp) in pCAMBIA‐1381z binary vector. The resulting recombinant vector was then introduced into Agrobacterium tumefaciens strain LBA4404 and transferred into Arabidopsis. Histochemical localization of GUS activity in the generated transgenic Arabidopsis lines was performed by incubating whole seedlings with X‐Gluc (5‐bromo‐4‐chloro‐3‐indolyl‐β‐d‐glucuronide) as described by Jefferson et al. (1987).

Evaluation of the transgenic lines

Seeds of wild‐type (WT) and three T3 lines were germinated on solid MS agar medium. The seeds were stratified for 2 days in the dark at 4 °C and then transferred to growth chamber at 22 °C under long‐day conditions (16 h light/8 h dark) with light intensity of 120–150 µmol/m2 sec and relative humidity of 70%. Thirty seedlings (10 days old) of WT and BRB‐overexpressing lines were transferred into new MS plates. Plates were placed vertically, and after seven days, lateral root length, lateral root number and estimated total lateral root length were obtained by summing up lengths of all lateral roots. In addition, root fresh weight of bulked seedling was recorded. At the reproductive stage, drought stress was applied by withholding water supply for 15 days and the number of survived plants out of 50 investigated plants was recorded. In parallel, 50 plants of each line were grown under well‐watered conditions, and plant height, branch number, silique number per plant and seed yield per plant data were recorded at the maturity stage. The entire experiment was repeated two more times with four replicates in each experiment for statistical analysis.

Transcriptome sequencing and analysis

Total RNAs were extracted from whole plants of two independent Arabidopsis mutant lines and WT plants in two biological replicates using the EASYspin Plus kit (Aidlab) according to the manufacturer’s instructions and reverse transcribed using the Superscript III reverse transcription kit (Invitrogen) according to the manufacturer’s instructions. The cDNA libraries were pair‐end sequenced on an Illumina Hiseq 4000 platform according to the methods described by Dossa et al. (2017b). The raw data were processed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to filter out adapters and low‐quality sequences. Thereafter, the clean reads were mapped to the reference genome of Arabidopsis (TAIR10) using HISAT (Kim et al., 2015). The gene expression level for each sample expressed as fragments per kilobase of transcript per million fragments mapped (FPKM) was estimated employing the RSEM package v1.3.0 (Li and Dewey, 2011). Comparison of the gene expression levels between samples helped to identify the differentially expressed genes with the following parameters: Fold change ≥ 2, Probability ≥ 0.8, false discovery rate‐adjusted P value < 0.05 (Tarazona et al., 2011). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for the DEGs were performed using clusterProfiler version 3.8.

Statistical analysis

All the data were analysed with the R software (www.r‐project.org). One‐way analysis of variance was performed by comparing each transgenic line to the wild‐type plants. This was followed by Tukey HSD test for mean comparison.

Conflict of interest

They declare no conflict of interest.

Author contribution

Komivi Dossa, Xiurong Zhang and Jun You conceived and designed this project. Komivi Dossa, Rong Zhou, Donghua Li, Aili Liu, Lu Qin, Marie A. Mmadi, Ruqi Su, Yujuan Zhang, Jianqiang Wang and Yuan Gao conducted the experiments and collected data. Komivi Dossa performed data analysis and drafted the manuscript. Xiurong Zhang and Jun You supervised the study, provided funding and technical support, and revised the drafts of the manuscript. All authors have read and approved the final version of this manuscript.

Figure S1 Frequency distribution of root traits in 327 Sesamum indicum accessions.

Figure S2 Pearson correlation analysis of the different root‐related traits. Total root length (RL), root surface area (SA), root volume (RV), root number (RN), shoot dry weight (SDW), root dry weight (RDW) and root–shoot ratio (RSR).

Figure S3 Manhattan plots of genome‐wide association studies using the mixed model for seven root‐related traits. The top significant trait‐associated loci were marked with a red star.

Figure S4 Cumulative contribution (Phenotypic variance explained = PVE) of the top 5 significant SNPs associated with the different root traits. Total root length (RL), root surface area (SA), root volume (RV), root number (RN), shoot dry weight (SDW), root dry weight (RDW) and root–shoot ratio (RSR).

Figure S5 Alignment of the 5’‐UTR reverse complemented sequences from 25 sesame genotypes of Hap1 group and 25 sesame genotypes from Hap2 group. The blue box shows the conserved motif (AACACACAC) between Hap2 and Hap1 groups, while the red box shows the duplicated motif only observed in Hap2 group.

Figure S6 Gene structure of BRB from various species.

Figure S7 Alignment of the coding sequences of the BRB gene from various species.

Figure S8 qRT‐PCR analysis of the auxin‐related genes differentially expressed between wild‐type plants and BRB‐overexpressing Arabidopsis lines. Data are represented as mean ± SD of three biological replicates.

Click here for additional data file.

Table S1 Full list of the 327 accessions used in this study, their origin and their breeding status.

Table S2 Summary of descriptive statistics of the seven root traits investigated in this study. CV: coefficient of variation; H2: broad‐sense heritability; SD: standard deviation, Min: minimum value; Max: maximum value.

Table S3 SNPs significantly associated with root traits in Sesamum indicum. PVE: phenotypic variance explained, LG: linkage group. The QTL names in red are those significantly associated with various root traits.

Table S4 List of genes located within QTLs detected for root traits through GWAS in sesame.

Table S5 List of genes containing significant SNPs in coding or regulatory regions and their functional annotation.

Table S6 Candidate genes identified for root traits in sesame and their homologs in Arabidopsis thaliana.

Table S7 The 11 significant SNPs located in BRB (SIN_1025576).

Table S8 Association analysis of root number (RN) and root dry weight (RDW) traits in 197 F2:3 sesame population. The markers highlighted in red are those significantly associated with RN an RDW.

Table S9 Genes co‐expressed with BRB (SIN_1025576) and their annotation. The four transcription factors were highlighted in red.

Table S10 List of the differentially expressed genes identified between BRB‐overexpressing line 1 (BRB‐OE1) and wild‐type (WT) plants.

Table S11 List of the differentially expressed genes identified between BRB‐overexpressing line 2 (BRB‐OE2) and wild‐type (WT) plants

Table S12 Differentially expressed genes related to phytohormones between BRB‐overexpressing lines (BRB‐OE) and wild‐type (WT) plants.

Table S13 Sequences of the different primers used in this study.

Click here for additional data file.