Response to early drought stress and identification of QTLs controlling biomass production under drought in pearl millet

Pearl millet plays a major role in food security in arid and semi-arid areas of Africa and India. However, it lags behind the other cereal crops in terms of genetic improvement. The recent sequencing of its genome opens the way to the use of modern genomic tools for breeding. Our study aimed at identifying genetic components involved in early drought stress tolerance as a first step toward the development of improved pearl millet varieties or hybrids. A panel of 188 inbred lines from West Africa was phenotyped under early drought stress and well-irrigated conditions. We found a strong impact of drought stress on yield components. This impact was variable between inbred lines. We then performed an association analysis with a total of 392,493 SNPs identified using Genotyping-by-Sequencing (GBS). Correcting for genetic relatedness, genome wide association study identified QTLs for biomass production in early drought stress conditions and for stay-green trait. In particular, genes involved in the sirohaem and wax biosynthesis pathways were found to co-locate with two of these QTLs. Our results might contribute to breed pearl millet lines with improved yield under drought stress.


Introduction
Pearl millet (Pennisetum glaucum (L.) R. Br., syn. Cenchrus americanus (L.) Morrone) was domesticated in the central Sahel region of Africa [1]. It is adapted to dry and hot climates and it plays a major role in food security in arid and semi-arid areas of Africa and India. However, there is still a wide scope to increase its adaptation to environmental constraints limiting its production. Pearl millet is the staple crop of more than 90 million farmers. Grains are glutenfree and rich in proteins and essential micronutrients such as iron or zinc. It is consumed directly as couscous or porridge and also used to make flour to produce flat bread or to be introduced in bread production to reduce the imports of wheat in some African countries. Moreover, the aerial biomass is an important source of fodder for animals. However, pearl millet lags behind the other cereal crops in terms of genetic improvement and its average yields are still low. The recent sequencing of its genome [2] opens the way to tap on its large genetic diversity to breed varieties and hybrids adapted to current and predicted future climatic constraints [3].
Drought is one of the major factors limiting pearl millet production. In sub-Saharan Africa, pearl millet is mostly grown in areas characterized by low rainfall and sandy soils with very little organic matter that have therefore low water retention ability. The climate is characterized by a long dry season and a short rainy season where most of the rain-fed agriculture is concentrated [4]. Pearl millet is usually sown after or just before the first rain of the rainy season. Because of the rain pattern, pearl millet can face drought stress at early stages if the first rains of the season are distant from each other. Another important drought stress period is the grain filling stage. Current climate models predict that the rain pattern in sub-Saharan Africa will be more variable from one year to the other and will involve more extreme events [4]. Episodes of drought are therefore expected to be more frequent during the rainy season. Altogether, climate change is expected to lead to a reduction in pearl millet yield in the area [4]. While terminal drought stress was the focus of several studies in pearl millet [5][6][7][8] little is known about the impact of drought episodes during the vegetative phase.
The aim of this study was to characterize the impact of early drought episodes on pearl millet in field conditions and to identify genetic components contributing to tolerance to this stress.

