Charlotte Jones1
Jose De Vega2
Margaret Worthington3
Ann Thomas1Dagmara Gasior1
John Harper1
John Doonan1Yuan Fu1
Maurice Bosch1
Fiona Corke1
Jacobo Arango4Juan Andres Cardoso5Juan de la Cruz Jimenez6
Ian Armstead1*
Narcis Fernandez-Fuentes1*- 1Institute of Biological, Environmental and Rural Sciences, Aberystwyth University, Aberystwyth, United Kingdom
- 2Earlham Institute, Norwich, United Kingdom
- 3Department of Horticulture, University of Arkansas, Fayetteville, AR, United States
- 4International Center for Tropical Agriculture, Cali, Colombia
- 5Tropical Forage Program, Alliance Biodiversity-CIAT, Cali, Colombia
- 6School of Agriculture and Environment, Faculty of Science, The University of Western Australia, Crawley, WA, Australia
Brachiaria (Trin.) Griseb. (syn. Urochloa P. Beauv.) is a C4 grass genus belonging to the Panicoideae. Native to Africa, these grasses are now widely grown as forages in tropical areas worldwide and are the subject of intensive breeding, particularly in South America. Tolerance to abiotic stresses such as aluminum and drought are major breeding objectives. In this study, we present the transcriptomic profiling of leaves and roots of three Brachiaria interspecific hybrid genotypes with the onset of water stress, Br12/3659-17 (gt-17), Br12/2360-9 (gt-9), and Br12/3868-18 (gt-18), previously characterized as having good, intermediate and poor tolerance to drought, respectively, in germplasm evaluation programs. RNA was extracted from leaf and root tissue of plants at estimated growing medium water contents (EWC) of 35, 15, and 5%. Differentially expressed genes (DEGs) were compared between different EWCs, 35/15, 15/5, and 35/5 using DESeq2. Overall, the proportions of DEGs enriched in all three genotypes varied in a genotype-dependent manner in relation to EWC comparison, with intermediate and sensitive gt-9 and gt-18 being more similar to each other than to drought tolerant gt-17. More specifically, GO terms relating to carbohydrate and cell wall metabolism in the leaves were enriched by up-regulated DEGs in gt-9 and gt-18, but by down-regulated DEGs in gt-17. Across all genotypes, analysis of DEG enzyme activities indicated an excess of down-regulated putative apoplastic peroxidases in the roots as water stress increased. This suggests that changes in root cell-wall architecture may be an important component of the response to water stress in Brachiaria.
Introduction
For healthy development, growth and reproduction plants need sufficient water. However, due to their sessile nature, plants often encounter unfavorable environmental conditions during their life cycles and water stress is a major environmental factor that limits crop growth and yield. Around one third of the planet is arid to semi-arid, with periodic drought affecting most of the rest of the landmass. As climate changes, more areas are being affected by water stress and for longer periods, and this poses major challenges for global agriculture. Therefore, understanding whole plant and molecular mechanisms that influence responses to water stress is of significant interest to plant scientists and breeders in seeking to maintain and improve crop yields.
Brachiaria (Trin.) Griseb. (syn. Urochloa P.Beauv.) is a C4 grass genus belonging to the Panicoideae (Renvoize et al., 1996). This genus includes several species which are important as agricultural grasses, notably B. decumbens Stapf., B. brizantha (Hochst. ex A. Rich.), B. humidicola (Rendle) Schweick and B. ruziziensis (R. Germ. and C. M. Evrard). These grasses, native to Africa, are now widely grown in the form of individual species and hybrids, as forage grasses in tropical areas worldwide (Keller-Grein et al., 1996). A number of Brachiaria species (referred to collectively as Brachiaria from this point forward) have been the subject of intensive breeding efforts. It is estimated that resultant forage varieties, many of them developed at the International Centre for Tropical Agriculture (CIAT), Colombia, now cover an area of 25 million hectares of agricultural land in Latin America1 and a further 99 million hectares in Brazil (Jank et al., 2014). In addition, opportunities for expanding the use of Brachiaria in Africa and Asia are currently being explored (Maass et al., 2015; Bui Van and Ba, 2019). Particular aims of breeding for these grasses have been to maintain and improve forage quality while increasing tolerances to abiotic stresses such as aluminum (acid soils) and drought, in addition to disease and pest resistance.
Much of the published research on the agriculturally important Brachiaria cultivars has focused either on agronomy and physiological evaluations of trait performance and response to abiotic and biotic stresses and exploring the genetics and cell biology of apomixis. In terms of abiotic stress, one of the reasons for the widespread use of Brachiaria cultivars as forage grasses is that they are often considered to be able to maintain productivity and ground cover under water-limited conditions (Petter et al., 2013; Pizarro et al., 2013; Cheruiyot et al., 2018) and a number of studies have recorded physiological responses of Brachiaria genotypes in response to drought. Observations include that hybrid Brachiaria cultivar Mulato II (B. ruziziensis x B. brizantha) manifests a ‘water saving’ strategy under imposed water-limitation in comparison to Napier grass (Pennisetum purpureum). This strategy includes closure of stomata, leaf rolling and reduced transpiration rates at relatively high soil moisture contents (Cardoso et al., 2015). In addition to these physiological responses, in a comparison of the B. brizantha cultivars Marandu and BRS Piatã, increased production of roots at lower soil levels and increased leaf senescence were also observed in response to water stress (Santos et al., 2013). A further study comparing five different Brachiaria species (B. brizantha, B. decumbens, B. mutica, B. humidicola, and B. dictyoneura) has indicated that differences in overall growth rates, root distributions, osmotic adjustments and timings of stomatal closures could all contribute to variations in drought tolerances (Guenni et al., 2002, 2004). Thus, it is established that there exists appreciable physiological variation in response to water-limitation within Brachiaria, which can be exploited in improving drought tolerance.
Despite the widespread importance of Brachiaria in tropical agriculture, it can still be considered an ‘orphan crop’ in terms of the resources for molecular genetics, biological and genomic analyses. Much of the published work has either focused on identifying the apomixis locus and linked markers which may be useful in plant breeding (Tohme et al., 1996; Pessino et al., 1997, 1998; Zorzatto et al., 2010; Thaikua et al., 2016; Worthington et al., 2016; Worthington et al., 2019) or generating improved understanding of the molecular phylogeny of the species group (recent references include Pessoa-Filho et al., 2015, 2017; Triviño et al., 2017; Kuwi et al., 2018; Tegegn et al., 2019). Again, because of the modest resources available for research on tropical forages, only recently have comprehensive RNAseq-based gene expression studies for Brachiaria been published. These consist of comparative leaf transcriptomics of two highly divergent B. humidicola genotypes (Vigna et al., 2016) and differential gene expression in a B. decumbens genotype exposed to aluminum (Salgado et al., 2017). However, recently, a major landmark in Brachiaria research has been the release of the first draft genome, derived from a diploid B. ruziziensis accession. This release was accompanied by a comprehensive gene annotation and a study of differential gene expression in response to aluminum in both B. ruziziensis and B. decumbens (Worthington et al., 2020).
While many studies on differential gene expression in response to water stress have been published on grass species, particularly the major cereals, there is only a limited body of knowledge on water stress-related changes in patterns of gene expression in forage grasses (Foito et al., 2009; Liu and Jiang, 2010; Pan et al., 2016; Zhao et al., 2016; Zhu et al., 2017; Ji et al., 2018; Dinkins et al., 2019; Fradera-Sola et al., 2019), the majority of which focus on temperate C3 forages. Drought tolerance of breeding selections is routinely evaluated in the CIAT forage breeding program, including measuring water extraction under progressive drying soil conditions – for which variation exists across apomictic Brachiaria germplasm. As described previously, the orphan status of Brachiaria in terms of molecular genetic evaluation means that little is known about gene expression responses linked to the onset of water stress in this genus. Our motivation for the present study was to begin to address this lack of information. Thus, using an experimental system which allowed for progressive sampling of both leaves and roots in a drying growing medium, we have undertaken a bioinformatic comparison of patterns of differential gene expression between three Brachiaria genotypes with different tolerance to water stress in CIAT evaluations, in order to examine the nature and conservation of these responses.
Materials and Methods
Plant Material and Growth Conditions
Three Brachiaria hybrid breeding selections developed at CIAT, which have shown contrasting responses to an imposed drought condition, were used in this study (Supplementary Table 1). These included Br12/3659-17 (gt-17), Br12/2360-9 (gt-9), and Br12/3868-18 (gt-18) previously characterized as having good, intermediate and poor tolerance to drought, respectively, in evaluations conducted in the CIAT breeding program. Plants and seeds used in this study were obtained directly from the CIAT tropical forage breeding program.
The three hybrids used were developed from an interspecific recurrent selection program focused on developing improved apomictic Brachiaria cultivars by crossing a synthetic, fully sexual breeding population with a non-inbred apomictic tester (Miles, 2007). The recurrent selection population was developed by crossing a sexually reproducing synthetic autotetraploid accession of B. ruziziensis with nine apomictic tetraploid accessions of B. decumbens and B. brizantha and recombining their sexual progeny during nine cycles of open pollination between 1992 and 2011. All the wild accessions used to initiate the CIAT Brachiaria breeding program are publicly available in the CIAT Genebank2. The three breeding selections chosen for this experiment are apomictic progeny of sexual selections from the ninth cycle of recurrent selection with the apomictic tester B. decumbens CIAT 606 (cv. Basilisk).
One seed per accession was sown in John Innes potting compost and germinated and grown in a climate-controlled cabinet under a photoperiod of 12 h per day with 60% relative humidity at 25°C constant temperatures. When the plants were large enough, they were split and left to grow in the control cabinet. After establishment, the four most similar looking clones from each accession were split again into three further clones, and nine clones per accession were chosen for the experiment (three sampling points and three replicates per sampling point). The roots of the nine clones were cleaned of soil, placed in vermiculite, and watered with a nutrient solution consisting of: [μM], NH4NO3 [500], KNO3 [300], Ca[NO3] [200], NaH2PO4 [5], MgCl2 [90], MgSO4 [60], FeCl3 [5], Na-EDTA [5], H3B03 [6], MnSO4 [1], ZnSO4 [1], CuSO4 [0.2], Na2MoO4 [1], Na2SiO3 [5], NaCl [55], adjusted to pH 4.2 using 1M HCl. This nutrient solution, mimicking the pH of acid soils, was developed at CIAT to encourage root growth in Brachiaria (based on Wenzl et al., 2003, 2006). When the plants had recovered from transplantation, they were watered to full capacity and then, with no further watering, the fall in estimated growing medium water content (EWC) was measured using a HH2 Delta-T meter (AT Delta-T devices, Cambridge, United Kingdom). Leaf and root material were sampled at 35% (full capacity), 15% and 5% EWC (c. 15 days). Leaves were sampled directly into liquid nitrogen and root material was washed briefly in distilled water to remove vermiculite and then placed in liquid nitrogen and stored at –80°C prior to extraction of RNA.
Leaf Relative Water Content
Estimations of leaf relative water content (RWC) were measured at 35, 15, 5, and 1% EWC (RNA extraction was not carried out at 1%). RWC estimations were carried out on additional replicates of the same genotypes prepared and grown identically to the plants used for RNA extraction, as follows. Three leaves from each replicate were removed, an 8 cm mid-section was cut from each leaf, and the fresh weight (FW) was measured. This excised section was then placed in a 50 ml tube containing 5 ml water, capped and left at 4°C for 24 h. After this period, the leaf sections were blotted and turgid weight (TW) was measured. The sections were then dried for 24 h at 80°C for the dry weight (DW). RWC was calculated as (FW-DW/TW-DW) × 100.
RNA Extraction and cDNA Library Construction and Sequencing
RNA extraction, library construction and sequencing was carried out on each replicate independently (i.e., tissues and samples were not pooled). RNA was extracted from the root and leaf material, using a hot phenol technique (Ougham and Davies, 1990), and suspended in 100 μl of sterile 0.5M TE buffer. Quantification of the RNA samples was performed on a NanoDrop® 1000 (Thermo, Waltham, MA, United States) at a wavelength of 260 nm. RNA sample quality was evaluated using the A260 nm/A280 nm wavelength ratio and by direct observation on a 1% agarose gel. Sequencing libraries were constructed and sequenced at the DNA sequencing Center, Brigham Young University, Provo, UT, United States using Illumina kits for either Poly-A selected or Ribo-Zero rRNA removal. Illumina sequencing was performed using a HiSeqTM 2000 platform according to the manufacturer’s instructions (Illumina, San Diego, CA, United States).
RNAseq Processing and Quality Control and Mapping
Prior to mapping, raw reads were processed using Trimmomatic v.0.33 (Bolger et al., 2014) to remove adapters using the following parameters (optimized after several run tests): ILLUMINACLIP:TruSeq3-PE-2.fa LEADING:15 SLIDINGWINDOW:4:15 MINLEN:30 HEADCROP:12, and the quality of resulting trimmed and cleaned reads was assessed using FastQC v.0.11 (Wingett and Andrews, 2018). Reads were then mapped to the assembly version of the Brachiaria genome (Worthington et al., 2020) using the splice-aware mapper Hitsat2 v.2.0.0 (Kim et al., 2015).
Pre-Processing and Quantification of Transcripts
Prior to calling of differentially expressed genes (DEGs) a pre-processing filtering was performed to remove potential artifacts and assess the quality of the replicates. Count matrices were derived from bam files above using the GenomicFeatures and GenomicAlignments R libraries. Transcripts with a count lower than one in any sample were discarded. We applied the regularized logarithm transformation (rlog) as implemented in the DESeq2 package to decrease the variance among gene expression values (Love et al., 2014) and then calculated a distance matrix between samples and performed a principle component analysis (PCA) to quantify experimental covariates and batch effects among samples and replicates (Fisher et al., 2004).
Estimating the Completeness of Transcriptomes
The transcriptome in each sample was assessed for its completeness as a measure of quality of the sequencing. Clean reads were mapped to the reference genome and assembled using StringTie v1.1.0 (Pertea et al., 2015) using default parameters. The completeness of each transcriptome was assessed using BUSCO (Simao et al., 2015) on the early_release plantdb set, composed of 1440 core genes.
Identification of Differentially Expressed Genes
Quantification of transcripts was done using Salmon (Patro et al., 2017) using precomputed mapping files (bam files) generated as described above using the –ValidateMappings –gcBias and –numBootstraps set to 1000 to improve the quantification. Derived counts were used as inputs to call DEGs using DESeq2 across three EWC-point comparisons: 35% vs. 15% (35/15), 35% vs. 5% (35/5) and 15% vs. 5% (15/5) for the three genotypes independently both in shoot and root tissues. Genes with a log2 fold change (LFC) above one and a false discovery rate (FDR) of ≤5% were considered as DEGs.
Differentially expressed genes were categorized according to pattern of expression and up- or down-regulation over the three sampling point comparisons and assigned an expression category. Each expression category was defined by a three-letter code: the first letter indicates whether the DEGs contributing to the enrichment of the indicated GO terms were up- (u) or down- (d) regulated at comparison point 35/15, the second letter at 35/5 and the third letter at 15/5. The letter n at position one, two or three indicates that the GO term was not significantly enriched at that comparison point. E.g., GO terms in expression category d-d-n were associated with down-regulated DEGs at 35/15 and 35/5 but were not-significantly enriched at 15/5 by either up- or down-regulated DEGs; GO terms in expression category n-n-u were associated with up-regulated DEGs at 15/5 but were not-significantly enriched at either 35/15 or 35/5.
Functional Annotation of Differentially Expressed Genes
The reference genome was functionally re-annotated using Blast2GO 5.25 (Pro) (Conesa and Gotz, 2008) as a prior step before computing Gene Ontology (GO) term enrichments. The functional annotation was done as follows: BLAST searches were performed on the nr database (release March 2019) using BLASTx command from ncbi-blast-2.2.28 + release (Camacho et al., 2009) at an E-value cut-off of 1 × 10–6 and selecting the top 20 hits. InterPro searches were performed using InterProScan v.5.18-57 (Jones et al., 2014) on TIGRFAM (Haft et al., 2013), PFAM (Finn et al., 2016), SMART (Letunic et al., 2006), PANTHER (Mi et al., 2013), and Gene3d (Lees et al., 2010) databases.
Identification of GO Terms and Mapping of Enzyme Codes to KEGG Pathways
Gene Ontology term enrichment was calculated for each set of DEGs associated with each three letter expression category (see section “Identification of differentially expressed genes”) for each genotype, using Blast2GO 5.25 (Pro) (Conesa and Gotz, 2008) and a 5% FDR cut-off threshold. GO terms were subsequently grouped into putative functional hierarchies visualized at https://www.ebi.ac.uk/QuickGO/slimming.
The complete lists of GO terms were further filtered according to either of two criteria in order to focus on the best-supported enrichments: (1) GO terms which were associated with differential gene expression in all three genotypes, or (2) GO terms that were within the top 10% of most significantly enriched GO terms within the original 5% FDR (Fisher’s exact test p-value cut-off thresholds of –log10 p = 10.4 and p = 8 for leaf and root, respectively).
Mapping of enzyme codes to KEGG pathways (Kanehisa and Goto, 2000) was accomplished using the relevant module contained within Blast2GO.
Results
Relative Water Content
Supplementary Figure 1 illustrates the change in RWC for the three hybrid Brachiaria genotypes as the EWC of the medium decreased. Major changes in RWC for all 3 genotypes only occurred after the 5% EWC point was reached. No significant changes in RWC were observed with decreasing EWC between 35 and 5% for gt-9 and gt-18, however, gt-17 did show a significant decrease in RWC over the same range (P < 0.01). These results indicated the onset of water-stress for the genotypes at around 5% EWC.
Pre-processing, Mapping, and Quality of Sequencing and Replicates
RNA libraries were processed and sequenced in a single batch yielding an average of c. 11M reads per replicate for both leaf and root tissue, with a maximum and minimum of 13.2 and 10.1M reads for any individual replicate. The drop-off rate upon trimming and quality control (discarding low quality and non-paired reads) was less than 1%. The mapping rate of the retained reads to the reference genome (Worthington et al., 2020) varied between 76 and 58% for leaf tissue and 72 and 44% for the root tissue (Supplementary Table 2). From the BUSCO analysis, the completeness of the transcriptome across all samples was estimated at an average of c. 77% in terms of complete and partial gene coverage. This compares with 85% for the Brachiaria reference genome (Worthington et al., 2020) (Supplementary Table 3).
The variability of gene expression based on normalized counts among the replicates is illustrated in the form of density and principal component analysis (PCA) plots (Supplementary Figures 2, 3). The density plots of the rlog transformation of normalized counts both in leaf and root samples were very homogenous and, thus, indicate little variability among replicates (Supplementary Figures 2A–F). This is also reflected in the PCA plots (Supplementary Figure 3), which show that within each EWC category the different replicates tended to cluster together, with the first two principle components accounting for 62–87% of the total variability. The exception to this was for gt17, which showed less tight clustering of the replicates at 15% EWC (leaf; Supplementary Figure 3B) and 35% EWC (root; Supplementary Figure 3E). In terms of the first principle component, the analysis showed that the distance between samples was greatest for the 35 and 5% EWC sampling points and that the 15% sampling point was intermediate between 35 and 5%. It would appear, therefore, that the change in EWC was the main driver explaining the overall variability among samples.
Distribution of Differentially Expressed Genes Over Genotypes and Comparison Points
Leaf and root transcriptomes were analyzed for the presence of DEGs across three EWC-point comparisons, 35/15, 35/5, and 15/5 for the three genotypes independently. The total numbers of up- and down-regulated DEGs at any stage and at each of the comparison points is given in Table 1. On average across genotypes at any time, 2,898 and 2,595 DEGs were identified in leaves and root, respectively, evenly distributed between up- and down-regulated DEGs. The total number of DEGs associated with the individual genotypes was quite variable with c. twice as many DEGs identified in the leaf tissue for gt-18 as compared to gt-9 and c. three-times as many DEGs identified in gt-9 as compared to gt-17 in the root tissue. For the comparison points, the highest numbers of DEGs were associated with the 35/5 comparison in both leaves and roots; and the smallest number with the 35/15 comparison in leaves and the 15/5 comparison in roots. Out of the total of 35,196 gene models included in the analysis, 5,821 were significantly differentially up- or down-regulated in any of the genotypes at any of the comparison points in the leaf tissue and 5,322 in the root tissue, with 2,118 specific gene models being present as DEGs in both leaves and roots (Supplementary Figure 4). The relative distributions of DEGs between and among the three genotypes in leaves and roots is described in Supplementary Figure 5 and the complete lists of DEGs associated with genotypes and comparison points for both leaf and root is given in Supplementary File 1.
Table 1. The total number of up- and down-regulated differentially expressed genes at all comparison points and at the individual comparison points for each genotype and tissue.
The overall patterns of differential gene expression varied between genotypes and between leaves and roots (Figure 1). In the leaves, while all the genotypes were similar in having high proportions of their DEGs present in just the 35/5 comparison stage, differences between the genotypes were also apparent. For example, for gt-9 and gt-18 between c. 10 and 22% of their up- and down-regulated DEGs were detected in both the 35/15 and 35/5 comparison stages, whereas the equivalent figures for gt-17 were < 5%. Also, almost 50% of the DEGs for gt-17 were identified only in the 15/5 comparison stage; for gt-9 and gt-18 the equivalent figures were <15% and <5%, respectively. For roots, gt-9 and gt-18 showed more even distributions of DEGs across all of the comparison stages. However, this was in contrast to gt-17, in which 72% of the DEGs were present only in the 35/5 comparison stage. In summary, patterns of differential gene expression were similar in gt-9 and gt-18 and this differed from the pattern seen in gt-17. Overall, for all genotypes, the patterns of differential expression were different between leaves and roots.
Figure 1. Relative proportions of differentially expressed genes (DEGs) across genotypes (gt-) 9, 17 and 18 and combined estimated water content (EWC) comparison points (expression patterns). The percentage proportions of up- and down-regulated DEGs identified for each individual expression pattern are indicated by the colored columns. Horizontal bars beneath the x axis indicate the individual EWC comparison points (indicated at the end of the x axis) in which the genes were differentially expressed for each expression pattern; red, differentially expressed for that EWC comparison point; gray, not differentially expressed for that EWC comparison point. E.g., the first expression pattern indicates the percentage proportion of the total number of up- and down-regulated DEGs for each genotype that were differentially expressed for the 35% vs. 15% comparison point (red) but were not differentially expressed for the 35% vs. 5% and 15% vs. 5% comparison points (gray).
Association of DEGs With Enriched GO Terms
The DEGs identified within each category were analyzed for association with GO terms on an individual genotype basis using a 5% FDR. Across the three genotypes, a total of 1,210 significantly enriched GO terms (referred to as GO terms from this point forward) were identified from the leaf DEGs and 856 from the root DEGs. For the individual genotypes the total numbers of GO terms were 280, 617, and 722 for leaves and 420, 212 and 467 for roots for gt-9, gt-17 and gt-18, respectively (a complete lists of GO terms for individual genotypes for leaf and root are provided in Supplementary Files 2, 3 respectively). Of the total numbers of GO terms associated with leaf tissue, 137 and 85 were enriched for all three genotypes in leaves and roots, respectively. The number of DEGs identified for each comparison stage, the number of associated enriched GO terms and the percentage proportion of the DEGs which contribute to these enriched GO terms is summarized in Table 2.
Table 2. Summary descriptions of the numbers and proportions of differentially expressed genes (DEGS) according to genotype (gt) and expression category.
GO Terms Associated With Leaf DEGs
Of the 137 GO terms associated with leaf transcriptomes of all three genotypes, 67 could be assigned to two main GO term hierarchies relating to chloroplast/photosynthetic metabolism, and carbohydrate/cell wall metabolism (Groups 2 and 8 in Supplementary File 4). For chloroplast/photosynthetic metabolism (25 Cellular Compartment [CC] and 9 Biological Process [BP] GO terms), the vast majority of GO terms were associated with down-regulated DEGs. Gt-9 and gt-18 started showing significant GO terms association with down-regulation during the earliest comparison (35/15), i.e., d-d-n and d-d-d, whereas for gt-17 significant association of GO terms with down-regulation was not detected until the 15/5 comparison (n-d-n and n-d-d). Thus, as a trend, down-regulation of the genes associated with the chloroplast/photosynthetic metabolism GO terms occurred earlier in gt-9 and gt-18 than in gt-17.
The second major hierarchy, consisting of carbohydrate/cell wall metabolism-related GO terms (24 BP, 8 Molecular Function [MF] and 1 CC GO terms) showed a major difference in the direction of regulation of DEGs associated with the same GO terms between genotypes. The great majority of GO terms for gt-9 and gt-18 were associated with up-regulated expression categories (mainly u-u-n and n-u-n) whilst all of the gt-17 GO terms were associated only with down-regulated expression categories (n-d-n and n-d-d), indicating major differences in the associated metabolic processes between genotypes at these point comparisons. Of the remaining GO categories, focusing on more specific terms, down-regulated categories were associated with carotenoid, cysteine and pyruvate metabolism, nitrate assimilation and response to light stimulus as well as glyoxysome and stromule CC terms. Up and down regulated GO categories were associated with alpha-amino acid, carboxylic acid and malonyl CoA biosynthetic metabolic processes. Only two GO terms were exclusively up-regulated, associated with the hexosamine pathway.
Besides the terms described above, a further 63 GO terms were not detected in all three genotypes but had highly significant p-values (10% most significant p-values within the 5% FDR) in at least one of the genotypes (Supplementary File 6). Nine of these were exclusively down-regulated in gt-9 and gt-18 and could be associated with chloroplast/photosynthetic metabolism and a further 12 were up- and down-regulated in gt-17 and gt-18 and were related to organelle compartment and organo-nitrogen/phosphate and carbohydrate metabolic processes. The remaining 42 GO terms were exclusively up-regulated, all were present in gt-17. Two were also present in gt-18 though these were not within the most significant 10% of p-values. These latter GO terms were associated with ribosome metabolism and location, translation and nucleotide binding. All were present in expression category n-n-u and 14 were also identified in n-u-u.
GO Terms Associated With Root DEGs
A total of 85 GO terms were associated with all three genotypes from the root data, with the largest GO hierarchy consisting of nine GO terms (Supplementary File 5). Again, there was a noticeable difference between genotypes in terms of the expression categories in which the majority of GO terms were represented, with gt-9 and gt-18 being associated with GO terms across a number of the expression categories and gt-17 being associated with GO terms, predominantly, from only two expression categories, n-u-n and n-d-n. Thirty-seven GO terms were mostly down-regulated and fell into six main GO hierarchies with a single unconnected GO term. These were associated with peroxide activities, cellular detoxification, nitrate, phosphate and serine family amino acid metabolism and plant cell walls. Thirty-eight GO terms showed a degree of up- and down-regulation falling into 10 small GO hierarchies with five unconnected GO terms. The most specific GO descriptions were associated with membrane transport, galactose and glutamine family amino acid metabolism and metal ion binding. The remaining eight GO terms were exclusively up-regulated and fell into two small hierarchies relating to xanthine catabolism and inosine monophosphate (IMP) salvage.
A further 63 GO terms were not detected in all three genotypes but were within the 10% most significant p values within the 5% FDR (Supplementary File 7). Forty-eight of these were primarily down regulated and contained within eight GO hierarchies and two unconnected GO terms. Six of the down regulated GO hierarchies were from gt-9, in expression categories d-n-n and d-d-n and were associated with nucleotide, energy and amino acid metabolism and a single GO term with nitrate transport. A further two down-regulated hierarchies were in gt-18 and associated with the cytoskeleton and the terpenoid biosynthetic process, with an additional unconnected GO term for DNA replication. Seven primarily up-regulated GO terms in two hierarchies were in gt-9 and associated with glucan and beta-glucosidase activity and a further four GO terms in gt-17 were associated with transmembrane transport.
KEGG Metabolic Pathways and Enzyme Codes
Differentially expressed genes associated with enzyme codes were mapped onto KEGG metabolic pathways. A total of 441 enzyme codes could be mapped to 134 pathways across all three genotypes, though with very uneven distribution (Supplementary Files 8, 9 for KEGG pathways and enzymes respectively). The five most frequently occurring pathways (excluding biosynthesis of antibiotics) were starch and sucrose metabolism, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, purine metabolism and galactose metabolism (Table 3). In terms of overall trends, DEGs contributing to starch and sucrose metabolism, galactose metabolism and glutathione metabolism were up-regulated markedly more frequently than down-regulated and DEGs contributing to phenylpropanoid metabolism, glycolysis/gluconeogenesis and cysteine and methionine metabolism showed the opposite trend. However, these overall figures can hide differences between the genotypes, particularly in relation to differences between gt-17 and gts-9/18 and DEGs from the starch and sucrose metabolism pathway; aligning with the observations from enriched GO terms, enzyme activities associated with carbohydrate and cell wall metabolism were more frequently down-regulated in leaves of gt-17 and up-regulated in gts-9/18. This can be seen specifically in terms of the most frequently occurring enzyme code associated within the starch and sucrose metabolism pathway, ec:2.4.1.12 [cellulose synthase (UDP-forming); Table 3].
Table 3. KEGG pathways and enzyme codes associated with up- and down-regulated differentially expressed genes (DEGs) responding to increasing water stress in Brachiaria.
The single most frequently occurring enzyme code was ec:1.11.1.7 (peroxidase; phenylpropanoid pathway), which was represented c. twice as often as the second most frequently occurring enzyme code, ec:2.5.1.18 (glutathione transferase) and showed no obvious difference between genotypes. As well as occurring most frequently, ec:1.11.1.7 was notable in that it occurred c. three times more often in the root as compared to the shoot and was down-regulated c. 2.5 times more often than it was up-regulated. Both ec:1.11.1.7 and ec:2.5.1.18-type activities would be predicted to be associated with the response to the presence of reactive oxygen species though with activities likely to be manifested in different cellular locations (apoplastic and intracellular, respectively; Dixon et al., 2009; Podgorska et al., 2017).
Discussion
In this study, we have compared the gene expression, in response to the onset of water stress, of three Brachiaria hybrid genotypes that have previously been evaluated in terms of ability to extract water from progressively drying soil in a container-based assay. In that study gt-17 showed the greatest, gt-9 intermediate and gt-18 the least water extraction. In the experiments reported here, our focus has not been directly to suggest specific up- or down-regulated genes which have significant effects on these or other physiological responses to increasing water stress and discriminate the genotypes under study, but more to describe overall gene-expression responses, interpreted through GO and pathway analyses, in both leaves and roots. Clearly, experimental designs of this kind are comparing widely spaced ‘snapshots’ of gene expression at assay points rather than approximating continuous monitoring. Also, identified differences between tissues and genotypes do not necessarily indicate overall presences or absences of particular kinds of biological responses, just their detections at defined assay points. However, the outcomes of these studies can be valuable for indicating variability between genotypes and so point to areas which may be fruitful for further research focus in the contexts of both plant biology and, particularly in the context of this study, the exploitation of Brachiaria genetic variability for forage grass improvement.
Drought Tolerant gt-17 Hybrid Has Distinctive Patterns of Gene Expression and GO Terms
Differentially expressed genes were identified by quantitative comparisons at three stages, 35/15, 35/5, and 15/5, representing significant and progressive changes in gene-expression profiles relating to the early stages of response to water stress. Thus, gene-models could be characterized as to their differential expression into six categories (Figure 1). From this analysis, it was apparent that overall patterns of gene expression differed both between leaves and roots and between genotypes. In particular, in the leaves, drought tolerant gt-17 tended to have more DEGs at lower EWCs relative to the genotypes with less drought tolerance, gt-9 and gt-18. For the roots, the pattern of appearance of DEGs was also distinctive for gt-17. In this latter case, the majority of DEGs were detected during the 35/5 comparison, indicating a preponderance of differential expression patterns which only reached significance over the whole analysis period, i.e., not significant within the 35/15 (early) and 15/5 (late) periods when considered separately (Figure 1).
The overall trends seen in Figure 1 can be compared with those generated when just the DEGs which contribute to the enrichment of GO terms are included (c. 70% of the total DEGs for leaves and c. 62% of the total DEGs for roots; Figure 2). The trends are very similar when considering all DEGs and only those DEGS which contribute to the significant enrichment of GO terms. Thus, the DEGs at the different comparison stages are contributing proportionally to biological processes as indicated through GO terms. However, while the relative proportions across comparison periods were consistent, there was some variation in terms of genotype as indicated in Table 4. Particularly, (a) proportionally fewer of the total number of up- and down-regulated leaf DEGs in gt-9 contribute to the enrichment of GO terms compared to gt-17 and gt-18 and (b) for gt-18, only 34% of the up-regulated root DEGs contribute to GO term enrichment whereas the equivalent figure for down-regulated DEGs was 72%. The equivalent comparisons for gt-9 and gt-17 were far more even.
Figure 2. Relative proportions of differentially expressed genes (DEGs) contributing to the enrichment of GO terms across genotypes (gt-) 9, 17, and 18 and combined estimated water content (EWC) comparison points (expression patterns). The percentage proportions of up- and down-regulated DEGs identified for each individual expression pattern are indicated by the colored columns. Horizontal bars beneath the x axis indicate the individual EWC comparison points (indicated at the end of the x axis) in which the genes were differentially expressed for each expression pattern; red = differentially expressed for that EWC comparison point; gray = not differentially expressed for that EWC comparison point. E.g., the first expression pattern indicates the percentage proportion of the total number of up- and down-regulated DEGs for each genotype that were differentially expressed for the 35% vs. 15% comparison point (red) but were not differentially expressed for the 35% vs. 5% and 15% vs. 5% comparison points (gray).
Table 4. The relative proportions of differentially expressed genes (DEGs) which contributed to the significant enrichment of GO terms according to genotype, organ and direction of regulation.
Overall differences in differential gene expression itself and the contribution of DEGs to GO terms indicate that there can be substantial variation between genotypes and tissues. This genotype-dependent manifestation of the variation in the overall patterns of DEGs and GO terms can be seen when looking at particular biological examples, as illustrated in Figure 3. The data shown focuses on the GO terms predominantly associated with down-regulation relating to chloroplast/photosynthetic metabolism in the leaves and peroxidase metabolic processes and cellular detoxification in the roots. While the responses are not identical in leaves and roots, the genotype-specific trend is the same with the down-regulation of the indicated biological process commencing (i.e., reaching significance in terms of enrichment) earlier for gt-9 and gt-18 than for gt-17. In fact, when all GO terms are considered (Supplementary Files 2, 3) a total of 1611 enriched GO terms are associated with DEGs up- or down-regulated in the 35/15 comparison period for gts-9 and gt-18 but only 46 for gt-17.
Figure 3. GO terms identified in all three genotypes relating to photosynthetic metabolism in the leaves and peroxide and cellular detoxification processes in the root. GO Terms (A) relative positions of GO terms linked through hierarchies visualized at https://www.ebi.ac.uk/QuickGO/slimming. (B) GO category, BP, biological process; CC, cell compartment; MF, molecular function. (C) GO description (1MP, metabolic process; 2CP, catabolic process). Expression Categories: Expression patterns of differentially expressed genes contributing to enrichment of GO terms. d, down-regulated; u, up-regulated; n, not significantly down- or up-regulated; x-y-z, d, u, or n relating to the 35/15 (x), 35/5 (y), and 15/5 (z) sampling point comparisons. The first column for each expression category indicates whether a GO term(s) was significantly enriched for gt-9 (blue), the second column for gt-17 (orange) and the third for gt-18 (green).
An additional example of the distinct responses of gt-17 is indicated by the leaf GO hierarchies for which the terminal GO terms were GO:0030244 (cellulose biosynthetic process) and GO:0016760 (cellulose synthase [UDP-forming]) (Figure 4). These were predominantly associated with up-regulation at the first comparison point (35/15) for gt-9 and gt-18 and down-regulation for gt-17 over the 35/5 comparison point; for most of these gt-17 GO terms, there was association with further down-regulation during the 15/5 comparison point. The GO terms associated with up-regulation in this hierarchy for gt-9 (69 DEGs) and gt-18 (217 DEGs) were enriched by a total of 234 individual gene models, 52 of which were common to both genotypes (representing 75% of the total for gt-9). The down-regulated gene models for gt-17 for this hierarchy contained 154 gene models, 23 of which were also in the up-regulated DEGs for gt-9 and gt-17. Thus, the DEGs being differentially regulated in opposite directions in gts-9/18 and in gt-17 had some overlap but were, in the main, different subsets of genes relating to the same GO terms.
Figure 4. GO terms identified in all three genotypes relating to carbohydrate and cell wall metabolism in the leaf. Leaf GO Terms (A) relative positions of GO terms linked through hierarchies visualized at https://www.ebi.ac.uk/QuickGO/slimming. (B) GO category, BP, biological process; CC, cell compartment; MF, molecular function. (C) GO description (MP, metabolic process; CP, catabolic process; BsP, biosynthetic process; act., activity). Expression Categories: Expression patterns of differentially expressed genes contributing to enrichment of GO terms. d, down-regulated; u, up-regulated; n, not significantly down- or up-regulated; x-y-z, d, u or n relating to the 35/15 (x), 35/5 (y), and 15/5 (z) sampling point comparisons. The first column for each expression category indicates whether a GO term(s) was significantly enriched for gt-9 (blue), the second column for gt-17 (orange) and the third for gt-18 (green).
KEGG pathway analysis of these sets of up- (gts-9/18) and down- (gt-17) regulated gene models mapped onto 16 and 19 enzyme codes, respectively, in the Starch and Sucrose Metabolism pathway including 13 gene models in both of gts-9/18 and gt-17 assigned to ec:2.4.1.12 (cellulose synthase [UDP-forming]) but regulated in the opposite direction (Supplementary Figure 6). Additionally, the down-regulated gt-17 gene model set contained nine annotations describing ‘starch synthase’ and four describing ‘sucrose-phosphate synthase’, indicating down-regulation of key enzymes of energy metabolism, additional to those more directly involved in cellulose synthesis. Thus, while in gt-17 GO terms relating to photosynthesis and carbohydrate metabolism are all associated with down-regulation during the same comparison points, for gt-9 and gt-18, the association with down-regulation of photosynthesis-related GO terms is accompanied by association with up-regulation of the cell-wall and carbohydrate-related GO terms (Figures 3, 4). Again, this suggests underlying differences between gts-9/18 and gt-17 in their metabolic responses to the progression of water limitation. Similarly, one particular hierarchy of GO terms, relating to ribosomal metabolism, was present for gt-17 but completely absent for gts-9/18 (Figure 5) with significant enrichment for all GO terms in the 15/5 comparison period and for five of the terms also in the 35/5 comparison period. A total of 313 up-regulated gene models contributed to the enrichment of these GO terms, of which 82 were annotated as ribosomal proteins. There was also up-regulation of 15 ribosomal protein genes from gt-18 in the u-n-n, u-n-u and u-u-n expression categories, though this was not sufficient for significant enrichment of the related GO terms. Only the ribosomal protein genes were DEGs from gt-9 and all were down-regulated. The basic function of ribosomal proteins is protein synthesis as part of the ribosomal complex. However, the potential role of ribosomal proteins in conferring abiotic stress, and particularly drought tolerance, has been reported (Moin et al., 2017). Thus, this represents another interesting contrast between gt-17 and gts-9/18.
Figure 5. GO terms identified relating to ribosomal metabolism and location in the leaf. Leaf GO Terms (A) relative positions of GO terms linked through hierarchies visualized at https://www.ebi.ac.uk/QuickGO/slimming. (B) GO category, BP, biological process; CC, cell compartment; MF, molecular function. (C) GO description (MP, metabolic process; BsP, biosynthetic process). Expression Categories: Expression patterns of differentially expressed genes contributing to enrichment of GO terms. d, down-regulated; u, up-regulated; n, not significantly down- or up-regulated; x-y-z, d, u, or n relating to the 35/15 (x), 35/5 (y), and 15/5 (z) sampling point comparisons. The first column for each expression category indicates whether a GO term(s) was significantly enriched for gt-9 (blue), the second column for gt-17 (orange) and the third for gt-18 (green). Only gt-17 had significantly enriched GO terms.
As referred to earlier, it has been observed that gt-17 can be more drought tolerant than gt-9 or gt-18 in container-based trials conducted as a part of the CIAT breeding program. Clearly, extrapolating from the present growth-room/artificial-media based experiment to a container-based trial, let alone to a field situation has to be done with extreme caution. However, the results reported here do indicate that genetic variation between the three genotypes has significant and measurable impacts on patterns of gene expression in response to water stress. The fact that such variation should exist is not surprising and there are numerous examples of gene-expression studies and QTL analyses for different aspects of drought tolerance in forage and other grasses (e.g., Alm et al., 2011; Merewitz et al., 2014; Baldoni et al., 2016; Jiang et al., 2017; Liang et al., 2017; Zhu et al., 2017; Cheruiyot et al., 2018; Hu et al., 2018; Ji et al., 2018; Kumar et al., 2018; Ling et al., 2018; Ma et al., 2018; Chaichi et al., 2019; Dinkins et al., 2019; Fradera-Sola et al., 2019). However, in the context of Brachiaria breeding, these genotypes were selected from a cultivar development program and the results suggest selectable variation with measurable impacts.
EC:1.11.1.7 – Peroxidase (Phenylpropanoid Biosynthesis Pathway) Down-Regulated Genes Are Over-Represented in the Root in All Genotypes
One interesting observation which was not genotype specific, was that gene models which could be assigned to ec:1.11.1.7 (peroxidase) were strongly over-represented in down-regulated DEGs in the root, in comparison to other enzyme classes as a whole (Table 3 and Supplementary File 9) and other enzyme classes within the phenylpropanoid biosynthesis pathway (Table 5 and Supplementary Figure 7). In a similar study (Fradera-Sola et al., 2019) focusing on the C3 forage grass perennial ryegrass (Lolium perenne), the same enzyme class was also over-represented in down-regulated expression categories in root tissue and, to a lesser extent, in the leaves. Peroxidase genes and activities (though not necessarily ec:1.11.1.7) have been reported as being differentially expressed in a number of studies (Csiszár et al., 2012; Ranjan et al., 2012; Hu et al., 2018; Chaichi et al., 2019) pertaining to drought stress, in some cases suggesting that relatively higher levels of peroxide transcripts are detected in more drought tolerant genotypes under drought stress (Merewitz et al., 2014). Peroxidases play pivotal roles in both cytosolic and apoplastic responses to reactive oxygen species and, in the apoplast, specifically lignification and elasticity of secondary cell walls (reviewed in Podgorska et al., 2017; Meents et al., 2018). In the present study it is particularly the down-regulation of apoplastic peroxidases that is implicated which might suggest weakening of the mechanical properties of the root cell wall as water stress increases. This may allow newly divided cells to expand more readily under water-limited conditions (Tenhaken, 2015), but may also reflect that root cell walls more generally (rather than just at the expansion zone at the root tip) may develop different mechanical properties in response to the drying of the medium in which they are growing.
Table 5. The number of differentially expressed genes (DEGs) associated with enzyme codes from the KEGG Phenylpropanoid Biosynthesis pathway according to genotype, organ and direction of regulation.
Conclusion
We have undertaken a study comparing differential gene-expression, in relation to the onset of water stress, across three Brachiaria hybrid genotypes with contrasting physiological responses to water deficit and identified a range of GO terms linked to biological responses in leaves and roots. In leaves, GO term hierarchies relating to photosynthetic metabolism and carbohydrate metabolism were well supported as were antioxidant and peroxidase activities in roots. The two most striking aspects of the results were, (a) both overall proportions of DEGs and specific related GO terms enriched in all three genotypes varied in a genotype-dependent manner in relation to expression categories, with gt-9 (intermediate water stress tolerance) and gt-18 (susceptible to water stress) being more similar to each other than gt-17 (tolerant of water stress), and (b) the GO terms relating to carbohydrate and cell wall metabolism in the leaves were enriched by up-regulated DEGs for gt-9 and gt-18, but down-regulated DEGs for gt-17 at the EWC points assayed. Additionally, across all genotypes, analysis of enzyme activities relating to DEGs using the KEGG database, indicated a large excess of down-, as compared to up-regulated DEGs with likely apoplastic peroxidase activities in the roots as water stress increased. This suggests that changes in root cell-wall architecture may be an important component of the response to water stress in these genotypes. In taking this work forward, it will be particularly interesting to see if overall time-course patterns of gene expression can be correlated with response to water stress in a wider range of Brachiaria genotypes. If so, it may be possible to identify the genetic basis by which the syndrome of metabolic events, which manifest with increasing water stress are initiated and enacted.
Data Availability Statement
The data presented in this study can be found in the European Nucleotide Archive repository under the study accession number PRJEB41722 (https://www.ebi.ac.uk/ena/browser/view/PRJEB41722). Additional data that support the results are included within the article and its additional files. Other relevant materials are available from the corresponding authors on reasonable request. Plant material is publicly available at the CIAT Genebank (genebanks.org/resources/crops/forages-grass/).
Author Contributions
JAC, JC, MW, and JA conducted the initial experiments to identify breeding selections with good, intermediate, and poor tolerance to drought to be used in this study. NF-F and IA conceived and designed the experiments and computational analyses. CJ, AT, DG, JH, YF, and NF-F performed the experiments. JDV, MW, JD, MB, FC, JA, JAC, and JC contributed expertise and reagents. CJ, YF, AT, IA, and NF-F analyzed the data. CJ, NF-F, and IA drafted the initial version. NF-F and IA finalized the drafting with help from all authors. All the authors have read and approved this manuscript.
Funding
This work was funded by a UK BBSRC Newton Fund Postdoctoral Twinning Award BBS/OS/NW/000011 and UK BBSRC BBS/E/W/0012844A.
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.
Acknowledgments
We would like to thank Colin Sauze, Aberystwyth University, for his assistance with the computing infrastructure.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.637956/full#supplementary-material
Footnotes
- ^ https://ciat.cgiar.org/what-we-do/forages-livestock/
- ^ https://www.genebanks.org/resources/crops/forages-grass/
References
Alm, V., Busso, C. S., Ergon, A., Rudi, H., Larsen, A., Humphreys, M. W., et al. (2011). QTL analyses and comparative genetic mapping of frost tolerance, winter survival and drought tolerance in meadow fescue (Festuca pratensis Huds.). Theoret. Appl. Genet. 123, 369–382. doi: 10.1007/s00122-011-1590-z
Baldoni, E., Bagnaresi, P., Locatelli, F., Mattana, M., and Genga, A. (2016). Comparative Leaf and Root Transcriptomic Analysis of two Rice Japonica Cultivars Reveals Major Differences in the Root Early Response to Osmotic Stress. Rice 9:25. doi: 10.1186/s12284-016-0098-1
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
Bui Van, L., and Ba, N. (2019). Yield and Nutritional Value of Brachiaria humidicola Grass Planted in Thua Thien Hue Province, Vietnam. Int. J. Res. Agricult. Sci. 6, 40–45.
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: architecture and applications. BMC Bioinformatics 10:421. doi: 10.1186/1471-2105-10-421
Cardoso, J. A., Pineda, M., Jimenez, J. D., Vergara, M. F., and Rao, I. M. (2015). Contrasting strategies to cope with drought conditions by two tropical forage C-4 grasses. AoB Plants 7:lv107. doi: 10.1093/aobpla/plv107
Chaichi, M., Sanjarian, F., Razavi, K., and Gonzalez-Hernandez, J. L. (2019). Analysis of transcriptional responses in root tissue of bread wheat landrace (Triticum aestivum L.) reveals drought avoidance mechanisms under water scarcity. PLoS One 14:212671. doi: 10.1371/journal.pone.0212671
Cheruiyot, D., Midega, C. A. D., Van den Berg, J., Pickett, J. A., and Khan, Z. R. (2018). Genotypic responses of Brachiaria grass (Brachiaria spp.) accessions to drought stress. J. Agronomy 17, 136–146. doi: 10.3923/ja.2018.136.146
Conesa, A., and Gotz, S. (2008). Blast2GO: A comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008:619832. doi: 10.1155/2008/619832
Csiszár, J., Gallé, Á, Horváth, E., Dancsó, P., Gombos, M., Váry, Z., et al. (2012). Different peroxidase activities and expression of abiotic stress-related peroxidases in apical root segments of wheat genotypes with different drought stress tolerance under osmotic stress. Plant Physiol. Biochem. 52, 119–129. doi: 10.1016/j.plaphy.2011.12.006
Dinkins, R. D., Nagabhyru, P., Young, C. A., West, C. P., and Schardl, C. L. (2019). Transcriptome Analysis and Differential Expression in Tall Fescue Harboring Different Endophyte Strains in Response to Water Deficit. Plant Genome 12:71.
Dixon, D. P., Hawkins, T., Hussey, P. J., and Edwards, R. (2009). Enzyme activities and subcellular localization of members of the Arabidopsis glutathione transferase superfamily. J. Exp. Bot. 60, 1207–1218. doi: 10.1093/jxb/ern365
Finn, R. D., Coggill, P., Eberhardt, R. Y., Eddy, S. R., Mistry, J., Mitchell, A. L., et al. (2016). The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 44, D279–D285. doi: 10.1093/nar/gkv1344
Fisher, L. D., van Belle, G., Heagerty, P., and Lumley, T. (2004). Biostatistics: A Methodology for the Health Sciences. New Jersey: John Wiley and Sons.
Foito, A., Byrne, S. L., Shepherd, T., Stewart, D., and Barth, S. (2009). Transcriptional and metabolic profiles of Lolium perenne L. genotypes in response to a PEG-induced water stress. Plant Biotechnol. J. 7, 719–732. doi: 10.1111/j.1467-7652.2009.00437.x
Fradera-Sola, A., Thomas, A., Gasior, D., Harper, J., Hegarty, M., Armstead, I., et al. (2019). Differential gene expression and gene ontologies associated with increasing water-stress in leaf and root transcriptomes of perennial ryegrass (Lolium perenne). PLoS One 14:e0220518. doi: 10.1371/journal.pone.0220518
Guenni, O., Baruch, Z., and Marin, D. (2004). Responses to drought of five Brachiaria species. II. Water relations and leaf gas exchange. Plant Soil 258, 249–260. doi: 10.1023/b:plso.0000016555.53797.58
Guenni, O., Marin, D., and Baruch, Z. (2002). Responses to drought of five Brachiaria species. I. Biomass production, leaf growth, root distribution, water use and forage quality. Plant Soil 243, 229–241. doi: 10.1023/a:1019956719475
Haft, D. H., Selengut, J. D., Richter, R. A., Harkins, D., Basu, M. K., and Beck, E. (2013). TIGRFAMs and Genome Properties in 2013. Nucleic Acids Res. 41, D387–D395. doi: 10.1093/nar/gks1234
Hu, L., Xie, Y., Fan, S. J., Wang, Z. S., Wang, F. H., Zhang, B., et al. (2018). Comparative analysis of root transcriptome profiles between drought-tolerant and susceptible wheat genotypes in response to water stress. Plant Sci. 272, 276–293. doi: 10.1016/j.plantsci.2018.03.036
Jank, L., Barrios, S. C., do Valle, C. B., Simeao, R. M., and Alves, G. F. (2014). The value of improved pastures to Brazilian beef production. Crop Pasture Sci. 65, 1132–1137. doi: 10.1071/cp13319
Ji, Y., Chen, P. L., Chen, J., Pennerman, K. K., Liang, X. Y., Yan, H. D., et al. (2018). Combinations of Small RNA, RNA, and Degradome Sequencing Uncovers the Expression Pattern of microRNA-mRNA Pairs Adapting to Drought Stress in Leaf and Root of Dactylis glomerata L. Int. J. Mol. Sci. 19:19103114. doi: 10.3390/ijms19103114
Jiang, Y. W., Wang, X. C., Yu, X. Q., Zhao, X. W., Luo, N., Pei, Z. Y., et al. (2017). Quantitative Trait Loci Associated with Drought Tolerance in Brachypodium distachyon. Front. Plant Sci. 8:811. doi: 10.3389/fpls.2017.00811
Jones, P., Binns, D., Chang, H. Y., Fraser, M., Li, W., McAnulla, C., et al. (2014). InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240. doi: 10.1093/bioinformatics/btu031
Kanehisa, M., and Goto, S. (2000). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30.
Keller-Grein, G., Maass, B. L., and Hanson, J. (1996). “Natural variation in Brachiaria and existing germplasm collections,” in Brachiaria: biology, agronomy, and improvement, eds J. W. Miles, B. L. Maass, C. B. do Valle, and V. Kumble (Cali: International Center for Tropical Agriculture), 16–35.
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317
Kumar, J., Gunapati, S., Kianian, S. F., and Singh, S. P. (2018). Comparative analysis of transcriptome in two wheat genotypes with contrasting levels of drought tolerance. Protoplasma 255, 1487–1504. doi: 10.1007/s00709-018-1237-x
Kuwi, S. O., Kyalo, M., Mutai, C. K., Mwilawa, A., Hanson, J., Djikeng, A., et al. (2018). Genetic diversity and population structure of Urochloa grass accessions from Tanzania using simple sequence repeat (SSR) markers. Braz. J. Bot. 41, 699–709. doi: 10.1007/s40415-018-0482-8
Lees, J., Yeats, C., Redfern, O., Clegg, A., and Orengo, C. (2010). Gene3D: merging structure and function for a Thousand genomes. Nucleic Acids Res. 38, D296–D300. doi: 10.1093/nar/gkp987
Letunic, I., Copley, R. R., Pils, B., Pinkert, S., Schultz, J., and Bork, P. (2006). SMART 5: domains in the context of genomes and networks. Nucleic Acids Res. 34, D257–D260. doi: 10.1093/nar/gkj079
Liang, J. J., Chen, X., Deng, G. B., Pan, Z. F., Zhang, H. L., Li, Q., et al. (2017). Dehydration induced transcriptomic responses in two Tibetan hulless barley (Hordeum vulgare var. nudum) accessions distinguished by drought tolerance. BMC Genomics 18:775. doi: 10.1186/s12864-017-4152-1
Ling, H., Xie, Y., Fan, S. J., Wang, Z. S., Wang, F. H., Zhang, B., et al. (2018). Comparative analysis of root transcriptome profiles between drought-tolerant and susceptible wheat genotypes in response to water stress (vol 272, pg 276, 2018). Plant Sci. 276, 162–162. doi: 10.1016/j.plantsci.2018.08.010
Liu, S., and Jiang, Y. (2010). Identification of differentially expressed genes under drought stress in perennial ryegrass. Physiol. Plant 139, 375–387. doi: 10.1111/j.1399-3054.2010.01374.x
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Ma, Y. H., Yu, X. X., Yu, Z., Sun, F. C., Li, X. D., and Li, X. L. (2018). RNA-Seq of Agropyron mongolicum Keng in response to drought stress. Grassland Sci. 64, 3–15. doi: 10.1111/grs.12176
Maass, B. L., Midega, C. A. O., Mutimura, M., Rahetlah, V. B., Salgado, P., Kabirizi, J. M., et al. (2015). Homecoming of Brachiaria: Improved Hybrids Prove Useful for African Animal Agriculture. East Afr. Agricult. Forestry J. 81, 71–78. doi: 10.1080/00128325.2015.1041263
Meents, M. J., Watanabe, Y., and Samuels, A. L. (2018). The cell biology of secondary cell wall biosynthesis. Ann. Bot. 121, 1107–1125. doi: 10.1093/aob/mcy005
Merewitz, E., Belanger, F., Warnke, S., Huang, B. R., and Bonos, S. (2014). Quantitative Trait Loci Associated with Drought Tolerance in Creeping Bentgrass. Crop Sci. 54, 2314–2324. doi: 10.2135/cropsci2013.12.0810
Mi, H., Muruganujan, A., and Thomas, P. D. (2013). PANTHER in 2013: modeling the evolution of gene function, and other gene attributes, in the context of phylogenetic trees. Nucleic Acids Res. 41, D377–D386. doi: 10.1093/nar/gks1118
Miles, J. W. (2007). Apomixis for Cultivar Development in Tropical Forage Grasses. Crop Sci. 47, S–238–S–249. doi: 10.2135/cropsci2007.04.0016IPBS
Moin, M., Bakshi, A., Madhav, M. S., and Kirti, P. B. (2017). Expression Profiling of Ribosomal Protein Gene Family in Dehydration Stress Responses and Characterization of Transgenic Rice Plants Overexpressing RPL23A for Water-Use Efficiency and Tolerance to Drought and Salt Stresses. Front. Chem. 5:97–97. doi: 10.3389/fchem.2017.00097
Ougham, H. J., and Davies, T. G. E. (1990). LEAF DEVELOPMENT IN LOLIUM-TEMULENTUM - GRADIENTS OF RNA COMPLEMENT AND PLASTID AND NONPLASTID TRANSCRIPTS. Physiol. Plant. 79, 331–338.
Pan, L., Zhang, X., Wang, J., Ma, X., Zhou, M., Huang, L., et al. (2016). Transcriptional Profiles of Drought-Related Genes in Modulating Metabolic Processes and Antioxidant Defenses in Lolium multiflorum. Front. Plant Sci. 7:519. doi: 10.3389/fpls.2016.00519
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419. doi: 10.1038/nmeth.4197
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Pessino, S. C., Evans, C., Ortiz, J. P. A., Armstead, I., Do Valle, C. B., and Hayward, M. D. (1998). A genetic map of the apospory-region in Brachiaria hybrids: identification of two markers closely associated with the trait. Hereditas 128, 153–158. doi: 10.1111/j.1601-5223.1998.00153.x
Pessino, S. C., Ortiz, J. P. A., Leblanc, O., doValle, C. B., Evans, C., and Hayward, M. D. (1997). Identification of a maize linkage group related to apomixis in Brachiaria. Theoret. Appl. Genet. 94, 439–444. doi: 10.1007/s001220050434
Pessoa-Filho, M., Azevedo, A. L. S., Sobrinho, F. S., Gouvea, E. G., Martins, A. M., and Ferreira, M. E. (2015). Genetic Diversity and Structure of Ruzigrass Germplasm Collected in Africa and Brazil. Crop Sci. 55, 2736–2745. doi: 10.2135/cropsci2015.02.0096
Pessoa-Filho, M., Martins, A. M., and Ferreira, M. E. (2017). Molecular dating of phylogenetic divergence between Urochloa species based on complete chloroplast genomes. BMC Genomics 18:516. doi: 10.1186/s12864-017-3904-2
Petter, F. A., Pacheco, L. P., Zuffo, A. M., Piauilino, A. C., Xavier, Z. F., dos Santos, J. M., et al. (2013). Performance of cover crop under the water deficit. Semina-Ciencias Agrarias 34, 3307–3319. doi: 10.5433/1679-0359.2013v34n6Supl1p3307
Pizarro, E. A., Hare, M. D., Mutimura, M., and Changun, B. (2013). Brachiaria hybrids: potential, forage use and seed yield. Tropical Grasslands – Forrajes Tropicales 1, 31–35.
Podgorska, A., Burian, M., and Szal, B. (2017). Extra-Cellular But Extra-Ordinarily Important for Cells: Apoplastic Reactive Oxygen Species Metabolism. Frontiers in Plant Science 8:353. doi: 10.3389/fpls.2017.01353
Ranjan, A., Pandey, N., Lakhwani, D., Dubey, N. K., Pathre, U. V., and Sawant, S. V. (2012). Comparative transcriptomic analysis of roots of contrasting Gossypium herbaceum genotypes revealing adaptation to drought. BMC Genomics 13:680. doi: 10.1186/1471-2164-13-680
Renvoize, S. A., Clayton, W. D., and Kabuye, C. H. S. (1996). “Morphology, taxonomy, and natural distribution of Brachiaria (Trin.) Griseb,” in Brachiaria: biology, agronomy, and improvement, eds J. W. Miles, B. L. Maass, C. B. do Valle, and V. Kumble (Cali: International Center for Tropical Agriculture), 1–15.
Salgado, L. R., Lima, R., dos Santos, B. F., Shirakawa, K. T., Vilela, M. D., Almeida, N. F., et al. (2017). De novo RNA sequencing and analysis of the transcriptome of signalgrass (Urochloa decumbens) roots exposed to aluminum. Plant Growth Regulation 83, 157–170. doi: 10.1007/s10725-017-0291-2
Santos, P. M., da Cruz, P. G., de Araujo, L. C., Pezzopane, J. R. M., do Valle, C. B., and Pezzopane, C. D. (2013). Response mechanisms of Brachiaria brizantha cultivars to water deficit stress. Revista Brasileira De Zootecnia-Brazilian Journal of Animal Science 42, 767–773. doi: 10.1590/s1516-35982013001100001
Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V., and Zdobnov, E. M. (2015). BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212. doi: 10.1093/bioinformatics/btv351
Tegegn, A., Kyalo, M., Mutai, C., Hanson, H., Asefa, G., Djikeng, A., et al. (2019). Genetic diversity and population structure of Brachiaria brizantha (A.Rich.) Stapf accessions from Ethiopia. African Journal of Range and Forest Science 36, 129–133. doi: 10.2989/10220119.2019.1573760
Tenhaken, R. (2015). Cell wall remodeling under abiotic stress. Frontiers in Plant Science 5:771. doi: 10.3389/fpls.2014.00771
Thaikua, S., Ebina, M., Yamanaka, N., Shimoda, K., Suenaga, K., and Kawamoto, Y. (2016). Tightly clustered markers linked to an apospory-related gene region and quantitative trait loci mapping for agronomic traits in Brachiaria hybrids. Grassland Science 62, 69–80. doi: 10.1111/grs.12115
Tohme, J., Palacios, N., Lenis, S., and Roca, W. (1996). “Applications of biotechnology to Brachiaria,” in Brachiaria: biology, agronomy, and improvement, eds J. W. Miles, B. L. Maass, C. B. do Valle, and V. Kumble (Cali: International Center for Tropical Agriculture).
Triviño, N. J., Perez, J. G., Recio, M. E., Ebina, M., Yamanaka, N., Tsuruta, S.-I., et al. (2017). Genetic Diversity and Population Structure of Brachiaria Species and Breeding Populations. Crop Science 57, 2633–2644. doi: 10.2135/cropsci2017.01.0045
Vigna, B. B. Z., de Oliveira, F. A., de Toledo-Silva, G., da Silva, C. C., do Valle, C. B., and de Souza, A. P. (2016). Leaf transcriptome of two highly divergent genotypes of Urochloa humidicola (Poaceae), a tropical polyploid forage grass adapted to acidic soils and temporary flooding areas. BMC Genomics 17:5. doi: 10.1186/s12864-016-3270-5
Wenzl, P., Arango, A., Chaves, A. L., Buitrago, M. E., Patiño, G. M., Miles, J., et al. (2006). A Greenhouse Method to Screen Brachiariagrass Genotypes for Aluminum Resistance and Root Vigor This work was performed at CIAT and funded by the Colombian Ministry of Agriculture and Rural Development, the Commission of Developmental Issues of the Austrian Academy of Sciences (Project No. 78), and the German Ministry of Economic Cooperation and Development (BMZ/GTZ 2000.7860.0-001.0). Crop Science 46, 968–973. doi: 10.2135/cropsci2005.07.0209
Wenzl, P., Mancilla, L., Mayer, J., Albert, R., and Rao, I. (2003). Simulating Infertile Acid Soils with Nutrient Solutions: The Effects on Brachiaria Species. Soil Science Society of America Journal 67, 1457. doi: 10.2136/sssaj2003.1457
Wingett, S. W., and Andrews, S. (2018). FastQ Screen: A tool for multi-genome mapping and quality control. F1000Res 7, 1338. doi: 10.12688/f1000research.15931.2
Worthington, M., Ebina, M., Yamanaka, N., Heffelfinger, C., Quintero, C., Zapata, Y. P., et al. (2019). Translocation of a parthenogenesis gene candidate to an alternate carrier chromosome in apomictic Brachiaria humidicola. BMC Genomics 20:41. doi: 10.1186/s12864-018-5392-4
Worthington, M., Heffelfinger, C., Bernal, D., Quintero, C., Zapata, Y. P., Perez, J. G., et al. (2016). A Parthenogenesis Gene Candidate and Evidence for Segmental Allopolyploidy in Apomictic Brachiaria decumbens. Genetics 203, 1117. doi: 10.1534/genetics.116.190314
Worthington, M., Perez, J. G., Mussurova, S., Silva-Cordoba, A., Castiblanco, V., Cardoso Arango, J. A., et al. (2020). A new genome allows the identification of genes associated with natural variation in aluminium tolerance in Brachiaria grasses. J Exp Bot 2020, eraa469. doi: 10.1093/jxb/eraa469
Zhao, P. C., Liu, P. P., Yuan, G. X., Jia, J. T., Li, X. X., Qi, D. M., et al. (2016). New Insights on Drought Stress Response by Global Investigation of Gene Expression Changes in Sheepgrass (Leymus chinensis). Frontiers in Plant Science 7:954. doi: 10.3389/fpls.2016.00954
Zhu, Y. Q., Wang, X., Huang, L. K., Lin, C. W., Zhang, X. Q., Xu, W. Z., et al. (2017). Transcriptomic Identification of Drought-Related Genes and SSR Markers in Sudan Grass Based on RNA-Seq. Frontiers in Plant Science 8:687. doi: 10.3389/fpls.2017.00687
Keywords: Brachiaria, drought, differentially expressed genes (DEGs), comparative transcriptomics, functional enrichment
Citation: Jones C, De Vega J, Worthington M, Thomas A, Gasior D, Harper J, Doonan J, Fu Y, Bosch M, Corke F, Arango J, Cardoso JA, de la Cruz Jimenez J, Armstead I and Fernandez-Fuentes N (2021) A Comparison of Differential Gene Expression in Response to the Onset of Water Stress Between Three Hybrid Brachiaria Genotypes. Front. Plant Sci. 12:637956. doi: 10.3389/fpls.2021.637956
Received: 04 December 2020; Accepted: 19 February 2021;
Published: 19 March 2021.
Edited by:
Raul Antonio Sperotto, Universidade do Vale do Taquari – Univates, BrazilReviewed by:
Yucheng Wang, Shenyang Agricultural University, ChinaMehanathan Muthamilarasan, University of Hyderabad, India
Linhai Wang, Oil Crops Research Institute, Chinese Academy of Agricultural Sciences, China
Copyright © 2021 Jones, De Vega, Worthington, Thomas, Gasior, Harper, Doonan, Fu, Bosch, Corke, Arango, Cardoso, de la Cruz Jimenez, Armstead and Fernandez-Fuentes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Narcis Fernandez-Fuentes, naf4@aber.ac.uk; Ian Armstead, ipa@aber.ac.uk