Transcriptome Analysis of Pennisetum glaucum (L.) R. Br. Provides Insight Into Heat Stress Responses

Pennisetum glaucum (L.) R. Br., being widely grown in dry and hot weather, frequently encounters heat stress at various stages of growth. The crop, due to its inherent capacity, efficiently overcomes such stress during vegetative stages. However, the same is not always the case with the terminal (flowering through grain filling) stages of growth, where recovery from stress is more challenging. However, certain pearl millet genotypes such as 841-B are known to overcome heat stress even at the terminal growth stages. Therefore, we performed RNA sequencing of two contrasting genotypes of pearl millet (841-B and PPMI-69) subjected to heat stress (42°C for 6 h) at flowering stages. Over 274 million high quality reads with an average length of 150 nt were generated, which were assembled into 47,310 unigenes having an average length of 1,254 nucleotides, N50 length of 1853 nucleotides, and GC content of 53.11%. Blastx resulted in the annotation of 35,628 unigenes, and functional classification showed 15,950 unigenes designated to 51 Gene Ontology terms. A total of 13,786 unigenes were allocated to 23 Clusters of Orthologous Groups, and 4,255 unigenes were distributed to 132 functional Kyoto Encyclopedia of Genes and Genomes database pathways. A total of 12,976 simple sequence repeats and 305,759 SNPs were identified in the transcriptome data. Out of 2,301 differentially expressed genes, 10 potential candidate genes were selected based on log2 fold change and adjusted p value parameters for their differential gene expression by qRT-PCR. We were able to identify differentially expressed genes unique to either of the two genotypes, and also, some DEGs common to both the genotypes were enriched. The differential expression patterns suggested that 841-B 6 h has better ability to maintain homeostasis during heat stress as compared to PPMI-69 6 h. The sequencing data generated in this study, like the SSRs and SNPs, shall serve as an important resource for the development of genetic markers, and the differentially expressed heat responsive genes shall be used for the development of transgenic crops.

Introduction

Temperature is a key physical parameter which affects the growth and development of a plant. According to the Intergovernmental Panel on Climate Change (IPCC), there has been an average increase of 4°C in global atmospheric temperature since the late 20th century (Porter et al., 2014). A plant undergoes a number of morphological, physiological, biochemical, and molecular changes during heat stress to ensure its survival (Wahid et al., 2007). These changes include reduction in chlorophyll content, changes in membrane fluidity and protein stability, production of reactive oxygen species (ROS), secondary signaling, and transcriptional changes. In some crops, the occurrence of heat stress during the flowering period leads to poor grain setting (Iqbal et al., 2009; Moriondo et al., 2011). Increased water stress due to heat throughout the growing cycle can reduce the crop yield (Lobell et al., 2013). Despite the impact of high temperature on plant growth and crop yield being quite profound, the underlying heat tolerance mechanisms are not clearly understood in many crops.