Plant material and trait measurements
A panel of 188 pearl millet (Pennisetum glaucum (L.) R. Br.) inbred lines developed at ICRISAT (Niger) from landrace and improved open-pollinated cultivars of West African origin was used in this study together with inbred line ICMB88004 as a control (S1 Table). Some of these inbred lines were previously characterized for grain quality [9].
Field trials were performed at the CNRA station (Centre National de Recherche Agronomique) of the Institut Sénégalais des Recherches Agricoles (ISRA) in Bambey, Senegal (14.42˚N, 16.28˚W) during the dry season in 2015 (planting: March 24 th , harvesting: June 16 th ) and 2016 (planting: March 10 th , harvesting: June 14 th ) on two adjacent fields. Field trials were conducted in collaboration with and with the permission of the ISRA. Trials did not involve endangered or protected species. Soils are deep sandy soil with low levels of clay and silt (12%) and organic matter (0.4%). Clay and silt content increase with soil depth from 10.2% in the 0 to 0.2 m layer to 13.3% in the 0.8 to 1.2 m layer. Trials were set up using an incomplete randomized blocks design (S1 and S2 Figs). A total of twenty plants were grown per variety in 2 rows of 3 m long with 0.30 m between plants and 0.8 m between rows. Fifteen days after sowing, a single plant per planting hole was conserved, so ten plants per variety for a given row. Fertilization (NPK) followed standard recommendation i.e. 150 kg ha -1 of NPK (15-15-15) after sowing. Urea was applied at 100 kg ha -1 at two stages, 50 kg ha -1 after thinning and 50 kg ha -1 30 days after sowing.
The flowering date corresponding to 50% of plants flowering (in days after sowing) was recorded. Height of the plant (in cm), panicle and stalk length (in cm), total number of leaves, aerial biomass (kg.plant -1 ), as well as stay-green trait (number of green leaves/total number of leaves) were measured on six plants per plot at maturity. Shoots and panicles were harvested and used to quantify total grain weight (g/plant), number of grains per plant, weight of 1000 grains (g) and aerial biomass (g dry mass/plant). These data were used to compute a harvest index (grain weight/biomass weight).
Plants were grown during the dry season (February to June) to fully control the amount of water provided (no rain during the experiments). Irrigation was performed twice per week with 30 mm water per irrigation. Early drought stress was applied by stopping irrigation three weeks after germination for four weeks. Irrigation was then resumed until the end of the growth cycle. Field dry-down was followed by measuring volumetric soil moisture to evaluate the fraction of transpirable soil water (FTSW) using Diviner probes (Sentek Pty Ltd). Canopy temperature was measured on two plants per plot 30 days after the start of the dry down period in both well-irrigated and drought-stress treatments using an Infrared thermometer (Quicktemp).
The spatial trend in the incomplete block design was evaluated using control plants (inbred line ICMB88004) and the phenotype of each line was corrected using the SpATS package in R (S3 Fig; [10]). Corrected data are provided as S2 Table. DNA isolation and genotyping by sequencing Genomic DNA was isolated from three-week old leaves as previously described [11]. A single plant was usually used per accession but for 68 inbred lines, 2 sister plants were used to control the procedure. Whole-genome genotyping was carried out using Genotyping-By-Sequencing (GBS) at the Genomic Diversity Foundation at Cornell University (Ithaca, USA). GBS libraries were prepared using restriction enzyme ApeKI [12]. Digested DNAs were ligated to barcoded adapter pairs. 96-plex libraries were sequenced using a HiSeq2500 and a NextSeq500 sequencing system (Illumina). Raw sequence data are available at Genbank under accession PRJNA492967. Adapter sequences were removed using the cutadapt software and low quality reads were filtered out with the pearl script Filter_Fastq_On_Mean_Quality.pl from South-GreenPlatform (Minimum mean quality allowed for a read = 30, Minimum length allowed for a read = 35). The reads were mapped on the pearl millet reference genome sequence using BWA [13,14]. Picard-tools-1.119 and Genome Analysis ToolKit (GATK-3.6 algorithm Uni-fiedGenotyper) softwares were used to create the vcf file. SNPs were filtered out based on missing percentage (<50% missing), homozygosity (inbreeding F > 0.5) and minor allele frequency (MAF � 5%).

Population structure analysis
A neighbor joining (NJ) phylogenetic tree was generated with TASSEL v5.2.39. Population structure was assessed with a random subset, with one SNP/100kb using the program sNMF. Five runs were performed for a given number of group (K) and the values with the smallest Cross-Entropy for each K were selected to generate the structure graph.

Genome wide association study (GWAS)
Marker-trait associations were established using pearl millet inbred lines phenotyped for 11 agromorphological traits under two conditions (early drought stress and well-watered) and in two years (2015 and 2016). For each year, we performed GWAS with a mixed model correcting for genetic relatedness (kinship matrix) using the R package GAPIT (Genome Association and Prediction Integrated Tool; [15,16]). The p-values obtained for both years were combined using a Fisher combining probability method [17]. We considered a False Discovery Rate (FDR, threshold = 0.1) in determining the significant SNPs.

Phenotypic variation under well-watered and early drought stress conditions
In order to evaluate pearl millet response to early drought stress, a panel of 188 inbred lines was characterized on well-watered (WW) and early drought stress (DS) conditions on two successive years in Bambey, Senegal. In the well-watered treatment, the fraction of transpirable soil water (FTSW) was generally above 40% along the soil profile (0-120 cm) across the experiment in both years, indicating that water was not limiting for plant growth (S4 Fig; [18]). In early drought stress treatments, FTSW was below 40% in 0-60 cm soil profiles at 40 days after sowing (DAS) in 2015 and at 49 days in 2016. This measure indicated efficient field dry-down and imposition of water limiting conditions at vegetative stage (S4 Fig). In 2015, water remained limiting for plant growth (FTSW below 40%) along the soil profile (0-120 cm) until around 110 DAS, extending the water-limiting period to the reproductive stage. In 2016, irrigation of the field at 49 DAS allowed increase in FTSW to around 40% at 55 DAS between 60-120 cm, although short periods of water limiting conditions appeared at reproductive stage (around 75 and 85 DAS). In both 2015 and 2016, canopy temperature measured 30 days after the start of the dry-down was significantly increased in the drought stress treatment as compared to the well-watered treatment (ANOVA, p < 0.001), indicating that plants were indeed subjected to drought stress (S5 Fig).
In order to evaluate the impact of early drought stress, a total of 11 agromorphological traits were measured at the seed maturation stage (Table 1). For most of the traits, a large range of variation was detected, with the coefficients of variation (CV) varying from 8.2% for staygreen under well-watered conditions to 98.6% for total grain weight under drought stress conditions. Most of the traits showed a normal distribution (S6 Fig) apart from flowering date that was flatter than a normal distribution (negative kurtosis). For each trait, each year and all conditions, variance analyses showed a strong significant genetic effect on phenotypic variation ( Table 2). A high correlation was found between results obtained in the two successive years both for well watered (WW) and drought stress (DS) conditions indicating that the response was comparable in both years (Fig 1). Analyses of relation between traits showed that some traits were highly correlated (Fig 2 and S2 Table) such as total grain number and total grain weight. Correlated traits under drought stress conditions were generally also correlated under well-watered conditions (correlation of correlation coefficients, r > 0.8, p<10 −12 , Fig 2C and  2D and S2 Table).
In order to evaluate drought effects on agromorphological traits, values measured under drought stress were divided by values measured under well-watered conditions for each trait in both years (Table 3). Drought stress led to a strong and very significant reduction in plant height, above ground biomass, total grain number and weight on both years (Table 3). Panicle length, stalk length and number of leaves were moderately reduced under drought stress. On the contrary, stay green increased slightly but significantly under drought stress on both years. However, the very low harvest index and the failure of some lines to set seeds suggest that the sink demand from grain was low and therefore the observed effect might not be related to a functional stay-green character, i.e. transition from carbon capture to nitrogen remobilization [19]. Flowering date and weight of 1000 grains were unevenly affected on the two years, showing either small or no variation under early drought stress. Harvest index (HI) was not significantly affected by drought on both years. Hence, early drought stress led to a reduction both in biomass and number of grains leading to reduced yield (but not in HI and grain weight). No change in individual grain weight was observed.
The loss in grain production under drought stress conditions was very variable among lines and years. Some lines did not produce grain under drought stress even if their panicle length and aboveground biomass were within the same range as yielding plants. This suggests that drought stress applied at the juvenile phase impacts seed production at later stages. On the  Response to early drought stress in pearl millet ratio for aerial biomass (BSA) and total grain weight (PGR) in both years. On the other hand, ICML-IS 11066 had low PGR ratios in both years. Altogether, our data indicate that response to early drought stress was variable in our panel but conserved across lines between the two experiments. This suggests that this population could be used to identify QTLs/genes involved in early drought tolerance through phenotype/genotype association analyses and could be a source of potential donor lines for drought tolerance.