Pearl millet [Pennisetum glaucum (L.) R. Br.], also known as Cenchrus americanus (L.) Morrone, is widely grown in the African and Indian subcontinents, since prehistoric times. The crop’s main center of diversity is known to be the Sahel zone of West Africa. Pearl millet is a C4 species having diploid number 2n = 14 and genome size around 1.79 Gb (Varshney et al., 2017) and is mostly grown under drought-prone semi-arid and arid tropics in the regions with 200–800 mm of annual rainfall. Pearl millet is grown in an area of about 31 million hectares worldwide in more than 30 countries with more than 90 million poor people depending on it for food and income (http://exploreit.icrisat.org/profile/Pearl%20Millet/178 accessed on 6 February 2022). Optimum temperature required for the growth of pearl millet is about 30–35°C (Maibam and Padaria, 2015). It is well known for its tolerance to extreme environmental conditions. Limited genomic resources are available as compared to other crop species.

Next-generation sequencing (NGS) based technology for the analysis of transcriptome is more powerful and accurate compared to the Sanger based EST sequencing, suppression subtractive hybridization, and hybridization based microarrays (Marioni et al., 2008; Agarwal et al., 2010). Moreover, over the last decade, several NGS platforms including the illumina are becoming more affordable and efficient for transcriptome sequencing. Additionally, transcriptome sequencing has emerged as an alternative core technology for the discovery and understanding of genes associated with desired traits, where full genome sequencing is not economically feasible especially in case of nonmodel plants. RNA-sequencing (RNA-Seq) has become a benchmark tool for whole transcriptome gene expression quantification and identification of differentially expressed genes (DEGs). It provides scope for the identification of probable candidate genes involved in abiotic and biotic stress tolerance and further for the development of molecular markers (Wang et al., 2009; Garber et al., 2011; Strickler et al., 2012). Third generation sequencing technology such as PacBio provides long/full length transcripts but still having high single error base rate (Sun et al., 2020).

Enormous amounts of ESTs generated from various transcriptomic studies and other genomic sequences are available in public databases for many model plant species (James et al., 2015). However, limited research emphasis has been given to the nonmodel crops including pearl millet, as evidenced by the presence of only 75,499 ESTs (https://www.ncbi.nlm.nih.gov/nuccore/?term=pearl+millet+ESTs, Accessed on 12th April 2022) of this crop in GenBank. It is now possible to assemble transcripts without the reference genome via de novo assembly using trinity (Haas et al., 2014) and/or one of the several other available software tools. Recently, genome wide expression profiling in various nonmodel plant species growing under abiotic or biotic stresses, for various biosynthetic pathways or developmental stages, has been carried out using RNA-sequencing. These plants include Raphanus sativus L. (Nie et al., 2016), Scrophularia ningpoensis Hemsl. (Pan et al., 2017), Camellia sinensis (L.) Kuntze (Guo et al., 2017), Lilium genotypes (Hu et al., 2017), Brassica rapa L. (Greenham et al., 2017), Sesamum indicum L. (Dossa et al., 2017), Asparagus officinalis L. (Li et al., 2017), and Agrostis stolonifera L. (Ma et al., 2017). Transcriptome studies performed for pearl millet have mostly covered drought and/or salinity stresses (Dudhate et al., 2018; Jaiswal et al., 2018; Shinde et al., 2018; Shivhare et al., 2020) and only two studies (Sun et al., 2020, 2021) have covered transcriptome of pearl millet under heat stress. Wherein they used only one variety to study the transcriptome, different parts of the plant viz. roots, leaves, and whole plant were used for RNA extraction. Given the limited number of studies available on transcriptome under heat stress in pearl millet, it would be appropriate to study the transcriptome of this crop under heat stress to gain novel insights into the underlying mechanism. For the said purpose, we chose two varieties having contrasting feature visa-a-vis their heat stress tolerance (heat tolerant genotype-841-B and heat sensitive genotype PPMI-69).

Flag leaf acts as an immediate source for panicle development during the reproductive stage in plants (Wang et al., 2011), RNA-sequencing of flag leaf subjected to heat stress during flowering stage was carried out using Illumina sequencing platform. De novo assembly resulted in 47,310 unigenes and further, functional annotation (gene ontology, corresponding metabolic pathways) was carried out. The aim of the present study was to unravel the gene pool responsible for conferring heat tolerance to pearl millet. Being one among a few transcriptome reports of pearl millet in response to heat stress, the data presented here will be a primary source of information for the research on genomics and functional genomics in this orphan crop.

Materials and Methods

Plant Materials, Heat Stress Treatment, RNA Isolation, and cDNA Library Construction

Seeds of two contrasting genotypes, 841-B and PPMI-69, were collected from Division of Genetics, Indian Agricultural Research Institute (ICAR-IARI), New Delhi, India. The genotype 841-B has higher tolerance to heat stress compared to PPMI-69 (Sankar et al., 2013, 2014). The tolerance of the genotypes was evaluated by membrane stability index (MSI) and Malondialdehyde assay (Heath and Packer, 1968). For heat stress treatment, seeds of both the genotypes were surface-sterilized and sown in plastic pots (10 inches) filled with vermiculite and grown under glasshouse (temperature 32 ± 2°C, relative humidity 70–80%, under day length of 12 h), at National Phytotron facility, IARI, New Delhi. 10 seeds per genotype were grown with one seedling per pot. At flowering stage (55 days after sowing), heat stress was applied in a growth chamber at a temperature of 42°C, relative humidity of 70–80%, and normal light conditions for different time intervals (30 min and 6 h). Plants grown at 32–34°C under normal light conditions in the glasshouse served as control. Different plant samples used in the study were given independent identity numbers (Supplementary Table S1). For RNA extraction, one flag leaf per plant was collected from each plant sample (both untreated and treated) respectively in biological triplicates and was immediately frozen in liquid nitrogen before storing at −80°C. RNA isolation was carried out in three biological replicates using TRIzol reagent (ThermoFisher Scientific, United States) and purified using NucleoTrap mRNA mini kit (Macherey-Nagel, Germany). DNA contamination was removed using TURBO DNase (Ambion, United States) according to the manufacturer’s instructions. The RNA quality was assessed using the 2,100 Bioanalyzer (Agilent Technology, United States). One µg of the total RNA from each sample was used to purify poly-A containing mRNA molecules using Oligotex mRNA mini kit (Qiagen, Germany) as described by the manufacturer. Four independent RNA-seq libraries were constructed using TruSeq^® Stranded mRNA Library Prep Kit (Illumina, United States) according to the manufacturer’s instructions (Supplementary Table S1). The RNA libraries thus constructed were sequenced using Illumina Hiseq platform.

Determination of Physio-Biochemical Characteristics of Plants

Malondialdehyde (MDA) content of both the genotypes was estimated as described previously (Heath and Packer, 1968). Briefly, 0.5 g of leaf tissue from each genotype was taken and homogenized in 10% trichloroacetic acid (TCA) and 0.65% thiobarbituric acid (TBA). The homogenate was incubated at 95°C for 30 min and allowed to cool to room temperature (∼25°C) followed by centrifugation at 10,000 × g for 10 min. The supernatant was collected and its absorbance was measured spectrophotometrically (Shimadzu UV–vis Spectrophotometer UV-1800, Japan) at 532 and 600 nm. The MDA equivalent was calculated in nMg⁻¹ fresh weight as MDA = [(A₅₃₂-A₆₀₀)]/155] × 100. For the estimation of membrane stability index (MSI), the leaf samples were washed with double distilled water (DDW) to remove surface contamination and 10 leaf discs were taken in sealed vials containing 10 ml of DDW separately, followed by incubation at 4°C for 24 h. The electrical conductivity (EC1) was recorded by using a digital conductivity meter (Blum and Ebercon, 1981). For the electrical conductivity (EC2), the samples were autoclaved at 121°C for 20 min, and allowed to cool down to room temperature. The membrane stability index was calculated as per the equation: MSI (%) = 1-(EC1/EC2) × 100.

De Novo Assembly of Flag Leaf Transcriptome

The next-generation sequencing run for whole transcriptome was performed using Paired end (PE) 2 bp × 150 bp library on Illumina HiSeq 2,500. Using FastQC tools (Andrews, 2014), quality check was performed for the raw data shown in Supplementary Table S2. Trimmomatic was used for preprocessing the raw reads to eliminate adapter sequences and poor quality reads. Trinity program was used for de novo assembly with the default parameters (Grabherr et al., 2011; Haas et al., 2014). The Cluster Database at High Identity with Tolerance (CD-HIT) program (Li and Godzik, 2006) was run to remove the similar short sequences based on 90% alignment coverage to longer sequence and produces a set of ‘nonredundant’ (nr) representative sequences and eliminating short redundant sequences. The sequences were clustered using TGICL tools (Pertea et al., 2003) with the default parameters to produce longer, more complete consensus sequences. Gene construction was carried out using EvidentialGene tools with the default parameters to retain the biologically significant transcripts.

Annotation of Transcriptome and Identification of Simple Sequence Repeats and Single-Nucleotide Polymorphisms

The transcriptome structural annotation was performed using TransDecoder (https://github.com/TransDecoder). The functional annotation was performed using BLAST + tools, with BLASTx using a translated nucleotide query (unigenes). Gene Ontology mapping was performed using Blast2GO (Conesa and Gotz, 2008), to specify all the annotated unigenes to various categories such as biological processes, molecular functions, and cellular components. Pathway mapping of unigenes was performed using KEGG database (Ogata et al., 1999). The unigene sequences were aligned to the Clusters of Orthologous Groups (COGs database) (Tatusov et al., 2000) to predict and classify proteins. PlantTFcat online tool (http://plantgrn.noble.org/PlantTFcat/) was used to identify the transcription factor in the generated data. SSRs were identified using MIcroSAtellite (MISA) identification tool. The Microsatellite search module (MISA) is available online for the public (http://pgrc.ipk-ga-tersleben.de/misa/). The SNPs were identified by using GATK best practice pipeline Version 4.1.2.0 (https://software.broadinstitute.org/gatk/best-practices/), the cleaned reads were mapped against the Transfuse. fasta file using BWA aligner (http://bio-bwa.sourceforge.net/). The alignment was performed in default mode. Picard tool was used to co-ordinate, sort, and remove duplicates from the aligned bam files. The GATK tool was used for processing the alignments and variant calling. SplitNCigarReads and HaplotypeCaller from GATK tools were used for reassigning mapping qualities and variant calling, respectively.

Expression Analysis

Fragments per kilobase of transcript per million mapped reads (FPKM) unit was used to calculate the expression level of unigenes. Read count for each unigenes were calculated and then converted to FPKM using the formula: (Read count X 10⁹)/(Sum of read count × Length). Differential gene expression was determined using DEGseq (Wang et al., 2010), in R package. The significant DEGs were filtered based on the adjusted p-value <0.005 and log ratio >1 and −1 between the samples. Heatmaps of the significant genes were generated with the heatmap package (Perry, 2018), in R package. Using R package pvclust (Suzuki and Shimodaira 2006), hierarchical clustering was performed with 1,000 bootstrap replications.

Quantitative RT-PCR

Total RNA was extracted (as described in Section 2.1) from the treated and control flag leaves to study the differential expression patterns under heat stress (42°C for 30 min and 6 h) of the few selected genes. A total of 10 genes (PgDnaJ, PgGST, PgNAC67, PgTIL, PgEXP, PgHd1, PgLTP, PgUCP1, PgUCP2, and PgUCP3) involved in heat stress response were selected for differential expression analysis, from the generated pearl millet transcriptome data based on the log2 fold change ≥2 (for upregulated transcripts) and adjusted p-values. cDNA was synthesized for studying the expression of selected DEGs with the samples of flag leaves originally used for RNA-seq, using SuperScript^® III First-Strand Synthesis System (Invitrogen, United States). qRT-PCR was performed using LightCycler^® 480 System (Roche, Switzerland) and KAPA SYBR^® FAST qPCR Kits (Kapa Biosystems, United States) was used as reaction components. Gene specific primers were designed using PrimerQuest Tool (Integrated DNA Technologies (IDT), United States) (Supplementary Table S3). PCR program was set as: 94°C for 5 min once followed by 40 cycles at 94°C for 10 s, 60°C for 10 s, and 72°C for 10 s. PgActin was used as an endogenous control to normalize all the data (Sankar et al., 2021). 2^−ΔΔCt method was used to calculate the relative fold change (Livak and Schmittgen 2001) among the treatments. The significance levels were calculated using two-tailed unpaired t-test.

Results and Discussion

Determination of Physio-Biochemical Properties of P. glaucum Genotypes

The detection of higher content of Malondialdehyde and lower membrane stability index in genotype PPMI-69 compared to that of 841-B indicated that PPMI-69 genotype is susceptible to heat stress compared to 841-B (Supplementary Figures S1,S2). The heat susceptibility of PPMI-69 genotype compared to that of 841-B has been previously reported (Sankar et al., 2013, 2014).

Illumina Sequencing and Raw Data Pre-Processing

Pear millet, despite being an important source of food and fodder in both arid and semi-arid regions of the world, has not been widely explored as a reservoir of heat tolerant genes. There are some genomic and very few transcriptomic studies exploring the genetic potential of pearl millet as a source of heat stress responsive genes. In the year 2016, Berthouly-Salazar et al. (2016) performed RNA sequencing of different pearl millet populations to explore their climatic adaptability. Varshney et al. (2017) reported whole genome sequencing of pearl millet using shortgun and BAC cloning approaches. Sun et al. (2020) performed Pacbio full-length RNA sequencing of pearl millet under heat and drought stress. Sun et al. (2021) reported root transcriptome of pearl millet under heat stress. Huang et al. (2021) studied the transcriptional changes in pearl millet leaves under heat stress. In this study, the whole transcriptome sequencing was performed using Paired end (PE) 2 bp × 150 bp library on Illumina HiSeq 2,500. The sequencing run produced a total raw data of 288, 876, 956 reads, details are given in Supplementary Table S2. After the removal of low quality sequences, ambiguous bases and adapter sequences by Trimmomatic tool (Bolger et al., 2014), a total of 274, 721, 009 high quality clean reads, containing 39, 782, 593,275 nucleotides (nts) were generated having an average length of 150 nt and GC content of 57.17%. The sequencing data has been deposited to NCBI (National Centre for Biotechnology Information), Sequence Read Archive (SRA) database under the accession number SRP151237.

De Novo Assembly of Pearl Millet Flag Leaf Transcriptome

Using the Trinity program based on the de Bruijn graph algorithm, we performed de novo transcriptome assemblies using their default K-mer sizes. The analysis generated 147,934 contigs having the mean length of 1,059 nt and N50 length of 1,526 nt (Supplementary Table S4), which is significantly higher than the average length (725 bp, N50 1,014 bp) reported by Varshney et al. (2017) and lesser than average length (3,102 bp, N50 3,302 bp) reported by Sun et al. (2020). In order to reduce the assembled contig numbers, CD-HIT software was used for grouping and estimating similarity and dissimilarity of nucleotide sequences, which resulted in the number of contigs being reduced from 147,934 to 129,893 due to the removal of redundant sequences. TGICL tools were further used to retrieve longer and complete contigs, as a result, 129,893 contigs were processed into 109,001 contigs, with N50 length of 1,649 nt. In order to retain the biologically significant contigs, EvidentialGene tools was used and these contigs were assembled in a nonredundant manner. 47,310 high quality unigenes were generated, with a total length of 59, 323, 119 nt, a mean length 1,254 nt, N50 length of 1,853 nt, and GC content 53.11% (Supplementary Table S4). The use of assembly tools (CD-HIT, TGICL, and EvidentialGene) led to the improvement of N50 values, compared to raw assembly. To find the read usage in the assembly, we aligned all the 274.721 million reads to 47,310 unigenes using Bowtie two software tool (Langmead and Salzberg, 2012), 72.62% of reads were aligned to the assembled transcripts (Supplementary Table S5).

Structural and Functional Annotation and Classification of P. glaucum Transcriptome

The transcriptome structural annotation analysis was performed using TransDecoder tool (https://github.com/TransDecoder/TransDecoder/wiki-Date of access 17th January 2020). Out of 47,310 unigenes analyzed, 29,919 (63.24%) were found to be coding sequences, in which 11,893 unigenes (25.13%) were detected with ORFs (Open Reading Frame) (Supplementary Table S6). Functional annotation of all the assembled unigenes were compared to the NCBI nonredundant (nr) protein database (Blastx program) (Altschul et al., 1990) with a cut-off value of 1.0 E-06. A total of 35,628 unigenes (75.31% of all unigenes) were annotated against the nr protein database while the remaining 11,682 (24.69%) were not annotated (Table 1). Based on the E-value distribution of Blastx results, 81.34% of the unigenes showed E-value < 1.0E-45 while 18.66% of unigenes had E-value in the range of 1.0E-0.6 to 1.0E-45 (Figure 1A). 83.35% of the aligned unigenes showed more than 80% of similarity distribution (Figure 1B). Blast top hits analysis showed that 57.07% of the annotated sequences correspond to Setaria italica (L.) P. Beauvois, followed by Sorghum bicolor (L.) Moench (6.61%), Cronobacter sakazakii (4.72%), Zea mays L (4.33%), and Oryza sativa japonica (2.74%) (Figure 2). A similar trend in the distribution of BLAST top hits was observed in pearl millet subjected to heat and drought stress previously (Sun et al., 2020). The 30 top-hit species based on nr annotation are shown in Figure 2. Based on Gene Ontology (GO), 15,950 unigenes were designated into three GO classes i.e., 23 biological processes, 14 cellular components, and 14 molecular functions as shown in Figure 3A. Transcriptional sequences for cellular process, biological regulation, establishment of localization, localization, pigmentation, response to stimulus, and metabolic process, among others were significantly enriched under the biological process. Within the cellular component, unigene sequences for the cell, cell part, organelle, organelle part, and macromolecular complex were identified as highly enriched categories. The major proportion of unigenes was assigned to binding, catalytic activity, and transporter categories under molecular function. Moreover, Clusters of Orthologous Groups (COGs) analysis showed that 13,786 unigenes (29.14% of all unigenes) were allocated to 23 COGs categories (Table 1; Figure 3B). Of these, maximum unigenes fall under the category of unknown functions (4,353), followed by a large number of unigenes falling under the categories of signal transduction mechanisms (1,776), transcription (1,105), carbohydrate metabolism and transport (1,095), and post-translational modification, protein turnover, chaperone functions (961). Minimum unigenes were observed to fall under the categories of extracellular structures (9) and cell motility (3). The assigned function of unigenes showed GO (44.77%) and COGs (29.14%) classifications, representing a broad range of cellular transcripts in pearl millet. We used PlantTFcat tool (http://plantgrn.noble.org/PlantTFcat/) and identified 3,841 unigenes associated with the plant transcription factors (Figure 4).

TABLE 1

TABLE 1. Summary of annotations against publicly available databases.

FIGURE 1

FIGURE 1. Characteristics of sequence homology of unigenes against the NR database. (A)- E-value distribution of BLAST hits for each unigene. (B) Similarity distribution of BLAST hits of each unigene.

FIGURE 2

FIGURE 2. Distribution of the top BLAST hits in different species.

FIGURE 3

FIGURE 3. (A)- GO (Gene Ontology) classification of the transcriptome. (B)- COGs (Clusters of Orthologous Groups) classification.

FIGURE 4

FIGURE 4. Summary of annotation of all unigenes.

Identification of Heat-Responsive Genes Involved in Biological Pathways During Flowering

To identify the potential heat responsive genes and understand their role in various biological pathways, KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis was performed with a cut-off E-value of 1.0E-05. In total, 4,255 unigenes (11.94% of total unigenes) were categorized into 132 KEGG pathways (Table 1; Figure 4). The most represented pathways were the ones related to “biosynthesis of antibiotics” (10.86%), “purine metabolism” (6.46%), “starch and sucrose metabolism” (3.45%), “pyrimidine metabolism” (3.2%), “phenylpropanoid biosynthesis” (2.70%), and “glycolysis/gluconeogenesis” (2.63%). A total of 147 transcripts of starch and sucrose metabolism (3.45%) and 112 transcripts of “glycolysis/gluconeogenesis” (2.63%) were identified. Our study will provide better understanding of molecular mechanisms that are prevailing during the flowering stage of pearl millet under heat stress. Interestingly, this analysis shall help in specifying pathways related to synthesis and turnover of compounds, which have favorable effects in grain filling and yield.

Identification of Simple Sequence Repeats and Single-Nucleotide Polymorphisms

Assembled transcriptome of P. glaucum was used for the identification of SSRs. Based on the criteria with a minimum of (5–10) repetitions of mono to hexa-nucleotide motifs, MISA software was used to search the SSR markers in all unigenes (Thiel et al., 2003). A total of 12,976 SSRs from 10,294 unigenes were detected, of which 2,116 unigenes had more than one SSR (Table 2). Among all the identified SSRs, 50.88% fall under tri-nucleotide repeats, followed by mono-nucleotides (31.88%) and di-nucleotides (13.7%), but tetra-nucleotide, penta-nucleotide, and hexa-nucleotide were represented only as a small fraction (Figure 5).

TABLE 2

TABLE 2. Statistics related to SSRs obtained from the transcriptome of Pennisetum glaucum.

FIGURE 5

FIGURE 5. Distribution of SSRs observed in this study. X-axis represents SSRs repeat motif type. Y-axis represents number of SSRs.

SSRs are simple motifs of nucleotides (1–10 nucleotides), which may occur as tandem or interspersed repeats and are abundant within the genome of prokaryotes and eukaryotes (Vieira et al., 2016). Genetic variability for heat stress tolerance in P. glaucum is still unexplored. Therefore, mining the SSR markers from P. glaucum would be utilized by breeders to develop heat stress tolerant crops. Several studies show that SSRs are not distributed randomly along the genome. For example, in case of Arabidopsis thaliana (L.) Heynh. rice (Lawson and Zhang, 2006) and Gossypium raimondii Ulbr. (Zou et al., 2012), it has been reported that the occurrence of GC-rich trinucleotides SSRs were frequent in exon regions, whereas distribution of AT-rich trinucleotides SSRs were found throughout the genome (coding sequences, untranslated regions, introns, and intergenic spaces) (Temnykh et al., 2001). SSRs are codominant, multi-allele genetic markers that are highly reproducible and transferable among related species (Mason, 2015). As a result, it has been the most widely used marker for genotyping and other breeding purposes (Senthilvel et al., 2008). Identification of new SSRs will provide the necessary impetus to the research community interested in genotyping, genetic mapping, and genetic diversity studies in various Pennisetum species (Senthilvel et al., 2008). A total of 305,759 SNPs were observed in the P. glaucum transcriptome data. Analysis of SNPs showed that 64.76% (581,800/898,460) of the nucleotide changes were transitions, while 35.24% (316,660/898,460) were transversions. The observed transition: transversion (Ts/Tv) ratio is 1.84. The ratio of transition to transversion is expected to be more than one due to the mutational processes in plant genome. The ratio is lower than the estimates from maize (2.5) and Arabidopsis (2.4). In our studies, an average of one variant was observed on every 5,116 bases. The identified SNPs may be utilized as a genomic resource for P. glaucum improved by mining alleles of genes and genome assisted breeding for future genome-wide association studies.

Differential Gene Expression Analysis and Validation of Heat Stress-Responsive Genes

FPKM (Fragments per kilobase per million fragments) unit was used to calculate the expression levels of genes in P. glaucum flag leaves transcriptome. Each sample reads were aligned separately to 47,310 unigenes. 75.02% reads got aligned (out of 72, 512, 870) for genotype 841-B control and 68.25% reads got aligned (out of 67, 530, 958) for genotype 841-B 6 h, 72.84% reads got aligned (out of 69, 132, 264) for genotype PPMI-69 control and 74.21% reads got aligned (out of 65, 544, 917) for genotype PPMI-69 6 h. After filtering, based on the adjusted p-value, less than 0.005 and (-1< log2 fold change >1) there were 850 DEGs identified between 841-B control and 841-B 6 h, of which 494 genes were upregulated and 356 genes were downregulated. Among 1,934 DEGs identified between PPMI-69 control and PPMI-69 6 h, 964 genes were upregulated and 970 genes were downregulated (Figure 6B).

FIGURE 6

FIGURE 6. Representation of heat stress responsive genes. (A)- Heatmaps of the significant differentially expressed genes with hierarchical clustering. (B)- Venn diagram of the differential expressed genes under heat stress conditions in 841-B and PPMI-69 genotypes.

Differential gene expression patterns across the treatments and comparison of differential expression patterns between the two different genotypes, were analyzed by cluster analysis with the hierarchical clustering method (Figure 6A) and the DEGs were visualized by volcano plots (Figure 7). The hierarchical cluster analysis of two different genotypes with heat stress treatment shows that the differential gene expression pattern of the transcriptome readily differentiates P. glaucum genotype PPMI-69 6 h from others (841-B control, 841-B 6 h, and PPMI-69 control) (Figure 6A), possibly indicating the variation in differential gene expression occurring between the genotypes in response to heat stress. The maximum number of DEGs (1,934 genes) was observed between PPMI-69 6 h and PPMI-69 control (Table 3). These DEGs between the samples 841-B 6 h and PPMI-69 6 h are grouped into two clusters (Figure 6A), representing that the genes involved in thermotolerance have different level of differential gene expression. In general, resemblance in the pattern of DEGs was observed for the samples 841-B 6 h, 841-B control, and PPMI-69. This suggests that 841-B 6 h has better ability to maintain homeostasis during heat stress as compared to PPMI-69 6 h. A comparative analysis between the differentially expressed genes in 841-B and PPMI-69 genotypes (Figure 6B) was performed to identify the common DEGs in both the genotypes under heat stress and those that are unique to each genotype of the P. glaucum. A total of 2,301 genes were differentially expressed, 20.99% (483 genes) of which were shared common by both the genotypes. 15.95% (367 genes) of the DEGs were unique to 841-B genotypes and 63.06% (1,451 genes) of the DEGs were unique to PPMI-69. However, the differential expression of genes between the two genotypes cannot be attributed exclusively to the treatments alone, as the two genotypes have inherently different levels of heat stress tolerance and therefore could have different gene expression profile. The sequencing data generated in this study is being utilized for mining and validation of heat stress responsive genes in different varieties of pearl millet including 841-B and PPMI-69 varieties (Sankar et al., 2021).

FIGURE 7

FIGURE 7. Volcano plots displaying the differentially expressed genes. (A)- 841-B genotype. (B)- PPMI-69 genotype. Red color represents upregulated genes, green color indicates downregulated genes, while nonsignificant genes are shown as black dots.

TABLE 3

TABLE 3. Statistics related to the number of significantly enriched differentially expressed genes (DEGs).

In order to investigate the temporal expression patterns by qRT-PCR analysis, 10 target genes were selected based on the fold changes and adjusted p-values. The selected genes displayed varying patterns in response to different durations of heat stress (Figure 8). The validation of selected 10 genes (seven known and three uncharacterized) by qRT-PCR deciphered their significant role in heat stress management in P. glaucum. These genes might play essential role(s) in the amelioration of heat stress in P. glaucum.

FIGURE 8

FIGURE 8. Real time PCR validation of 10 target genes. (A)- In genotype 841-B. (B)- In genotype PPMI-69. Two tailed unpaired t-test was used to calculate p value, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001.

Heat stress induces change in protein conformation and increase improper folding of the native proteins. As a result, accumulation of many heat shock proteins is triggered to counterbalance the negative effect of heat stress in plants (Schlesinger, 1990). Among the many heat shock proteins, DnaJ or Hsp40 is known to play significant role in plant development, signal transduction, and abiotic stresses response, either by itself or in association with Hsp70 (Yang et al., 2010; Shen et al., 2011; Sankar et al., 2021). PgDnaJ expression show significant upregulation, about 15 folds in 841-B in response to 30 min heat stress, and 16 folds in PPMI-69 (Figure 8). However, this result indicates the important role of DnaJ in maintaining cellular protein homeostasis during heat stress (42°C). Plasma membrane acts as the primary sensor of the cell when the cell is subjected to heat stress (Murata and Los, 1997). Consequently, membrane properties undergo a number of changes in its composition, ion concentration, and ion channels in response to heat stress (Saidi et al., 2010; Prasad, 1996). Peroxidation of lipid membrane is one of the most important change that occurs in the cell in response to various stresses (Thompson et al., 1998). MDA content is a direct indicator of lipid peroxidation. These byproducts of lipid peroxidation, such as electrophiles or xenobiotics, are detoxified by reactive oxygen species (ROS)-scavenging enzymes such glutathione S-transferases (GST), superoxide dismutases, catalases, and ascorbate peroxidases (APX). The expression of PgGST in our study was significantly upregulated (7 folds) in response to heat stress in 841-B but significantly downregulated in PPMI-69 (Figure 8), indicating its role in the heat stress tolerance pathway. This result correlates with MDA content (Supplementary Figure S2), as it was observed to be higher in PPMI-69 compared with 841-B. NAC transcription factors (TFs) (NAM, No apical meristem; ATAF, Arabidopsis transcription activation factor and CUC, Cup-shaped cotyledon) play important role in plant growth and development and in regulating response to abiotic or biotic stresses (Souer et al., 1996; Nuruzzaman et al., 2013). Among these NAC genes, NAC67 has a role in imparting tolerance to multiple abiotic stress such as drought, salt, and cold stresses. NAC67 has been reported to be involved in conferring tolerance to abiotic stresses in rice (Rahman et al., 2016). However, so far, there has been no report indicating its role in heat stress tolerance. In this study, PgNAC67 expression was found to be significantly upregulated (27 folds) in 841-B in response to 30 min and six folds in PPMI-69 in response to 6 h of heat stress (Figure 8). It shows that PgNAC67 plays a role in heat stress response in P. glaucum. Temperature-induced lipocalins (TILs), a plasma membrane protein, have an important role in basal and acquired thermotolerance in plant. TILs alleviate the heat induced lipid peroxidation in membrane. PgTIL expression was observed to be significantly upregulated (nine folds) in response to 6 h of heat stress in 841-B but significantly downregulated in PPMI-69 (Figure 8). This might be the reason why 841-B is able to maintain a low MDA content (Supplementary Figure S2) as compared to PPMI-69. PgEXP expression shows five folds significant upregulation in 841-B in response to 6 h of heat stress (Figure 8). Association of expression genes and heat stress tolerance in some plants has been reported (Xu et al., 2007; Kuluev et al., 2016). Overexpression of the EXP1 gene exhibits low electrolyte leakage, decrease in membrane lipid peroxidation but higher chlorophyll content, net photosynthetic rate, relative water content, and activity of antioxidant enzyme in transgenic plants (Xu et al., 2014). Heading 1 (Hd1) and early heading 1 (Ehd1) are mainly known for the regulation of flower development and flowering, leading to either induction or suppression corresponding to the particular photoperiod (Endo-Higashi and Izawa, 2011). Environmental factors such as day length and abiotic and biotic stresses regulate the expression of these genes. Previous studies show that inhibition of early heading 1 (Ehd1) in response to drought stress delays flowering in rice (Cho et al., 2017). In our studies, PgHd1 expression shows significant downregulation in both the genotypes in response to 30 min and 6 h heat stress (Figure 8). PgLTP expression shows significant (19 folds) upregulation in response to 6 h of heat stress in PPMI-69 but significantly downregulated in 841-B (Figure 8), indicating the involvement of high activity with regard to transfer of lipid molecules in the cell. This shows active regulation of membrane fluidity in PPMI-69 in response to heat stress. Lipid Transfer Protein (LTP) are reported to be involved with growth and development, response to abiotic and biotic stresses but their functions remain unclear (Debono et al., 2009; Guo et al., 2013; Yu et al., 2013; Pan et al., 2016). Moreover, LTPs has the ability to facilitate the transfer of phospholipids between membranes in vitro (Kader 1996). Among the differentially expressed contigs, uncharacterized ORFs share maximum proportion in our study. Uncharacterized genes with a predicted protein domain associated with zinc fingers (ZnF), ribonuclease (RN), and chaperone were validated for their expression during heat stress. Some of the uncharacterized genes were expressed uniquely either in 841-B or PPMI-69. For example, uncharacterized gene PgUCP1, with the predicted zinc fingers (ZnF) domain is significantly upregulated (34 folds) in 841-B and (10 folds) in PPMI-69 in response to heat stress (Figure 8). ZnF is known for the involvement in multiple stress response, but their exact molecular mechanism and their interaction is yet to be deciphered (Vij and Tyagi 2008; Bogamuwa and Jang, 2016; Maibam et al., 2019). On the contrary, uncharacterized gene PgUCP2, with the predicted ribonuclease (RN) domain is significantly downregulated in 841-B but significantly upregulated (29 folds) in PPMI-69 in response to heat stress (Figure 8). It has been reported that loss of ribonuclease function in Arabidopsis enhances its heat stress tolerance (Nguyen et al., 2015). However, uncharacterized gene PgUCP3, with the predicted chaperone domain show considerable upregulation (4 folds) in 841-B and is downregulated in PPMI-69 in response to heat stress. These uncharacterized genes could likely represent the important genes involved in imparting variation in thermotolerance among different genotypes. Furthermore, detailed investigation of these uncharacterized genes is required for understanding the role in response to stress.

Conclusion

This study investigated the transcriptome profile of pearl millet flag leaves in response to heat stress. In this study, high quality 47,310 unigenes were generated and annotated. This data will provide the foundation for research on gene expression, genomics, and functional genomics in pearl millet improvement program. Furthermore, the SSRs obtained in this study shall facilitate the research on genotyping, and diversity studies of this important crop. The candidate genes, whose expression patterns were validated by qRT-PCR, included seven genes known to have role in heat stress in Pennisetum and/or other crops and three uncharacterized genes whose role is yet to be established in heat stress. All these genes along with the pool of differentially expressed DEGs between two genotypes comprise important resource that can be explored for their effective utilization in the development of transgenic crops tolerant to heat stress.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

AM: data curation and writing—original draft preparation, SL: formal analysis and writing—review and editing, KG: writing—review and editing, SN, AS, MS, SS, and MD: formal analysis, and JP: conceptualization, resources, project administration, writing—review and editing, and funding acquisition.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Award of fellowship to the first author and third author by the Department of Biotechnology, Government of India, and by PG School-IARI, respectively, is duly acknowledged. The authors are grateful to the Director, ICAR-National Institute for Plant Biotechnology (ICAR-NIPB), New Delhi, for providing necessary facilities. The Director, ICAR-Indian Agricultural Research Institute (ICAR-IARI), is thanked for facilitating the use of National Phytotron Facility for growing seedlings. Thanks are also due to Dr. Tara Satyavathi Chellapilla, Division of Genetics, IARI, New Delhi, for providing the seed material.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.884106/full#supplementary-material

Abbreviations

nt, nucleotides; h, hours; BLAST, Basic Local Alignment Search Tool; SSRs, simple sequence repeats; mm, milli metre; ESTs, expressed sequence tags.

References

Agarwal, A., Koppstein, D., Rozowsky, J., Sboner, A., Habegger, L., Hillier, L. W., et al. (2010). Comparison and Calibration of Transcriptome Data from RNA-Seq and Tiling Arrays. BMC Genomics 11, 383. doi:10.1186/1471-2164-11-383

PubMed Abstract | CrossRef Full Text | Google Scholar

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic Local Alignment Search Tool. J. Mol. Biol. 215, 403–410. doi:10.1016/S0022-2836(05)80360-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrews, S. (2014). FastQC: A Quality Control Tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

Google Scholar

Berthouly-Salazar, C., Thuillet, A.-C., Rhoné, B., Mariac, C., Ousseini, I. S., Couderc, M., et al. (2016). Genome Scan Reveals Selection Acting on Genes Linked to Stress Response in Wild Pearl Millet. Mol. Ecol. 25 (21), 5500–5512. doi:10.1111/mec.13859

PubMed Abstract | CrossRef Full Text | Google Scholar

Blum, A., and Ebercon, A. (1981). Cell Membrane Stability as a Measure of Drought and Heat Tolerance in Wheat1. Crop Sci. 21, 43–47. doi:10.2135/cropsci1981.0011183x002100010013x

CrossRef Full Text | Google Scholar

Bogamuwa, S., and Jang, J.-C. (2016). Plant Tandem CCCH Zinc finger Proteins Interact with ABA, Drought, and Stress Response Regulators in Processing-Bodies and Stress Granules. PLoS One 11, e0151574. doi:10.1371/journal.pone.0151574

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics 30, 2114–2120. doi:10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, L.-H., Yoon, J., and An, G. (2017). The Control of Flowering Time by Environmental Factors. Plant J. 90, 708–719. doi:10.1111/tpj.13461

PubMed Abstract | CrossRef Full Text | Google Scholar

Conesa, A., and Götz, S. (2008). Blast2GO: A Comprehensive Suite for Functional Analysis in Plant Genomics. Int. J. Plant Genomics 2008, 1–12. doi:10.1155/2008/619832

CrossRef Full Text | Google Scholar

Debono, A., Yeats, T. H., Rose, J. K. C., Bird, D., Jetter, R., Kunst, L., et al. (2009). Arabidopsis LTPG Is a Glycosylphosphatidylinositol-Anchored Lipid Transfer Protein Required for Export of Lipids to the Plant Surface. Plant Cell 21, 1230–1238. doi:10.1105/tpc.108.064451

PubMed Abstract | CrossRef Full Text | Google Scholar

Dossa, K., Li, D., Wang, L., Zheng, X., Liu, A., Yu, J., et al. (2017). Transcriptomic, Biochemical and Physio-Anatomical Investigations Shed More Light on Responses to Drought Stress in Two Contrasting Sesame Genotypes. Sci. Rep. 7, 8755. doi:10.1038/s41598-017-09397-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Dudhate, A., Shinde, H., Tsugama, D., Liu, S., and Takano, T. (2018). Transcriptomic Analysis Reveals the Differentially Expressed Genes and Pathways Involved in Drought Tolerance in Pearl Millet [Pennisetum Glaucum (l.) R. Br]. PLoS One 13, e0195908. doi:10.1371/journal.pone.0195908

PubMed Abstract | CrossRef Full Text | Google Scholar

Endo-Higashi, N., and Izawa, T. (2011). Flowering Time Genes Heading Date 1 and Early Heading Date 1 Together Control Panicle Development in rice. Plant Cel Physiol 52, 1083–1094. doi:10.1093/pcp/pcr059

PubMed Abstract | CrossRef Full Text | Google Scholar

Garber, M., Grabherr, M. G., Guttman, M., and Trapnell, C. (2011). Computational Methods for Transcriptome Annotation and Quantification Using RNA-Seq. Nat. Methods 8, 469–477. doi:10.1038/nmeth.1613

PubMed Abstract | CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length Transcriptome Assembly from RNA-Seq Data without a Reference Genome. Nat. Biotechnol. 29, 644–652. doi:10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenham, K., Guadagno, C. R., Gehan, M. A., Mockler, T. C., Weinig, C., Ewers, B. E., et al. (2017). Temporal Network Analysis Identifies Early Physiological and Transcriptomic Indicators of Mild Drought in Brassica Rapa. Elife 6. eLife. doi:10.7554/eLife.29655

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, C., Ge, X., and Ma, H. (2013). The rice OsDIL Gene Plays a Role in Drought Tolerance at Vegetative and Reproductive Stages. Plant Mol. Biol. 82, 239–253. doi:10.1007/s11103-013-0057-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, F., Guo, Y., Wang, P., Wang, Y., and Ni, D. (2017). Transcriptional Profiling of Catechins Biosynthesis Genes during tea Plant Leaf Development. Planta 246, 1139–1152. doi:10.1007/s00425-017-2760-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2014). De Novo transcript Sequence Reconstruction from RNA-Seq Using the Trinity Platform for Reference Generation and Analysis. Nat. Protoc. 8, 1494–1512. doi:10.1038/nprot.2013.084

PubMed Abstract | CrossRef Full Text | Google Scholar

Heath, R. L., and Packer, L. (1968). Photoperoxidation in Isolated Chloroplasts. Arch. Biochem. Biophys. 125, 189–198. doi:10.1016/0003-9861(68)90654-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Z., Tang, B., Wu, Q., Zheng, J., Leng, P., and Zhang, K. (2017). Transcriptome Sequencing Analysis Reveals a Difference in Monoterpene Biosynthesis between Scented Lilium 'Siberia' and Unscented Lilium 'Novano'. Front. Plant Sci. 8, 1351. doi:10.3389/fpls.2017.01351

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D., Sun, M., Zhang, A., Chen, J., Zhang, J., Lin, C., et al. (2021). Transcriptional Changes in Pearl Millet Leaves under Heat Stress. Genes 12, 1716. doi:10.3390/genes12111716

PubMed Abstract | CrossRef Full Text | Google Scholar

Iqbal, M. M., Goheer, M. A., and Khan, M. A. (2009). Climate-change Aspersions on Food Security of Pakistan. Sci. Vis. 15, 15–24.

Google Scholar

Jaiswal, S., Antala, T. J., Mandavia, M. K., Chopra, M., Jasrotia, R. S., Tomar, R. S., et al. (2018). Transcriptomic Signature of Drought Response in Pearl Millet (Pennisetum Glaucum (L.) and Development of Web-Genomic Resources. Sci. Rep. 8, 3382. doi:10.1038/s41598-018-21560-1

PubMed Abstract | CrossRef Full Text | Google Scholar

James, D., Tarafdar, A., Biswas, K., Sathyavathi, T. C., Padaria, J. C., and Kumar, P. A. (2015). Development and Characterization of a High Temperature Stress Responsive Subtractive cDNA Library in Pearl Millet Pennisetum Glaucum (L.) R.Br. Indian J. Exp. Biol. 53, 543–550.

PubMed Abstract | Google Scholar

Kader, J.-C. (1996). Lipid-transfer Proteins in Plants. Annu. Rev. Plant Physiol. Plant Mol. Biol. 47, 627–654. doi:10.1146/annurev.arplant.47.1.627

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuluev, B., Avalbaev, A., Mikhaylova, E., Nikonorov, Y., Berezhneva, Z., and Chemeris, A. (2016). Expression Profiles and Hormonal Regulation of Tobacco Expansin Genes and Their Involvement in Abiotic Stress Response. J. Plant Physiol. 206, 1–12. doi:10.1016/j.jplph.2016.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods 9, 357–359. doi:10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson, M., and Zhang, L. (2006). Distinct Patterns of SSR Distribution in the Arabidopsis thaliana and rice Genomes. Genome Biol. 7, R14. doi:10.1186/gb-2006-7-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S.-F., Zhang, G.-J., Zhang, X.-J., Yuan, J.-H., Deng, C.-L., and Gao, W.-J. (2017). Comparative Transcriptome Analysis Reveals Differentially Expressed Genes Associated with Sex Expression in Garden asparagus (Asparagus Officinalis). BMC Plant Biol. 17, 143. doi:10.1186/s12870-017-1091-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., and Godzik, A. (2006). Cd-hit: a Fast Program for Clustering and Comparing Large Sets of Protein or Nucleotide Sequences. Bioinformatics 22, 1658–1659. doi:10.1093/bioinformatics/btl158

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 25, 402–408. doi:10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobell, D. B., Hammer, G. L., McLean, G., Messina, C., Roberts, M. J., and Schlenker, W. (2013). The Critical Role of Extreme Heat for maize Production in the United States. Nat. Clim Change 3, 497–501. doi:10.1038/nclimate1832

CrossRef Full Text | Google Scholar

Ma, Y., Shukla, V., and Merewitz, E. B. (2017). Transcriptome Analysis of Creeping Bentgrass Exposed to Drought Stress and Polyamine Treatment. PLoS One 12, e0175848. doi:10.1371/journal.pone.0175848

PubMed Abstract | CrossRef Full Text | Google Scholar

Maibam, A., Vishwakarma, H., and Padaria, J. C. (2019). In Silico studies Predict Role of PgUCP1 from Pennisetum Glaucum in Heat Stress Tolerance. Indian J. Agric. Sci. 89, 1703–1707.

Google Scholar

Maibam, A., and Padaria, J. C. (2015). Pearl Millet: A Genetic Resource for Abiotic Tolerant Transgenics. Biot. Today 5, 21–24. doi:10.5958/2322-0996.2015.00003.4

CrossRef Full Text | Google Scholar

Marioni, J. C., Mason, C. E., Mane, S. M., Stephens, M., and Gilad, Y. (2008). RNA-seq: An Assessment of Technical Reproducibility and Comparison with Gene Expression Arrays. Genome Res. 18, 1509–1517. doi:10.1101/gr.079558.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Mason, A. S. (2015). SSR Genotyping. Methods Mol. Biol. (Clifton, N.J.) 1245, 77–89. doi:10.1007/978-1-4939-1966-6_6

PubMed Abstract | CrossRef Full Text | Google Scholar

Moriondo, M., Giannakopoulos, C., and Bindi, M. (2011). Climate Change Impact Assessment: The Role of Climate Extremes in Crop Yield Simulation. Climatic Change 104, 679–701. doi:10.1007/s10584-010-9871-0

CrossRef Full Text | Google Scholar

Mukesh Sankar, S., Tara Satyavathi, C., Singh, S. P., Singh, M. P., Bharadwaj, C., and Barthakur, S. (2014). Genetic Diversity Analysis for High Temperature Stress Tolerance in Pearl Millet [Pennisetum Glaucum (L.) R. Br]. Ind. J. Plant Physiol. 19, 324–329. doi:10.1007/s40502-014-0099-2

CrossRef Full Text | Google Scholar

Murata, N., and Los, D. A. (1997). Membrane Fluidity and Temperature Perception. Plant Physiol. 115, 875–879. doi:10.1104/pp.115.3.875

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, A. H., Matsui, A., Tanaka, M., Mizunashi, K., Nakaminami, K., Hayashi, M., et al. (2015). Loss of Arabidopsis 5′-3′ Exoribonuclease AtXRN4 Function Enhances Heat Stress Tolerance of Plants Subjected to Severe Heat Stress. Plant Cel Physiol 56, 1762–1772. doi:10.1093/pcp/pcv096

PubMed Abstract | CrossRef Full Text | Google Scholar

Nie, S., Li, C., Xu, L., Wang, Y., Huang, D., Muleke, E. M., et al. (2016). De Novo transcriptome Analysis in Radish (Raphanus Sativus L.) and Identification of Critical Genes Involved in Bolting and Flowering. BMC Genomics 17, 389. doi:10.1186/s12864-016-2633-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Nuruzzaman, M., Sharoni, A. M., and Kikuchi, S. (2013). Roles of NAC Transcription Factors in the Regulation of Biotic and Abiotic Stress Responses in Plants. Front. Microbiol. 4, 248. doi:10.3389/fmicb.2013.00248

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., and Kanehisa, M. (1999). KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 27, 29–34. doi:10.1093/nar/27.1.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, Y., Chen, D. X., Song, X. H., Zhang, X., and Li, L. Y. (2017). Transcriptome Characterization for Scrophularia Ningpoensis Based on High-Throughput Sequencing Technology and Related Genes for Synthesis of Terpenoid Compounds. Zhongguo Zhong Yao Za Zhi 42, 2460–2466. doi:10.19540/j.cnki.cjcmm.20170614.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, Y., Li, J., Jiao, L., Li, C., Zhu, D., and Yu, J. (2016). A Non-specific Setaria Italica Lipid Transfer Protein Gene Plays a Critical Role under Abiotic Stress. Front. Plant Sci. 7, 1752. doi:10.3389/fpls.2016.01752

PubMed Abstract | CrossRef Full Text | Google Scholar

Perry, M. (2018). Heatmaps : Flexible Heatmaps for Functional Genomics and Sequence Features. R package version 1.4.0.

Google Scholar

Pertea, G., Huang, X., au, F., Antonescu, V., Sultana, R., Karamycheva, S., et al. (2003). TIGR Gene Indices Clustering Tools (TGICL): a Software System for Fast Clustering of Large EST Datasets. Bioinformatics 19, 651–652. doi:10.1093/bioinformatics/btg034

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, J. R., Xie, L., Challinor, A. J., Cochrane, K., Howden, S. M., Iqbal, M. M., et al. (2014). “Food Security and Food Production Systems,” in in Climate Change 2014: Impacts, Adaptation, and Vulnerability. The, Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of Change, Intergovernmental Panel on Climate. Editors M. D. M. Field, C.B, V. R. Barros, D. J. Dokken, K. J. Mach, S. M. T. E. Bilir, M. Chatterjeeet al. (Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press), 485–533.

Google Scholar

Prasad, T. K. (1996). Mechanisms of Chilling-Induced Oxidative Stress Injury and Tolerance in Developing maize Seedlings: Changes in Antioxidant System, Oxidation of Proteins and Lipids, and Protease Activities. Plant J. 10, 1017–1026. doi:10.1046/j.1365-313x.1996.10061017.x

CrossRef Full Text | Google Scholar

Rahman, H., Ramanathan, V., Nallathambi, J., Duraialagaraja, S., and Muthurajan, R. (2016). Over-expression of a NAC 67 Transcription Factor from finger Millet (Eleusine Coracana L.) Confers Tolerance against Salinity and Drought Stress in rice. BMC Biotechnol. 16, 35. doi:10.1186/s12896-016-0261-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Saidi, Y., Peter, M., Finka, A., Cicekli, C., Vigh, L., and Goloubinoff, P. (2010). Membrane Lipid Composition Affects Plant Heat Sensing and Modulates Ca2+-dependent Heat Shock Response. Plant Signaling Behav. 5, 1530–1533. doi:10.4161/psb.5.12.13163

CrossRef Full Text | Google Scholar

Sankar, S. M., Satyavathi, C. T., Singh, S. P., Singh, M. P., et al. (2013). Genetic Variability and Association Studies in Pearl Millet for Grain Yield and High Temperature Stress Tolerance. Indian J. Dryland Agric. Res. Develop. 28, 71–76.

Google Scholar

Sankar, S. M., Satyavathi, T. C., Barthakur, S., Singh, S. P., Bharadwaj, C., and Soumya, S. L. (2021). Differential Modulation of Heat-Inducible Genes across Diverse Genotypes and Molecular Cloning of a sHSP from Pearl Millet [Pennisetum Glaucum (L.) R. Br.]. Front. Plant Sci. 1333.

Google Scholar

Schlesinger, M. J. (1990). Heat Shock Proteins. J. Biol. Chem. 265, 12111–12114. doi:10.1016/s0021-9258(19)38314-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Senthilvel, S., Jayashree, B., Mahalakshmi, V., Kumar, P. S., Nakka, S., Nepolean, T., et al. (2008). Development and Mapping of Simple Sequence Repeat Markers for Pearl Millet from Data Mining of Expressed Sequence Tags. BMC Plant Biol. 8, 119. doi:10.1186/1471-2229-8-119

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, L., Kang, Y. G. G., Liu, L., and Yu, H. (2011). The J-Domain Protein J3 Mediates the Integration of Flowering Signals in Arabidopsis. Plant Cell 23, 499–514. doi:10.1105/tpc.111.083048

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinde, H., Tanaka, K., Dudhate, A., Tsugama, D., Mine, Y., Kamiya, T., et al. (2018). Comparative De Novo Transcriptomic Profiling of the Salinity Stress Responsiveness in Contrasting Pearl Millet Lines. Environ. Exp. Bot. 155, 619–627. doi:10.1016/j.envexpbot.2018.07.008

CrossRef Full Text | Google Scholar

Shivhare, R., Asif, M. H., and Lata, C. (2020). Comparative Transcriptome Analysis Reveals the Genes and Pathways Involved in Terminal Drought Tolerance in Pearl Millet. Plant Mol. Biol. 103, 639–652. doi:10.1007/s11103-020-01015-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Souer, E., van Houwelingen, A., Kloos, D., Mol, J., and Koes, R. (1996). The No Apical Meristem Gene of Petunia Is Required for Pattern Formation in Embryos and Flowers and Is Expressed at Meristem and Primordia Boundaries. Cell 85, 159–170. doi:10.1016/s0092-8674(00)81093-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Strickler, S. R., Bombarely, A., and Mueller, L. A. (2012). Designing a Transcriptome Next-Generation Sequencing Project for a Nonmodel Plant Species1. Am. J. Bot. 99, 257–266. doi:10.3732/ajb.1100292

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, M., Huang, D., Zhang, A., Khan, I., Yan, H., Wang, X., et al. (2020). Transcriptome Analysis of Heat Stress and Drought Stress in Pearl Millet Based on Pacbio Full-Length Transcriptome Sequencing. BMC Plant Biol. 20, 1–15. doi:10.1186/s12870-020-02530-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, M., Lin, C., Zhang, A., Wang, X., Yan, H., Khan, I., et al. (2021). Transcriptome Sequencing Revealed the Molecular Mechanism of Response of Pearl Millet Root to Heat Stress. J. Agro. Crop Sci. 00, 1–6. doi:10.1111/jac.12496

CrossRef Full Text | Google Scholar

Suzuki, R., and Shimodaira, H. (2006). Pvclust: an R Package for Assessing the Uncertainty in Hierarchical Clustering. Bioinformatics 22, 1540–1542. doi:10.1093/bioinformatics/btl117

PubMed Abstract | CrossRef Full Text | Google Scholar

Tatusov, R. L., Galperin, M. Y., Natale, D. A., and Koonin, E. V. (2000). The COG Database: a Tool for Genome-Scale Analysis of Protein Functions and Evolution. Nucleic Acids Res. 28, 33–36. doi:10.1093/nar/28.1.33

PubMed Abstract | CrossRef Full Text | Google Scholar

Temnykh, S., DeClerck, G., Lukashova, A., Lipovich, L., Cartinhour, S., and McCouch, S. (2001). Computational and Experimental Analysis of Microsatellites in rice (Oryza Sativa L.): Frequency, Length Variation, Transposon Associations, and Genetic Marker Potential. Genome Res. 11, 1441–1452. doi:10.1101/gr.184001

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiel, T., Michalek, W., Varshney, R., and Graner, A. (2003). Exploiting EST Databases for the Development and Characterization of Gene-Derived SSR-Markers in Barley (Hordeum Vulgare L.). Theor. Appl. Genet. 106, 411–422. doi:10.1007/s00122-002-1031-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, J. E., Froese, C. D., Madey, E., Smith, M. D., and Hong, Y. (1998). Lipid Metabolism during Plant Senescence. Prog. Lipid Res. 37, 119–141. doi:10.1016/s0163-7827(98)00006-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshney, R. K., Shi, C., Thudi, M., Mariac, C., Wallace, J., Qi, P., et al. (2017). Pearl Millet Genome Sequence Provides a Resource to Improve Agronomic Traits in Arid Environments. Nat. Biotechnol. 35, 969–976. doi:10.1038/nbt.3943

PubMed Abstract | CrossRef Full Text | Google Scholar

Vieira, M. L. C., Santini, L., Diniz, A. L., and Munhoz, C. d. F. (2016). Microsatellite Markers: What They Mean and Why They Are So Useful. Genet. Mol. Biol. 39, 312–328. doi:10.1590/1678-4685-gmb-2016-0027

PubMed Abstract | CrossRef Full Text | Google Scholar

Vij, S., and Tyagi, A. K. (2008). A20/AN1 Zinc-finger Domain-Containing Proteins in Plants and Animals Represent Common Elements in Stress Response. Funct. Integr. Genomics. 8, 301–307. doi:10.1007/s10142-008-0078-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahid, A., Gelani, S., Ashraf, M., and Foolad, M. (2007). Heat Tolerance in Plants: An Overview. Environ. Exp. Bot. 61, 199–223. doi:10.1016/j.envexpbot.2007.05.011

CrossRef Full Text | Google Scholar

Wang, L., Feng, Z., Wang, X., Wang, X., and Zhang, X. (2010). DEGseq: an R Package for Identifying Differentially Expressed Genes from RNA-Seq Data. Bioinformatics 26, 136–138. doi:10.1093/bioinformatics/btp612

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Cai, J., Jiang, D., Liu, F., Dai, T., and Cao, W. (2011). Pre-anthesis High-Temperature Acclimation Alleviates Damage to the Flag Leaf Caused by post-anthesis Heat Stress in Wheat. J. Plant Physiol. 168, 585–593. doi:10.1016/j.jplph.2010.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Gerstein, M., and Snyder, M. (2009). RNA-seq: a Revolutionary Tool for Transcriptomics. Nat. Rev. Genet. 10, 57–63. doi:10.1038/nrg2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Tian, J., Belanger, F. C., and Huang, B. (2007). Identification and Characterization of an Expansin Gene AsEXP1 Associated with Heat Tolerance in C3 Agrostis Grass Species. J. Exp. Bot. 58, 3789–3796. doi:10.1093/jxb/erm229

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Q., Xu, X., Shi, Y., Xu, J., and Huang, B. (2014). Transgenic Tobacco Plants Overexpressing a Grass PpEXP1 Gene Exhibit Enhanced Tolerance to Heat Stress. PLoS ONE 9, e100792. doi:10.1371/journal.pone.0100792

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Qin, Y., Xie, C., Zhao, F., Zhao, J., Liu, D., et al. (2010). The Arabidopsis Chaperone J3 Regulates the Plasma Membrane H+-ATPase through Interaction with the PKS5 Kinase. Plant Cell 22, 1313–1332. doi:10.1105/tpc.109.069609

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, K., Soares, J. M., Mandal, M. K., Wang, C., Chanda, B., Gifford, A. N., et al. (2013). A Feedback Regulatory Loop between G3P and Lipid Transfer Proteins DIR1 and AZI1 Mediates Azelaic-Acid-Induced Systemic Immunity. Cel Rep. 3, 1266–1278. doi:10.1016/j.celrep.2013.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, C., Lu, C., Zhang, Y., and Song, G. (2012). Distribution and Characterization of Simple Sequence Repeats in Gossypium Raimondii Genome. Bioinformation 8, 801–806. doi:10.6026/97320630008801

PubMed Abstract | CrossRef Full Text | Google Scholar