Genotypes and population structure
In order to perform genetic association, we generated GBS data on the 188 lines and obtained a total of 3,168,971 unfiltered single nucleotide polymorphisms (SNPs). Filtering on quality led to 392,493 SNPs. The average SNP density was 2.5 per 10 kb.
SNPs were first used to assess the population structure in the panel. The cross-validation error reached a minimum for K = 1 suggesting only a single genetic group (S7 Fig). The unweighted neighbor-joining (NJ) tree constructed to illustrate the phylogenetic structure of the panel (Fig 3) displayed also a weak structuration signal. Altogether, this indicates that our panel has a negligible genetic structure.

Genome wide association study (GWAS)
Genetic association (GWAS) was conducted on 175 inbred lines presenting both genetic and phenotypic data. A total of thirteen plants were removed because their genotypic data did not pass filter quality tests. GWAS was performed on all phenotypic traits measured in two years (2015 and 2016) in WW and DS conditions. The most significant associations (p−value < 10 −5 ) detected in the GWAS study are listed in Table 4. The FDR for the associations ranged from 0.0001 to 1, with 10 associations having a FDR < 0.1.
In WW conditions, we identified four marker-trait associations (MTAs) for the stay-green character on chromosome 6 (Fig 4 and Table 4). These four MTAs corresponded to two different polymorphisms located at the same position. The stay-green MTA explained 12% of phenotypic variation in the panel in 2015 and a difference of 25% in phenotype was observed between the two alleles (Fig 5). The corresponding polymorphisms were mapped at position 111262313 on chromosome 6 of the pearl millet reference genome [2]. It corresponded to a peak of SNPs/stay green association (Fig 5). Interestingly, this locus falls within a predicted gene (Pgl_GLEAN_10-013220) encoding a potential uroporphyrin-III C-methyltransferase (UPM), an enzyme involved in sirohaem and cobalamin biosynthetic pathway. Sirohaem is an essential cofactor involved in inorganic sulfur (S) and nitrogen (N) assimilation in plants [20]. No significant association was found in DS conditions but as reported above it might be due to the failure of some lines to set seeds leading to reduced sink demand from grain preventing leaf senescence. Six significant associations were detected for biomass on chromosomes 2, 3 and 6 in DS conditions (Fig 4 and Table 4). However, one inbred line had high BSA under DS (Fig 6) and the associations were lost when this line was removed from the analysis. Three of these associated SNPs were very closely located at position 216,179,589, 216,179,591 and 216,179,592 on chromosome 3. These markers explained 19% of the phenotypic variation in the panel. The 3 SNPs were located between two predicted genes Pgl_GLEAN_10034145 and Pgl_GLEAN_10034146 encoding a putative 3-ketoacyl-CoA synthase and an unknown protein respectively. 3-ketoacyl-CoA synthases are involved in elongation of C24 fatty acids, an essential condensing step during wax and suberin biosynthesis [21]. Comparison of the genome sequence of pearl millet and other cereals revealed an expansion of gene families involved in wax and suberin biosynthesis that was proposed to have contributed to pearl millet adaptation to drought and heat stress [2]. Accordingly, changes in wax and/or suberin biosynthesis linked to a 3-ketoacyl-CoA synthase encoding gene might be responsible for quantitative changes in biomass under drought stress conditions as observed in our field trials.
The association on chromosome 2 was mapped at position 17,262,004. It explained 14% of the phenotypic variation in the panel (Fig 6). The SNP is located in the predicted gene Pgl_GLEAN_10005405 that encodes a putative chloroplastic threonine dehydratase. This enzyme catalyzes the formation of alpha-ketobutyrate from threonine. This catabolic reaction might change the amino acid metabolism in response to stress and facilitate photorespiration. Photorespiration is known to sustain CO 2 fixation capacity in conditions of water stress and to improve drought tolerance in cereals [22]. Two SNPs on chromosome 6 were significantly associated with biomass production in early drought stress conditions. The first one was mapped at position 64,437,341, between two predicted genes Pgl_GLEAN_10037360 and Pgl_GLEAN_10037359 encoding an unknown protein and a SPA1-related 3 protein homolog. The SPA proteins contain serine/threonine kinase-like and WD40 protein domains and are involved in signal transduction during photomorphogenesis [23]. It has been shown that SPA1 modulates MYC2-mediated ABA and JA responses [24] that are known regulators of plant response to stresses including drought. This association explained 14% of the phenotypic variation in the panel (Fig 6). The second SNP associated with biomass production under drought was located at position 157,567,763 between two predicted genes Pgl_GLEAN_10036945 and Pgl_GLEAN_10036946 encoding an unknown protein and a putative E3 SUMO-protein ligase respectively. E3 SUMO-protein ligases regulate protein sumoylation and have been involved in tolerance to a number of abiotic stresses including drought [24]. In Arabidopsis, E3 SUMO-protein ligase SIZ1 contributes to the accumulation of SUMO-protein conjugates induced by drought stress [25]. Moreover, SIZ1 was demonstrated to regulate growth and drought stress response [25] as well as abscisic acid signaling [26].

Discussion
As little is known about the impact of early drought episodes in pearl millet, we analyzed the response of a panel of 188 pearl millet inbred lines to drought stress during the juvenile phase for two years in field conditions in Senegal. Drought stress was applied from 21 days after germination (DAG) to 49 DAG. On drawback of such approach is that while the stress was applied at the same fixed time for all lines, we have a large diversity of flowering time in our panel and all plants did not therefore face stress at the exact same developmental stage. However, the earliest flowering genotypes in our panel flowered at 47 or 39 DAG in well-watered environment and at 46 or 48 DAG in drought environment depending on the year. Beginning of flowering therefore corresponded almost to the end of drought stress for the early flowering lines. The mean flowering time in our panel was 66 and 69 for the two years respectively for drought stress treatment, so between 17 to 20 days after the end of the stress. Consequently, the majority of the lines flowered after drought stress. Accordingly, we did not find any significant relation between flowering time and PGR (grain yield) or panicle length under stress in both years. One strategy to avoid such issue would have been to work with groups of lines having similar flowering time but this would have limited the number of lines available for analysis and therefore the potential genetic diversity usable for gene discovery.
We found that early drought stress led to a reduction in both grain and biomass production. Limited change in grain weight (as would be expected for a terminal drought stress) was observed only one year out of two suggesting that most pearl millet lines were able to adapt their biomass to water availability in order to produce a limited number of viable seeds. Yield loss in early drought stress conditions was very variable among lines. Some lines did not produce grain under drought stress even if their panicle length and aboveground biomass were within the same range as yielding plants. On the other hand, some lines did perform better in drought stress conditions. Nevertheless, our results clearly show that drought episodes during the vegetative phase of the plants can have a dramatic impact on yield even if water is available during the end of the cycle. By reducing water availability at this early stage, plants were facing drought stress during the juvenile phase, when the young panicle develops, deep inside the meristem. How this impacts seed formation later in the cycle when water supply remains to be deciphered. Response to early drought stress in pearl millet We exploited the phenotypic diversity observed in our panel to perform GWAS to uncover the genetic bases of tolerance to early drought stress in pearl millet. We identified 10 SNPs associated with stay-green and biomass phenotypes corresponding to five potential QTLs. As pearl millet biomass is increasingly used as fodder, biomass quantity and quality (of which stay-green is an important character) are becoming important targets for breeding. Moreover, changing temperature and photoperiod associated to climate change are known to impact senescence [27]. Regulation of stay-green is therefore important to design varieties better adapted to future climate in West Africa. Stay-green is linked to the remobilization of N from leaves for grain filling leading to leaf senescence [19]. Slower or inhibited remobilization leads to stay-green and while it is linked with drought stress tolerance in some crops (sorghum for example) this is not always the case [19]. Interestingly, the SNPs associated to stay-green are located in a gene encoding a potential uroporphyrin-III C-methyltransferase (UPM), an enzyme of the sirohaem and cobalamin biosynthetic pathway [28]. Sirohaem is an essential component of the plant sulphite and nitrite reductase enzymes [28,29]. It has been proposed that enhanced N uptake might contribute to stay-green and drought tolerance in sorghum [30]. Accordingly, increased nitrite reductase activity in tobacco leads to a stay-green phenotype [31]. We can hypothesize that polymorphisms in the UPM gene might impact sirohaem biosynthesis and as a consequence nitrite reductase activity leading to increased nitrogen assimilation and potentially to increased stay green. Further work will be required to test this hypothesis. No significant association was found for stay green in DS conditions but it might be due to the failure of some lines to set seeds leading to reduced sink demand from grain preventing leaf senescence. Interestingly, the harvest index of the germplasm we tested was also relatively low (around 0.15), suggesting that the sink strength of the panicles was limited in comparison to standard elite genotypes. This suggests that the stay-green we measured is not an indication of a drought stress tolerance but rather an indication of the balance of the remobilization process between the reproductive and vegetative organs.
Pearl millet is a tough crop that can survive and yield in harsh low water and high-temperature conditions [2,3,8]. As such it is a great model to identify mechanisms involved in drought and heat stress tolerance in cereals. Here we identified four loci associated with increased biomass production under early drought stress. These associations were dependent on the inclusion of one line with a strong BSA phenotype under drought stress so they will need to be confirmed. Genes encoding proteins involved in signal transduction pathways related to drought were found to co-localize with two of these loci while a third associated SNP was found in a gene encoding an enzyme involved in threonine catabolism. The last loci associated with increased biomass under drought stress contained a gene encoding an enzyme (3-ketoacyl-CoA synthase or KCS) that catalyzes the elongation of C24 fatty acids during both wax and suberin biosynthesis [20]. Cuticular waxes are the main transpiration barriers in leaves and a correlation between wax content and drought tolerance has been reported in many crops [32]. Suberin is a polymer made of aliphatic and aromatic compounds and confers hydrophobic characteristic to the walls of certain root cells (in particular in the endodermis), providing barrier properties against water diffusion [33]. It has been associated with drought tolerance in some species. For instance, the Arabidopsis enhanced suberin 1 (esb1) mutant, characterized by twofold increased root suberin content, has increased water-use efficiency and drought tolerance [34,35]. Suberin deposits have been observed in the endodermis and exodermis in Response to early drought stress in pearl millet pearl millet roots [36]. Interestingly, the gene families related to wax, cutin and suberin biosynthesis and transport are expanded in pearl millet compared to other cereal crops [2]. It has been proposed that this might have contributed to heat and drought tolerance in pearl millet [2,3]. The KCS gene in the locus associated to biomass production under early drought stress might be related to these tolerance mechanisms. The material we identified in this study will be instrumental to test this hypothesis and it will be interesting to analyze in future studies the relation between drought tolerance and/or suberin composition and content in lines contrasted for drought tolerance.
In conclusion, our results reveal new potential mechanisms regulating response to drought stress and the stay green trait in pearl millet and might contribute to breed pearl millet lines with improved yield and drought stress tolerance.