Takeshi Kuroha1,2*
Keisuke Nagai1
Yusuke Kurokawa1
Yoshiaki Nagamura3
Miyako Kusano4,5
Hideshi Yasui6Motoyuki Ashikari1
Atsushi Fukushima4*- 1Bioscience and Biotechnology Center, Nagoya University, Nagoya, Japan
- 2Department of Developmental Biology and Neurosciences, Graduate School of Life Sciences, Tohoku University, Sendai, Japan
- 3Genome Resource Unit, National Institute of Agrobiological Sciences, Tsukuba, Japan
- 4RIKEN Center for Sustainable Resource Science, Yokohama, Japan
- 5Graduate School of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Japan
- 6Faculty of Agriculture, Kyushu University, Fukuoka, Japan
To avoid low oxygen, oxygen deficiency or oxygen deprivation, deepwater rice cultivated in flood planes can develop elongated internodes in response to submergence. Knowledge of the gene regulatory networks underlying rapid internode elongation is important for an understanding of the evolution and adaptation of major crops in response to flooding. To elucidate the genetic and molecular basis controlling their deepwater response we used microarrays and performed expression quantitative trait loci (eQTL) and phenotypic QTL (phQTL) analyses of internode samples of 85 recombinant inbred line (RIL) populations of non-deepwater (Taichung 65)- and deepwater rice (Bhadua). After evaluating the phenotypic response of the RILs exposed to submergence, confirming the genotypes of the populations, and generating 188 genetic markers, we identified 10,047 significant eQTLs comprised of 2,902 cis-eQTLs and 7,145 trans-eQTLs and three significant eQTL hotspots on chromosomes 1, 4, and 12 that affect the expression of many genes. The hotspots on chromosomes 1 and 4 located at different position from phQTLs detected in this study and other previous studies. We then regarded the eQTL hotspots as key regulatory points to infer causal regulatory networks of deepwater response including rapid internode elongation. Our results suggest that the downstream regulation of the eQTL hotspots on chromosomes 1 and 4 is independent, and that the target genes are partially regulated by SNORKEL1 and SNORKEL2 genes (SK1/2), key ethylene response factors. Subsequent bioinformatic analyses, including gene ontology-based annotation and functional enrichment analysis and promoter enrichment analysis, contribute to enhance our understanding of SK1/2-dependent and independent pathways. One remarkable observation is that the functional categories related to photosynthesis and light signaling are significantly over-represented in the candidate target genes of SK1/2. The combined results of these investigations together with genetical genomics approaches using structured populations with a deepwater response are also discussed in the context of current molecular models concerning the rapid internode elongation in deepwater rice. This study provides new insights into the underlying genetic architecture of gene expression regulating the response to flooding in deepwater rice and will be an important community resource for analyses on the genetic basis of deepwater responses.
Introduction
Flooding, an environmental stress, affects the growth and development of most plants by limiting the exchange of gases such as oxygen and carbonic dioxide and reducing the light intensity; this reduces yield around the world. Some plants such as Oryza, Rumex, Rorippa, and Echinochloa genera can survive under this stressful condition by implementing specific strategies (Bailey-Serres and Voesenek, 2008). Flood-tolerant Oryza sativa (rice) developed two different strategies for surviving full and partial submergence (Voesenek and Bailey-Serres, 2015; Loreti et al., 2016). The ‘quiescent strategy’ allows flash-flood-tolerant rice to resume growth after the floodwater recedes (Bailey-Serres and Voesenek, 2010; Mickelbart et al., 2015). By the ‘escape strategy,’ the internodes of rice plants tolerant to deepwater flooding rapidly elongate as the water level rises (Nagai et al., 2010; Hattori et al., 2011). Consequently, deepwater rice can be grown in many flood-prone areas in South and Southeast Asia such as Bangladesh, Cambodia, and Thailand (Kende et al., 1998).
Many agronomically important traits, including tolerance to flash- and deepwater floods, tend to be governed by genes known as quantitative trait loci (QTL) (Quarrie et al., 1997; Price et al., 2002; Ashikari and Matsuoka, 2006; Salvi and Tuberosa, 2015). An Indian rice variety, FR13A that applies the quiescent strategy for submergence tolerance, harbors the major QTL Submergence 1 (SUB1) (Nandi et al., 1997; Siangliw et al., 2003). SUB1A, a member of the APETALA2 (AP2)/Ethylene Response Factor (ERF) superfamily, has been cloned as the gene responsible for SUB1 QTL (Xu et al., 2006). QTL analysis has also shown that the internode elongation trait in deepwater rice is regulated by three major QTLs (Tang et al., 2005; Hattori et al., 2007, 2008; Kawano et al., 2008). Subsequent work demonstrated that the QTL in chromosome 12, the major QTL controlling the deepwater response, harbors two transcription factors, SNORKEL1 and SNORKEL2 (SK1/2) (Hattori et al., 2009), members of the AP2/ERF superfamily. These findings suggest that deepwater rice acquired key regulators to adapt to different types of flood. However, our understanding of the gene regulatory networks underlying these traits in flood-tolerant rice in response to submergence remains far from complete and no other key regulators for submergence tolerance in deepwater rice have been identified.
In flash flood-tolerant rice, Rumex, Rorippa, and Echinochloa, large-scale transcriptome analysis using microarrays or RNA sequencing (RNA-seq) has been performed focusing on mechanism of flooding tolerance (Jung et al., 2010; Sasidharan et al., 2013; Van Veen et al., 2013; Nah et al., 2015). Progress in the acquisition of large-scale transcriptome data has greatly facilitated expression QTL (eQTL) analysis (Jansen and Nap, 2001; Kliebenstein, 2009; Druka et al., 2010; Cubillos et al., 2012a; Chitwood and Sinha, 2013); the gene expression level is considered a trait (expression trait or e-trait). eQTL studies of the model plant Arabidopsis thaliana (Decook et al., 2006; Kliebenstein et al., 2006; Keurentjes et al., 2007; West et al., 2007; Burow et al., 2009; Jimenez-Gomez et al., 2010; Cubillos et al., 2012b; Lowry et al., 2013), of crops (Potokina et al., 2008; Swanson-Wagner et al., 2009; Chen et al., 2010; Holloway et al., 2011; Bolon et al., 2014; Wang Y. et al., 2014; Ranjan et al., 2016), and of trees (Drost et al., 2010, 2015) mapped expression variations to both cis- and trans-acting regulatory polymorphisms: cis-eQTL (DNA polymorphisms within or close to the gene) and trans-eQTL (polymorphisms located anywhere in the genome). In rice, Wang et al. (2010) performed eQTL mapping of recombinant inbred lines (RILs) derived from two indica varieties, Zhenshan 97 and Minghui 63, to identify gene regulatory networks that contribute to an agriculturally important trait. They subsequently developed a new framework that integrates eQTL- and co-expression analysis to narrow the pool of candidate causal genes (Wang J. et al., 2014). An approach to integrate genetic and genomic data, called genetical genomics (Jansen and Nap, 2001), would be useful for identifying possible candidate genes and molecular regulatory mechanisms that play a role in the rapid internode elongation seen in deepwater rice as the water level rises. Gene expression profiling using structured populations segregated into RILs and derived from non-deepwater- and deepwater rice may yield additional insights into the etiology of deepwater responses.
To elucidate the genetic and molecular basis controlling this important trait, we performed microarray experiments in our eQTL analysis of the shoot samples from RIL populations of non-deepwater [Taichung 65 (T65)]- and deepwater rice (Bhadua) subjected to submergence. We identified 10,047 significant eQTLs and significant trans-eQTL hotspots on chromosomes 1, 4, and 12 that regulate variations in transcript abundance and exert strong effects on gene expression levels. We also evaluated the QTL for phenotypic traits (phQTL) of the same line. We discuss candidate genes located in the novel eQTL hotspots and candidate target genes of SK1/2. This study provides a useful resource for reconstructing the transcriptional regulatory network at the system-level, and for cloning genes and regulatory control mechanisms that underlie deepwater responses.
Materials and Methods
Plant Materials and Growth Conditions
Rice RILs from a cross between Taichung 65 (O. sativa L. ssp. japonica), and Bhadua (O. sativa L. ssp. indica; Bangladesh) were used in this study. Seed stocks for the RILs were provided by the National Bio-Resource Project (NBRP) in Japan. We used 85 RILs described by Nagai et al. (2014). The rice seeds were sterilized by boiling at 60°C for 10 min, allowed to germinate in water (30°C for 72 h), and then sown in June 2011 in perforated plastic pots (diameter 10 cm, height 12 cm) filled with soil. They were grown in the greenhouse under a natural light environment. During growth, the water level was beneath the soil surface (air conditions). One month after germination, eight-leaf-stage plants were put into the tank of 160 cm height and completely submerged in water for phenotypic evaluation or microarray analysis.
Phenotypic Evaluation of Total Internode Length, Plant Height, and Number of Elongated Internodes
Quantitative evaluation of the deepwater response of all RILs and their parental lines was as described elsewhere (Hattori et al., 2009; Nagai et al., 2014). Seven days after submergence treatment, we measured the total internode length (TIL), defined as the length between the uppermost and the basal node, the plant height (PH), and the number of elongated internodes (NEI) in the RILs (Hattori et al., 2009; Nagai et al., 2014) (Supplementary Figure S1A). The main culms of each plant were cut along the midline and the TIL was recorded. The NEI was defined as the number of internodes longer than 5 mm. The phenotype of each line was based on the average of no fewer than three plants.
Genotyping with Illumina’s BeadArray/BeadChip Technology
Genomic DNA was extracted from 85 individuals of the T65/Bhadua RILs and their parental lines using the isoplant method (Zhu et al., 1993) to create a linkage map (Nagai et al., 2014). Briefly, the extracted DNA samples were subjected to genotyping with single-nucleotide polymorphism (SNP) markers using an Affymetrix customized SNP array for rice (Kurokawa et al., 2016). SNPs were detected with the Golden Gate assay of BeadXpress (Illumina, Inc., San Diego, CA, United States). To build the linkage map we used MAPMAKER/EXP version 3.0 (Lander et al., 1987). Our RIL population was genotyped for 188 markers as shown in Supplementary Figure S2.
Microarray Analysis
We used the Agilent Rice Oligo Microarray (44k, custom-made; Agilent Technologies, Redwood City, CA, United States) with the one-color method (Yano et al., 2012). The design of the microarray platform was based on the curated annotation of the rice genome (Sakai et al., 2013) and contained, in addition to four probe-sets corresponding to SK1 and SK2 genes, an additional 56 probe-sets corresponding to hormone-related genes and micro RNAs. Since key genes for submergence response are drastically induced until 6 h after submergence treatment (Fukao et al., 2006; Xu et al., 2006; Hattori et al., 2009), we performed the microarray analysis with samples including elongating internodes extracted 6 h after submergence (Supplementary Figure S1B). After submergence treatment, the leaves were removed and total RNA was isolated from the region between 1 cm upper part from the uppermost and the basal node in the main culms using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, United States). Fluorescent probe labeling was with the Agilent Quick Amp Labeling Kit. Labeled cRNA was then fragmented and hybridized on 4 × 44k microarray slides using the Agilent Gene Expression Hybridization Kit in a hybridization oven (65°C for 17 h). All microarray slides were then scanned on an Agilent DNA microarray scanner (model G2505B with G2565BA). Scanned data were analyzed using the default settings of Feature Extraction software (Agilent Technologies). We used unreplicated data for each RIL. Annotation information was from the Rice Annotation Project Database (RAP-DB) (Sakai et al., 2013). All microarray data are in the Gene Expression Omnibus (GEO) (Barrett et al., 2013) and accessible via GEO series accession number GSE87702. Analysis of the microarray data was with R1 and the Bioconductor (Gentleman et al., 2004). The Agi4X44PreProcess package2 was used for pre-processing and normalization of the microarray data. For quality control and visualization we compared the box plots for raw intensity and background-corrected processed signals; this yielded proper normalized data (Supplementary Figure S3A). Using the annotation of probe-sets and rice gene annotation, multiple probe-set identifiers associated with a single gene were summarized to obtain a gene identifier (RAP Id) for the determination of the number of unique eQTLs. The expression value of a gene was averaged across multiple probe-sets mapped to the same gene.
eQTL and phQTL Mapping
For eQTL mapping, we performed standard quality control, pre-processing of multiple probe-sets assigned to a single gene, and variance filtering. We used 12,264 e-traits in further eQTL mapping to increase the statistical power (Bourgon et al., 2010). Normalized microarray data were used for eQTL mapping with the R package “eqtl” (Cubillos et al., 2012b). The package provides, for example, functions to automatically detect QTLs with supporting intervals from the results of interval mapping for the classification of cis- and trans acting eQTLs, and for the visualization of the genome-wide relationship between e-traits and eQTLs. Interval mapping was calculated by Haley-Knott regression (Haley and Knott, 1992) implemented in the R/qtl package (Arends et al., 2010). We classified the acting type of eQTLs using the R/eqtl::classify.qtl() function. The function estimates whether an eQTL is cis- or trans-acting entirely based on logarithm of odds (LOD) support interval definition. The trans-eQTLs were defined as those which do not co-localize the affected gene. The cis-acting eQTL was defined as those which contains the regulated gene within their LOD support interval. The support interval (si) was estimated by the accepted drop of LOD score from the LOD peak. In this study, we used si = 1.5, default settings in R/eqtl::define.peak() function to define QTLs with support interval.
For phQTL mapping, we used the R package “qtl” (Arends et al., 2010). LOD thresholds were obtained using 1,000 permutations for a significance level of 5%. A global permutation threshold (GPT) was used for determining a genome-wide threshold for statistically significant eQTLs. The significance level for the GPT was set at 0.05. To identify eQTL hotspots, clusters of e-traits co-mapping to the same genomic position, we used the R package “qtlhot” (Neto et al., 2012), which is based on quantile-based permutations.
Statistical Data Analysis
To identify differentially expressed genes (DEGs) between deepwater and non-deepwater phenotypes, we calculated log2-mean expression fold-change of each gene in deepwater type (TIL ≥ 20 cm; Bhadua and 44 RILs) compared to that of non-deepwater type (TIL < 20 cm; T65 and 41 RILs). For the calculation, we used the R package LIMMA (Smyth, 2005); it features false discovery rate (FDR) correction for multiple testing (Benjamini and Hochberg, 1995). Hierarchical cluster analyses (HCAs) and visualization of the expression level (log2 ratio) were performed in R. Euclidean distances and the complete linkage method were used for HCA. Principal component analysis (PCA) was applied for the expression data matrix (log10-transformed) with autoscaling using SIMCA-P + (v 14.0, Umetrics, Sweden) software. Gene Ontology (GO) term enrichment analysis was with BiNGO software (Maere et al., 2005) with Benjamini and Hochberg (1995) correction for multiple testing (FDR < 0.05).
Promoter Enrichment Analysis
Our promoter enrichment analysis was based on R-scripts called MotifEnrichment3 (Ichihashi et al., 2014; Chitwood et al., 2015). We obtained the 1,000 basepairs upstream of the ATG translational start site of genes with significant trans-eQTL hotspots on chromosomes 1, 4, and 12 using the Bioconductor package “BiomaRt” (Durinck et al., 2005, 2009). We employed 100 curated motifs (filename: element_name_and_motif_IUPACsupp.txt) that are based on the Gene Regulatory Information Server AtTFDB database4 (Yilmaz et al., 2011). The abundance of motif elements in promoters of genes with trans-eQTLs was compared with that of all considered genes using Fisher’s exact test and the Biostrings package (Pages et al., 2016).
Visualization of phQTLs and Transcript Profiling in a Chromosomal Context
We employed the customized version of QTLVizR5 for the visual inspection of phQTLs and for transcript profiling in a chromosomal context. The tool is a web-based application software originally developed by Maloof Lab at the University of California at Davis. The implementation is based on R/shiny6; it renders visualization highly dynamic. All data and R/shiny scripts are available upon request to the authors.
Results
Phenotypic Evaluation and Transcriptome Analysis
We performed phenotype analysis of the internode elongation under submergence and genome-wide transcript profiling of the internodes in the parents and 85 RILs, crosses between T65 (non-deepwater rice) and Bhadua (deepwater rice). First, we carried out phenotypic evaluation of TIL, PH, and NEI in the RILs (Supplementary Figure S1A). Based on 188 genetic markers (Supplementary Figure S2), the haplotypes and the physical genome position of the marker recombination bins for all RILs were investigated (Figure 1A). Then, we evaluated phQTLs using quantitative phenotype and haplotype data from the RILs (Supplementary Figures S4, S5). Composite interval mapping detected only two significant phQTLs for TIL and PH (LOD = 6.0 at pseudo-marker position “c1.loc40” on chromosome 1 and 11.4 at marker position “ad12009568” on chromosome 12) (Supplementary Figures S5A,B). Two significant phQTLs for NEI were located at marker position “P0770_1” on chromosome 1; on chromosome 12 they were at the same position for TIL and PH (Supplementary Figure S5C). The location of these two phQTLs for all phenotypic traits was consistent with the QTLs with respect to the internode elongation previously detected by mapping with T65/Bhadua (Kawano et al., 2008) and T65/C9285 (Hattori et al., 2007, 2008).
FIGURE 1. Overview of our genetic and transcriptomic data. (A) Haplotypes of 85 recombinant inbred lines (RILs) and their parental lines, T65, and Bhadua. Each row indicates an individual plant, columns correspond to one of the bins (188 markers) arranged in physical order based on the rice chromosomes [Os-Nipponbare-Reference-IRGSP-1.0 (Sakai et al., 2013)]. Red and blue identify the T65 and the Bhadua genotype, respectively. The white box represents missing genotypes. (B) Principal component analysis (PCA) of the expression data matrix of the population. The score scatter plot reveals a clear separation between individuals; the first principal component (PC1) accounts for 13.3% of the total explained variance. Orange stars and sky blue diamonds represent RILs with deepwater phenotype (TIL ≥ 20 cm) and non-deepwater phenotype (TIL < 20 cm), respectively. Brown and green rounded squares indicate T65 and Bhadua, respectively. The ellipse indicates Hotelling’s T2, with 95% confidence in the score scatter plots.
To assess the biological consistency of different microarray datasets from the same populations we applied PCA of the expression data matrix after appropriate pre-processing and normalization (Figure 1B and Supplementary Figure S3). The score scatter plot from PCA revealed a clear separation between T65- and Bhadua samples; the first principal component (PC1) accounted for 13.3% of the total explained variance. The plot of RILs indicates that the coordinates in non-deepwater type (TIL < 20 cm) and deepwater type (TIL ≥ 20 cm) are located opposite (Figure 1B). These results suggest that global expression pattern of deepwater rice is different from that of non-deepwater rice.
Identification and Genome-Wide Distribution of cis- and trans-eQTLs
Next, we extracted 12,264 e-traits from genome-wide transcript data by pre-processing and normalization. For eQTL mapping we applied interval mapping [Haley-Knott regression (Haley and Knott, 1992)] with R/qtl (Arends et al., 2010) and R/eqtl packages (Cubillos et al., 2012b). Subsequent mapping identified 10,047 significant eQTLs (FDR < 0.05) comprised of 2,902 cis-eQTLs (28.9%) and 7,145 trans-eQTLs (71.1%); 7,078 e-traits featured at least one eQTL (Supplementary Table S1). To determine the number of eQTLs for each e-trait we examined the genome-wide relationship between e-traits and eQTLs on a heatmap. This examination was based on the physical positions of the markers and their corresponding e-traits on the 12 rice chromosomes (Figure 2A). We also estimated the phenotypic contribution with the R-square based on analysis of variance (ANOVA) for individual eQTLs and their significant interactions for each e-trait (Figure 2B). We defined eQTLs based on a threshold of LOD significance calculated by permuting the phenotypes while maintaining the genotype across the RILs and a 1.5 LOD drop-support interval. On the 12 rice chromosomes the proportion of cis-eQTLs (Figure 2C) was homogeneously while the proportion of trans-eQTLs (Figure 2D) was heterogeneously distributed, especially on chromosomes 1 (18% of trans-eQTLs), 4 (14%), 9 (8%), and 12 (32%).
FIGURE 2. Summary of expression quantitative trait loci (eQTLs) and their location in the rice genome. We analyzed shoot tissue of germinating RIL seeds obtained by crossing T65 and Bhadua rice. (A) The eQTL heatmap for T65/Bhadua populations was significant at a 5% false discovery rate (FDR). The x-axis represents the genetic position of the detected eQTLs and the y-axis the physical position of expressed genes (e-traits). The 12 rice chromosomes are separated by dashed lines. The pseudo-colors indicate the LOD scores (from green to red in order of the increasing LOD score). The bar length along the x-axis indicates the eQTL support interval. (B) Distribution of the R-square based on analysis of variance (ANOVA), the proportion of phenotypic variation explained by the segregation of an individual QTL or the significant QTL interactions for all eQTLs. (C) Pie chart of the proportion of cis-eQTLs distributed on 12 rice chromosomes. (D) Pie chart of the trans-eQTLs distributed on 12 rice chromosomes.
Detection of cis-eQTLs for SK1/2
As described above, we detected phQTLs on end of chromosome 12 (Supplementary Figure S5), whose responsible gene would be the SK1/2 (Hattori et al., 2009). We next visualized LOD statistics of the phQTLs for TIL (Figure 3A, upper panel) and the DEGs, calculated from log2-fold-change between expressions of RILs with deepwater phenotype (TIL ≥ 20 cm) and non-deepwater phenotype (TIL < 20 cm) (Figure 3A, lower panel) along the physical genomic position on chromosome 12. The location of SK1/2 genes was consistent with that of phQTLs for TIL (Figure 3A). The peak of eQTLs for the expression of SK1/2 was detected as cis-eQTL in the region where their own genes were located (Figure 3B). Thus, the responsible genes of the cis-eQTL are at least SK1/2 themselves because these genes exist in Bhadua but not in T65 (Hattori et al., 2009).
FIGURE 3. Phenotypic QTL (phQTL) and eQTL for SK1/2. (A) LOD statistics of phQTLs for total internode length (TIL) on chromosome 12 are shown (upper panel). The differential expressions between non-deepwater- and deepwater-type RILs are plotted along the physical genome position in log2-fold change (bottom panel). Note that transcription factors SK1 and SK2 exist in Bhadua but not in T65 plants. The red horizontal line indicates LOD = 3. (B) LOD statistics of eQTLs for SK1/2 gene expressions. Vertical blue lines represent support interval for the detected QTLs.
The Identification of Significant eQTL Hotspots Revealed Novel Genomic Loci That May Be Responsible For Expression Variations between Non-deepwater- and Deepwater Rice
Groups of e-traits that co-map to the same genomic region are called eQTL hotspots (Brem et al., 2002; Schadt et al., 2003; Breitling et al., 2008). They are biologically meaningful because they may harbor key regulators. We identified the major significant eQTL hotspots using a quantile-based permutation approach called the NL-method (Neto et al., 2012); all permutation thresholds were calculated targeting a genome-wide error rate (GWER) < 0.05. This analysis revealed the major eQTL hotspots on chromosomes 1, 4, and 12 (Figure 4A and Supplementary Figure S6); some of their positions on these chromosomes are enriched for trans-eQTLs. The eQTL hotspot at pseudo-marker position “c12.loc25.5” on chromosome 12 was consistent with the location of SK1/2. The eQTL hotspots at “P0298” on chromosomes 1 and “c4.loc24” on chromosome 4 located at different position from phQTLs detected in this study (Supplementary Figure S5) and the three major phQTLs reported by others (Tang et al., 2005; Hattori et al., 2007, 2008; Kawano et al., 2008).
FIGURE 4. Expression quantitative trait loci hotspot size significance profile and functional assessment of genes with trans-eQTLs. (A) The major significant eQTL hotspots were derived with a quantile-based permutation approach (Neto et al., 2012). All permutation thresholds were calculated targeting a genome-wide error rate (GWER) < 0.05. The major hotspots on chromosomes 1, 4, and 12 were identified. Arrows indicate the (pseudo-)marker name of the maximum LOD peak. Black lines represent raw eQTL counts, red identifies counts smoothed over a 5 cM window, and blue the sliding LOD threshold approach. The horizontal dashed green line indicates 5% threshold of 195 that is for the number of e-traits with significant LOD at a position (Neto et al., 2012). (B) Venn diagram of the number of e-traits with trans-eQTL at the physical position of the peaks within the three eQTL hotspots. (C) Over-represented gene ontology (GO) terms among 638 genes with trans-eQTLs at physical pseudo-position “c12.loc25.5” on chromosome 12. We used BiNGO (Maere et al., 2005) (FDR < 0.05). The top three significantly enriched terms were “photosynthesis” (FDR = 8.8E-33), “photosynthesis, light reaction” (FDR = 1.3E-28), and “monosaccharide metabolic process” (FDR = 1.1E-20). (D) Over-represented GO terms among 55 genes with trans-eQTLs at physical position “P0298” on chromosome 1. The term “protein tetramerization” (GO-Id: 51262) was significantly enriched among the genes (FDR = 2.9E-2). Note that there are no significant GO terms in the analysis for genes with trans-eQTLs at physical pseudo-position “c4.loc24” on chromosome 4.
Candidate Target Genes of the Novel eQTL Hotspots
The hypothesis that the candidate target genes of the eQTL hotspots harbor trans-eQTLs at the physical position of each eQTL hotspot provides a clue for a potential regulatory network where the eQTL hotspots control changes in the transcript abundance of numerous genes in trans. This may result in phenotypic variations. First, we narrowed the pool down to 55, 56, and 638 e-traits that feature trans-eQTL at physical position “P0298” on chromosome 1, “c4.loc24” on chromosome 4, and “c12.loc25.5” on chromosome 12, respectively (Supplementary Table S2). Next, from these e-traits, we searched for DEGs, between expressions of deepwater phenotype and non-deepwater phenotype on physical position of RILs for each e-traits. We found that 23.6% (13/55 genes) of e-traits with trans-eQTL at “P0298,” 10.7% (6/56 genes) of e-traits with trans-eQTL at “c4.loc24,” and 30.1% (192/638 genes) of e-traits with trans-eQTL at “c12.loc25.5” appeared in groups of significantly up-regulated DEGs (Supplementary Table S2).
Analysis of the relationships among the three groups of e-traits at different physical positions revealed no e-trait sharing trans-eQTLs at “P0298” and “c4.loc24” (Figure 4B). On the other hand, 45.5% (25/55 genes) of e-traits with trans-eQTL at “P0298” and 10.7% (6/56 genes) of e-traits with trans-eQTL at “c4.loc24” are shared with trans-eQTLs at position “c12.loc25.5.” To assess whether the eQTL hotspots are associated with biological functions, we conducted GO (The Gene Ontology Consortium, 2014) term enrichment analysis of genes with trans-eQTL at physical (pseudo-)position “c12.loc25.5” on chromosome 12 and at “P0298” on chromosome 1 using BiNGO (Maere et al., 2005) (FDR < 0.05). The three top significantly enriched terms among genes with trans-eQTLs at physical pseudo-position “c12.loc25.5” on chromosome 12 were “photosynthesis” (FDR = 8.8E-33), “photosynthesis, light reaction” (FDR = 1.3E-28), and “monosaccharide metabolic process” (FDR = 1.1E-20) (Figure 4C and Supplementary Table S3B), indicating potential SK1/2 target genes. The term “protein tetramerization” (GO-Id: 51262) was significantly enriched among genes with trans-eQTLs at physical positions on chromosome 1 (FDR = 2.9E-2) (Figure 4D and Supplementary Table S3A). As for physical pseudo-position “c4.loc24” on chromosome 4, our analysis for genes with trans-eQTLs revealed no significant GO terms.
Candidate Causal Genes Located in the Novel eQTL Hotspots on Chromosomes 1 and 4
To obtain insights into the novel eQTL hotspot on chromosome 1 we identified 18 DEGs within ±5 cM of the eQTL hotspot located on chromosome 1 (Supplementary Table S4). Up-regulated DEGs included genes encoding the cupredoxin domain containing protein, arabinogalactan, the DREPP plasma membrane polypeptide family protein, and gibberellin 3β-hydroxylase (DWARF 18; D18) (Figure 5A). We also identified 27 DEGs within ±5 cM of the eQTL hotspot located on chromosome 4 (Supplementary Table S4). The five top ranking up-regulated genes with large fold-changes encoded monosaccharide transporter 1 (Figure 5B), the major sperm protein domain containing protein, glyceraldehyde 3-phosphate dehydrogenase, the plant invertase/pectinesterase inhibitor domain containing protein, and heavy metal transport/detoxification protein domain containing protein. We also identified a gene encoding O. sativa GROWTH-REGULATING FACTOR 12 (OsGRF12) as a significantly up-regulated DEG (see “Discussion”). Our list contained six down-regulated genes, including genes encoding phragmoplast-associated kinesin-related protein 1 and protein kinase- and ubiquitin domain containing proteins.
FIGURE 5. Comparison between differentially expressed genes (DEGs) and eQTL hotspots on chromosomes 1 and 4. (A) The differential expressions in non-deepwater- and deepwater-type RILs are plotted on the region of eQTL hotspots on chromosome 1 (genome position: 0–7.72 Mb) in log2-foldchanges. The red horizontal line indicates LOD = 3. (B) The differential expressions in non-deepwater- and deepwater-type RILs are plotted on the region of eQTL hotspots on chromosome 4 (genome position: 19.7–29.7 Mb) along the physical genome position in log2-foldchange. The red horizontal line indicates LOD = 3.
To further elucidate the gene regulatory networks that underlie expression variations between non-deepwater- and deepwater rice, promoters of the candidate target genes of the novel eQTL hotspots were investigated for over-represented motifs. Our analysis identified significant enrichment (Fisher’s exact test, p < 0.05) for an AP1 binding site (consensus sequence CCATTTTTAG) for candidate genes with trans-eQTLs at a physical position on chromosome 4 (Figure 6A); there were no significant motifs in 55 e-traits on chromosome 1. As for the promoters of candidate genes with trans-eQTL on chromosome 12, the enriched motifs were involved in ABA signaling including ABRE (ABA responsive element)- and ABF (ABRE-BINDING FACTOR) binding sites (Figure 6B). The second largest class of enriched motifs was related to MYB transcription factors, especially the MYB2 binding site motif [TAACT(G/C)GTT] and the MYB4 binding site motif [A(A/C)C(A/T)A(A/C)C]; the other enriched motif was the ARF (Auxin Response Factor) binding site motif (TGTCTC). Our observation suggests that these transcription factors play a role in the promotion of rapid internode elongation in deepwater rice.
FIGURE 6. Network visualization of enriched motifs and genes with trans-eQTLs hotspots. We performed promoter enrichment analysis for promoters of e-traits that feature trans-eQTL at physical position “c4.loc24” on chromosome 4 (A) and “c12.loc25.5” on chromosome 12 (B). Significantly enriched motifs were listed in the tables (p < 0.05). The enrichment of different types of motif is visualized in colors representing the motif class. The network represents the relationship between evaluated genes and the enriched motifs (yellow triangles). Network nodes (green) indicate genes. Motif classes are as follows: yellow green, AP1; magenta, ABA; black, ARF; and cyan, MYB. Edge thickness is proportional to the number of motifs that exist in the promoters of the respective genes. Note that the analysis showed no significant motifs in 55 e-traits on chromosome 1.
Discussion
Identification of eQTLs Controlling Gene Expression Networks Underlying Deepwater Responses
Transcript abundance, used as e-traits in genetical genomics studies, is highly heritable and genetically controlled (Jansen and Nap, 2001; Kliebenstein, 2009; Druka et al., 2010; Cubillos et al., 2012a; Chitwood and Sinha, 2013). Our aim was to map QTLs for gene expression that exist behind the regulatory control mechanisms underlying deepwater responses. Genetical genomics study using internode samples of RIL populations derived from non-deepwater (T65)- and deepwater rice (Bhadua) can be used as a model to gain a better understanding of the genotype-phenotype relationships elicited in response to submergence in monocot plants. Using genetic markers, microarray technology, and RIL population (Figure 1), we identified thousands of significant eQTLs (Figure 2).
Others (Brem et al., 2002; Schadt et al., 2003; Breitling et al., 2008) applied genetical genomics approaches to identify multiple e-traits that map to the same genomic locus; they are called eQTL hotspots. As such hotspots have been reported in plants such as Arabidopsis, rice, and soybean, it is expected that eQTL hotspots can be used for the identification of key regulators that elicit a wide range of downstream changes in gene expressions. However, the identification of significant biologically meaningful eQTL hotspots is hampered by undesirable variations from non-genetic factors such as uncontrolled environmental conditions and the physiological seed quality [for example, see (Breitling et al., 2008)]. To avoid the detection of spurious eQTL hotspots attributable to non-genetic variations we applied quantile-based permutation thresholds (Neto et al., 2012).
Detailed Evaluation of the eQTLs and the Novel trans-eQTL Hotspots to Infer Causal Regulatory Networks of Deepwater Responses
Our eQTL mapping highlighted that the three eQTL hotspots on chromosomes 1, 4, and 12 were possible candidate factors for internode elongation in deepwater rice (Figure 4). The eQTL hotspots on chromosome 12 co-localized with phQTL locating SK1/2 (Figures 3, 4), suggesting that SK1/2 are one of master regulator genes that controls the deepwater response. No significant phQTLs was observed on the locations of eQTL hotspots on chromosomes 1 and 4 (Figure 4 and Supplementary Figure S5). These inconsistencies could be explained by following two reasons. One reason is that resolution of the QTL analysis was not enough for detection of phQTL on chromosome 1 and 4; we detected small peaks of phQTL on both chromosome 1 and 4 although they were not statistically significant (Supplementary Figure S5). Especially, phQTL on chromosome 4 has been detected on QTL analysis for internode elongation focusing on response in early developmental stage (Nagai et al., 2012). Another reason is that these eQTL hotspots may contribute to other responses under submergence rather than internode elongation. Further QTL analysis with higher resolution is needed to evaluate these possibilities. On the other hand, no significant eQTL hotspot was detected at the location of phQTL at c1.loc40 on chromosome 1 (Figure 4 and Supplementary Figure S5). This imply that responsible gene of the phQTL at c1.loc40 could directly function on internode elongation rather than transcriptional regulation.
The trans-eQTL enrichment on chromosomes 1, 4, and 12 may indicate the location of a master regulator gene that controls the deepwater response. The D18 gene was detected as one of the DEGs in the eQTL hotspot on chromosome 1 (Figure 5A and Supplementary Table S4). D18 encodes O. sativa gibberellin 3 oxidase 2 (OsGA3ox2) known to control shoot elongation in rice by catalyzing the final step in the synthesis of bioactive gibberellins (GAs) (Itoh et al., 2001). It has been reported to be expressed in elongating internodes (Kaneko et al., 2003), and that D18 loss-of-function blocks deepwater-dependent internode elongation (Ayano et al., 2014). As the expression of D18 is significantly higher in deepwater phenotype than non-deepwater phenotype (Figure 5A), it is a strong candidate as the regulatory gene of the eQTL hotspot, suggesting that genes responsible for this eQTL hotspot may contribute to the specific gene expression resulting in internode elongation under submergence.
We identified a significant eQTL hotspot on chromosome 4 (Figure 5B). As one of the DEGs, we detected OsGRF12 in the eQTL hotspot on chromosome 4; its expression was higher in deepwater phenotype than non-deepwater phenotype (Figure 5B and Supplementary Table S4). Choi et al. (2004) reported that OsGRF12, a member of the plant-specific GRF transcriptional regulator family, expressed in internodes containing the intercalary meristem and in a part of the elongation zone, was induced by GA treatment. The overexpression of OsGRF10, a member of the GRF family with the highest similarity to OsGRF12, led to an increase in the number of internodes; a knock-down of both OsGRF3 and OsGRF4 resulted in a decrease in the internode length (Kuijt et al., 2014). Although its physical genome position is not consistent with the peak of phQTLs, it is possible that OsGRF12 is a candidate gene for the eQTL hotspot and it may play an important role in the transcriptional regulation of internode elongation under plant submergence.
We also looked for the candidate target genes of the three eQTL hotspots by searching for genes with trans-eQTLs at the physical position of each eQTL hotspot. Additional GO enrichment analysis of the candidate target genes of the novel eQTL hotspot on chromosome 12, locating SK1/2, showed that the GO functional categories photosynthesis, light reaction, monosaccharide metabolic process were significantly over-represented (Figure 4C). Thus, it is possible that SK1/2 mainly regulate these GO enriched genes for deepwater-rice specific submergence response. Similar to rice, in dicot plants, two Rumex species use ‘escape’ and ‘quiescence’ strategies upon submergence (Van Veen et al., 2013). Genome-wide transcriptome analysis by RNA-seq technology for comparison of the Rumex species suggested that light responsive signaling and plant hormone-related genes play a vital role in controlling the underwater elongation response (Van Veen et al., 2013). In our study, these genes were not enriched in the SK1/2 candidate target. Further approaches such as detailed RNA-seq analysis or identification of direct target of SK1/2 are needed to compare molecular mechanism of the submergence response in these monocot and dicot species.
Enhancement of Our Understanding of SK1/2-Dependent and Independent Pathways
Our results suggest that the downstream regulation of the eQTL hotspots on chromosomes 1 and 4 is independent, and that the target genes are partially regulated by SK1/2 (Figure 4B). The D18 gene, one of the candidate genes for eQTL hotspots on chromosome 1, harbors trans-eQTLs on chromosome 12 around the SK1/2 region (Supplementary Table S2C), suggesting that the expression of D18 is regulated by SK1/2. On the other hand, OsGRF12, one of the candidate genes for the eQTL hotspot on chromosome 4, may be expressed under the independent regulation of SK1/2 because it does not harbor trans-eQTLs at the physical position of the SK1/2 region (Supplementary Table S2C). Promoter enrichment analysis revealed that ABA, ARF, and MYB-related motifs were enriched in promoters of candidate target genes with the trans-eQTL hotspot on chromosome 12 (Figure 6B). In the AP2/ERF superfamily, SK1/2 belong to group VII including Arabidopsis hypoxia-responsive transcription factors such as RAP2.2 and RAP2.12 (Hattori et al., 2009). Given that RAP2.2 and RAP2.12 bind to the evolutionarily conserved motif consisting of the MYB binding site and the GC-rich region (Gasch et al., 2016), SK1/2 may transactivate these candidate target genes by binding to motif similar to the MYB-related motifs as well as to ABA and ARF.
Based on our observations, we present a hypothetical model of the signaling pathways leading to internode elongation in deepwater rice in response to submergence (Figure 7). When deepwater rice plants are submerged, ethylene accumulates, and this phytohormone induces SK1/2 expression via the rice ETHYLENE-INSENSITIVE 3 (EIN3)-like transcription factor (Hattori et al., 2009). SK1/2 promote the accumulation of GA and/or GA signal transduction, positively regulating elongation of the internodes and leaves. Based on our eQTL mapping findings we hypothesize the presence of 6 pathways for regulating the deepwater response (Figure 7).
FIGURE 7. Hypothetical model of signaling for internode elongation in submerged deepwater rice. Bold black arrows represent known pathways and dashed arrows the hypothetical pathways described in this study. Red and blue dashed arrows indicate the pathways regulated by genes in the eQTL hotspot on chromosome 1 and 4, respectively. Green dashed arrows indicate the pathway regulated by SK1/2 independent of the eQTL hotspots on chromosomes 1 and 4. Details are described in the “Discussion.”
The first pathway is regulated by genes in the eQTL hotspot on chromosome 1 and triggered by submergence and/or ethylene accumulation. The second is regulated by genes in the eQTL hotspot on chromosome 1 under direct or indirect regulation by SK1/2. The third is regulated by the D18 gene in the eQTL hotspot on chromosome 1 under direct or indirect regulation by SK1/2. The D18 gene product may directly regulate an increase in active GA. Since it is reported that the precursor level of active GAs is increased under submergence (Ayano et al., 2014), increase of D18 gene expression may contribute to acceleration of production of active GA. The fourth pathway is directly or indirectly regulated by SK1/2 and independent of the eQTL hotspots on chromosomes 1 and 4. The candidate 607 genes in this pathway may link between SK1/2 and submergence-induced internode elongation. The fifth is regulated by the gene in the eQTL hotspot on chromosome 4 under direct or indirect regulation by SK1/2. The sixth pathway is regulated by genes in the eQTL hotspot on chromosome 4 and triggered by submergence and/or ethylene accumulation. OsGRF12 may contribute to this pathway as a response gene of the eQTL hotspot on chromosome 4. Detailed molecular analyses such as regulation by GA or ethylene would help whether these candidate genes actually function in these pathways for internode elongation under submergence.
Conclusion
In this study, we present genome-wide eQTL mapping based on a microarray platform using RIL populations from T65- and Bhadua rice. Our analysis shows the existence of novel trans-eQTL hotspots on chromosomes 1, 4, and 12. The identification of cis- and trans-eQTLs yielded an array of candidate genes ready for further QTL analysis in structured populations with a deepwater response. Further QTL analyses with higher population size, marker density, and transcriptome under time-course submergence treatment will efficiently narrow down candidate genes for deepwater response. The association between trans-eQTLs and phQTLs in a given phenotype may provide a clue to candidate genes associated with a specific trait but does not necessarily yield candidates involved in the generation of the phQTLs. Such candidate genes will be useful for gene discovery and for studies on the deepwater response in rice. Integrating other omics data, e.g., proteomics and metabolomics (Fukushima et al., 2009; Fukushima and Kusano, 2014), with eQTL data from the same lines will greatly help to link hundreds of cellular molecules to eQTLs. Systematic and comprehensive studies may make possible the elucidation of complex regulatory networks underlying the genetic architecture of deepwater responses.
Availability of Supporting Data
All microarray data are available in the NCBI GEO database (Barrett et al., 2013) (accession no. GSE87702).
Author Contributions
TK, MA, and AF conceived and designed the experiments. TK, KN, YK, and YN performed the experiments. TK, MK, and AF analyzed the data. TK, KN, YK, YN, HY, MA, and AF contributed reagents/materials/analysis tools. TK and AF wrote the paper.
Conflict of Interest Statement
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
This work was supported by MEXT/JSPS KAKENHI (Grant Numbers 22119007, 24114001, 24114005, 26850024, 16K18565, 17K07663, 17H06473), SATREPS by JICA and JST, and Core Research for Evolutional Science and Technology by JST. We thank Ms. Midori Ito for technical assistance with RNA extraction, Ms. Yoko Niimi for RIL genotyping, Ms. Ritsuko Motoyama for technical support in our microarray analysis, Mr. Tomonori Noda for assistance with genotyping, Dr. Yasunori Ichihashi for guidance in promoter enrichment analysis, Ms. Ursula Petralia for editorial assistance, and Professor Julin Moloof for permitting us to re-use and modify QTLVizR (https://github.com/tiaho/QTL-Visualization). The preservation of the seed stocks of RILs were supported by the National Bio-Resource Project (NBRP) executed under the MEXT, Japan (http://www.shigen.nig.ac.jp/rice/oryzabase/top/top.jsp).
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017.01753/full#supplementary-material
FIGURE S1 | Schematic view of our experimental setup. Diagram showing how rice plants were submerged for measurement of the TIL and for RNA extraction. (A) Eight-leaf-stage plants of the RILs and their parents were grown in pots and completely submerged in water. For phenotype analysis, TIL was measured 7 days after submergence treatment and defined as the length between the uppermost and the basal node. (B) For microarray analysis, the RILs and their parents were submerged for 6 h. The region between 1 cm upper part from the uppermost and the basal node was used for RNA extraction. Leaves were removed from their nodes advanced to the extraction. Abbreviations: TIL, total internode length; PH, plant height; NEI, number of elongated internodes; and SAM, shoot apical meristem.
FIGURE S2 | Linkage map of T65/Bhadua RILs.
FIGURE S3 | Overview of genome-wide transcriptome data. (A) Results of normalization of the microarray data for the RILs and their parents. (B) Dendrogram of the RILs and their parents based on microarray data. Euclidean distances and the complete linkage method were used for hierarchical clustering.
FIGURE S4 | Statistical distribution of TIL, PH, and NEI among the RILs and their parents. (A) TIL, (B) PH, and (C) NEI are shown. White and black arrowheads represent the phenotypic values for T65 and Bhadua rice, respectively. Abbreviations: TIL, total internode length; PH, plant height; NEI, number of elongated internodes.
FIGURE S5 | Overview of phQTL data of the RILs. LOD statistics of phQTLs for TIL (A), PH (B), and NEI (C). Arrows point to the pseudo-marker or marker positions. The candidate regions of qRIE1 and qRIE12 (Kawano et al., 2008) are shown with blue lines. Abbreviations: TIL, total internode length; PH, plant height; NEI, number of elongated internodes; and RIE, rate of internode elongation. The red horizontal line indicates LOD = 3.
FIGURE S6 | Overview of eQTL hotspot size significance profile. The eQTL hotspots in the 12 rice chromosomes were derived with a quantile-based permutation approach (Neto et al., 2012). Data are presented essentially as shown in Figure 4A.
Footnotes
- ^ https://cran.r-project.org/
- ^ http://www.bioconductor.org/packages//2.10/bioc/html/Agi4x44PreProcess.html
- ^ https://bitbucket.org/jnmaloof/motifenrichment/
- ^ http://arabidopsis.med.ohio-state.edu/AtcisDB/bindingsites.html
- ^ https://github.com/tiaho/QTL-Visualization
- ^ http://shiny.rstudio.com/
References
Arends, D., Prins, P., Jansen, R. C., and Broman, K. W. (2010). R/qtl: high-throughput multiple QTL mapping. Bioinformatics 26, 2990–2992. doi: 10.1093/bioinformatics/btq565
Ashikari, M., and Matsuoka, M. (2006). Identification, isolation and pyramiding of quantitative trait loci for rice breeding. Trends Plant Sci. 11, 344–350. doi: 10.1016/j.tplants.2006.05.008
Ayano, M., Kani, T., Kojima, M., Sakakibara, H., Kitaoka, T., Kuroha, T., et al. (2014). Gibberellin biosynthesis and signal transduction is essential for internode elongation in deepwater rice. Plant Cell Environ. 37, 2313–2324. doi: 10.1111/pce.12377
Bailey-Serres, J., and Voesenek, L. A. (2008). Flooding stress: acclimations and genetic diversity. Annu. Rev. Plant Biol. 59, 313–339. doi: 10.1146/annurev.arplant.59.032607.092752
Bailey-Serres, J., and Voesenek, L. A. (2010). Life in the balance: a signaling network controlling survival of flooding. Curr. Opin. Plant Biol. 13, 489–494. doi: 10.1016/j.pbi.2010.08.002
Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. 41, D991–D995. doi: 10.1093/nar/gks1193
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300.
Bolon, Y. T., Hytenb, D. L., Orfa, J. H., Vancea, C. P., and Muehlbauer, G. J. (2014). eQTL networks reveal complex genetic architecture in the immature soybean seed. Plant Genome 7, 1–14. doi: 10.3835/plantgenome2013.08.0027
Bourgon, R., Gentleman, R., and Huber, W. (2010). Independent filtering increases detection power for high-throughput experiments. Proc. Natl. Acad. Sci. U.S.A. 107, 9546–9551. doi: 10.1073/pnas.0914005107
Breitling, R., Li, Y., Tesson, B. M., Fu, J., Wu, C., Wiltshire, T., et al. (2008). Genetical genomics: spotlight on QTL hotspots. PLOS Genet. 4:e1000232. doi: 10.1371/journal.pgen.1000232
Brem, R. B., Yvert, G., Clinton, R., and Kruglyak, L. (2002). Genetic dissection of transcriptional regulation in budding yeast. Science 296, 752–755. doi: 10.1126/science.1069516
Burow, M., Losansky, A., Muller, R., Plock, A., Kliebenstein, D. J., and Wittstock, U. (2009). The genetic basis of constitutive and herbivore-induced ESP-independent nitrile formation in Arabidopsis. Plant Physiol. 149, 561–574. doi: 10.1104/pp.108.130732
Chen, X., Hackett, C. A., Niks, R. E., Hedley, P. E., Booth, C., Druka, A., et al. (2010). An eQTL analysis of partial resistance to Puccinia hordei in barley. PLOS ONE 5:e8598. doi: 10.1371/journal.pone.0008598
Chitwood, D. H., Kumar, R., Ranjan, A., Pelletier, J. M., Townsley, B. T., Ichihashi, Y., et al. (2015). Light-induced indeterminacy alters shade-avoiding tomato leaf morphology. Plant Physiol. 169, 2030–2047. doi: 10.1104/pp.15.01229
Chitwood, D. H., and Sinha, N. R. (2013). A census of cells in time: quantitative genetics meets developmental biology. Curr. Opin. Plant Biol. 16, 92–99. doi: 10.1016/j.pbi.2012.11.003
Choi, D., Kim, J. H., and Kende, H. (2004). Whole genome analysis of the OsGRF gene family encoding plant-specific putative transcription activators in rice (Oryza sativa L.). Plant Cell Physiol. 45, 897–904. doi: 10.1093/pcp/pch098
Cubillos, F. A., Coustham, V., and Loudet, O. (2012a). Lessons from eQTL mapping studies: non-coding regions and their role behind natural phenotypic variation in plants. Curr. Opin. Plant Biol. 15, 192–198. doi: 10.1016/j.pbi.2012.01.005
Cubillos, F. A., Yansouni, J., Khalili, H., Balzergue, S., Elftieh, S., Martin-Magniette, M. L., et al. (2012b). Expression variation in connected recombinant populations of Arabidopsis thaliana highlights distinct transcriptome architectures. BMC Genomics 13:117. doi: 10.1186/1471-2164-13-117
Decook, R., Lall, S., Nettleton, D., and Howell, S. H. (2006). Genetic regulation of gene expression during shoot development in Arabidopsis. Genetics 172, 1155–1164. doi: 10.1534/genetics.105.042275
Drost, D. R., Benedict, C. I., Berg, A., Novaes, E., Novaes, C. R., Yu, Q., et al. (2010). Diversification in the genetic architecture of gene expression and transcriptional networks in organ differentiation of Populus. Proc. Natl. Acad. Sci. U.S.A. 107, 8492–8497. doi: 10.1073/pnas.0914709107
Drost, D. R., Puranik, S., Novaes, E., Novaes, C. R., Dervinis, C., Gailing, O., et al. (2015). Genetical genomics of Populus leaf shape variation. BMC Plant Biol. 15:166. doi: 10.1186/s12870-015-0557-7
Druka, A., Potokina, E., Luo, Z., Jiang, N., Chen, X., Kearsey, M., et al. (2010). Expression quantitative trait loci analysis in plants. Plant Biotechnol. J. 8, 10–27. doi: 10.1111/j.1467-7652.2009.00460.x
Durinck, S., Moreau, Y., Kasprzyk, A., Davis, S., De Moor, B., Brazma, A., et al. (2005). BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics 21, 3439–3440. doi: 10.1093/bioinformatics/bti525
Durinck, S., Spellman, P. T., Birney, E., and Huber, W. (2009). Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat. Protoc. 4, 1184–1191. doi: 10.1038/nprot.2009.97
Fukao, T., Xu, K., Ronald, P. C., and Bailey-Serres, J. (2006). A variable cluster of ethylene response factor-like genes regulates metabolic and developmental acclimation responses to submergence in rice. Plant Cell 18, 2021–2034. doi: 10.1105/tpc.106.043000
Fukushima, A., and Kusano, M. (2014). A network perspective on nitrogen metabolism from model to crop plants using integrated ’omics’ approaches. J. Exp. Bot. 65, 5619–5630. doi: 10.1093/jxb/eru322
Fukushima, A., Kusano, M., Redestig, H., Arita, M., and Saito, K. (2009). Integrated omics approaches in plant systems biology. Curr. Opin. Chem. Biol. 13, 532–538. doi: 10.1016/j.cbpa.2009.09.022
Gasch, P., Fundinger, M., Muller, J. T., Lee, T., Bailey-Serres, J., and Mustroph, A. (2016). Redundant ERF-VII transcription factors bind to an evolutionarily conserved cis-motif to regulate hypoxia-responsive gene expression in Arabidopsis. Plant Cell 28, 160–180. doi: 10.1105/tpc.15.00866
Gentleman, R. C., Carey, V. J., Bates, D. M., Bolstad, B., Dettling, M., Dudoit, S., et al. (2004). Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 5, R80. doi: 10.1186/gb-2004-5-10-r80
Haley, C. S., and Knott, S. A. (1992). A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. Heredity (Edinb) 69, 315–324. doi: 10.1038/hdy.1992.131
Hattori, Y., Miura, K., Asano, K., Yamamoto, E., Mori, H., Kitano, H., et al. (2007). A major QTL confers rapid internode elongation in response to water rise in deepwater rice. Breed. Sci. 57, 305–314. doi: 10.1270/jsbbs.57.305
Hattori, Y., Nagai, K., and Ashikari, M. (2011). Rice growth adapting to deepwater. Curr. Opin. Plant Biol. 14, 100–105. doi: 10.1016/j.pbi.2010.09.008
Hattori, Y., Nagai, K., Furukawa, S., Song, X. J., Kawano, R., Sakakibara, H., et al. (2009). The ethylene response factors SNORKEL1 and SNORKEL2 allow rice to adapt to deep water. Nature 460, 1026–1030. doi: 10.1038/nature08258
Hattori, Y., Nagai, K., Mori, H., Kitano, H., Matsuoka, M., and Ashikari, M. (2008). Mapping of three QTLs that regulate internode elongation in deepwater rice. Breed. Sci. 58, 39–46. doi: 10.1270/jsbbs.58.39
Holloway, B., Luck, S., Beatty, M., Rafalski, J. A., and Li, B. (2011). Genome-wide expression quantitative trait loci (eQTL) analysis in maize. BMC Genomics 12:336. doi: 10.1186/1471-2164-12-336
Ichihashi, Y., Aguilar-Martinez, J. A., Farhi, M., Chitwood, D. H., Kumar, R., Millon, L. V., et al. (2014). Evolutionary developmental transcriptomics reveals a gene network module regulating interspecific diversity in plant leaf shape. Proc. Natl. Acad. Sci. U.S.A. 111, E2616–E2621. doi: 10.1073/pnas.1402835111
Itoh, H., Ueguchi-Tanaka, M., Sentoku, N., Kitano, H., Matsuoka, M., and Kobayashi, M. (2001). Cloning and functional analysis of two gibberellin 3 beta -hydroxylase genes that are differently expressed during the growth of rice. Proc. Natl. Acad. Sci. U.S.A. 98, 8909–8914. doi: 10.1073/pnas.141239398
Jansen, R. C., and Nap, J. P. (2001). Genetical genomics: the added value from segregation. Trends Genet. 17, 388–391. doi: 10.1016/S0168-9525(01)02310-1
Jimenez-Gomez, J. M., Wallace, A. D., and Maloof, J. N. (2010). Network analysis identifies ELF3 as a QTL for the shade avoidance response in Arabidopsis. PLOS Genet. 6:e1001100. doi: 10.1371/journal.pgen.1001100
Jung, K. H., Seo, Y. S., Walia, H., Cao, P., Fukao, T., Canlas, P. E., et al. (2010). The submergence tolerance regulator Sub1A mediates stress-responsive expression of AP2/ERF transcription factors. Plant Physiol. 152, 1674–1692. doi: 10.1104/pp.109.152157
Kaneko, M., Itoh, H., Inukai, Y., Sakamoto, T., Ueguchi-Tanaka, M., Ashikari, M., et al. (2003). Where do gibberellin biosynthesis and gibberellin signaling occur in rice plants? Plant J. 35, 104–115. doi: 10.1046/j.1365-313X.2003.01780.x
Kawano, R., Doi, K., Yasui, H., Mochizuki, T., and Yoshimura, A. (2008). Mapping of QTLs for floating ability in rice. Breed. Sci. 58, 47–53. doi: 10.1270/jsbbs.58.47
Kende, H., Van Der Knaap, E., and Cho, H. T. (1998). Deepwater rice: a model plant to study stem elongation. Plant Physiol. 118, 1105–1110. doi: 10.1104/pp.118.4.1105
Keurentjes, J. J., Fu, J., Terpstra, I. R., Garcia, J. M., Van Den Ackerveken, G., Snoek, L. B., et al. (2007). Regulatory network construction in Arabidopsis by using genome-wide gene expression quantitative trait loci. Proc. Natl. Acad. Sci. U.S.A. 104, 1708–1713. doi: 10.1073/pnas.0610429104
Kliebenstein, D. (2009). Quantitative genomics: analyzing intraspecific variation using global gene expression polymorphisms or eQTLs. Annu. Rev. Plant Biol. 60, 93–114. doi: 10.1146/annurev.arplant.043008.092114
Kliebenstein, D. J., West, M. A., Van Leeuwen, H., Loudet, O., Doerge, R. W., and St Clair, D. A. (2006). Identification of QTLs controlling gene expression networks defined a priori. BMC Bioinformatics 7:308. doi: 10.1186/1471-2105-7-308
Kuijt, S. J., Greco, R., Agalou, A., Shao, J., t Hoen, C. C., Overnas, E., et al. (2014). Interaction between the GROWTH-REGULATING FACTOR and KNOTTED1-LIKE HOMEOBOX families of transcription factors. Plant Physiol. 164, 1952–1966. doi: 10.1104/pp.113.222836
Kurokawa, Y., Noda, T., Yamagata, Y., Angeles-Shim, R., Sunohara, H., Uehara, K., et al. (2016). Construction of a versatile SNP array for pyramiding useful genes of rice. Plant Sci. 242, 131–139. doi: 10.1016/j.plantsci.2015.09.008
Lander, E. S., Green, P., Abrahamson, J., Barlow, A., Daly, M. J., Lincoln, S. E., et al. (1987). MAPMAKER: an interactive computer package for constructing primary genetic linkage maps of experimental and natural populations. Genomics 1, 174–181. doi: 10.1016/0888-7543(87)90010-3
Loreti, E., Van Veen, H., and Perata, P. (2016). Plant responses to flooding stress. Curr. Opin. Plant Biol. 33, 64–71. doi: 10.1016/j.pbi.2016.06.005
Lowry, D. B., Logan, T. L., Santuari, L., Hardtke, C. S., Richards, J. H., Derose-Wilson, L. J., et al. (2013). Expression quantitative trait locus mapping across water availability environments reveals contrasting associations with genomic features in Arabidopsis. Plant Cell 25, 3266–3279. doi: 10.1105/tpc.113.115352
Maere, S., Heymans, K., and Kuiper, M. (2005). BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics 21, 3448–3449. doi: 10.1093/bioinformatics/bti551
Mickelbart, M. V., Hasegawa, P. M., and Bailey-Serres, J. (2015). Genetic mechanisms of abiotic stress tolerance that translate to crop yield stability. Nat. Rev. Genet. 16, 237–251. doi: 10.1038/nrg3901
Nagai, K., Hattori, Y., and Ashikari, M. (2010). Stunt or elongate? Two opposite strategies by which rice adapts to floods. J. Plant Res. 123, 303–309. doi: 10.1007/s10265-010-0332-7
Nagai, K., Kondo, Y., Kitaoka, T., Noda, T., Kuroha, T., Angeles-Shim, R. B., et al. (2014). QTL analysis of internode elongation in response to gibberellin in deepwater rice. AoB Plants 6, plu028. doi: 10.1093/aobpla/plu028
Nagai, K., Kuroha, T., Ayano, M., Kurokawa, Y., Angeles-Shim, R. B., Shim, J. H., et al. (2012). Two novel QTLs regulate internode elongation in deepwater rice during the early vegetative stage. Breed. Sci. 62, 178–185. doi: 10.1270/jsbbs.62.178
Nah, G., Im, J. H., Kim, J. W., Park, H. R., Yook, M. J., Yang, T. J., et al. (2015). Uncovering the differential molecular basis of adaptive diversity in three Echinochloa leaf transcriptomes. PLOS ONE 10:e0134419. doi: 10.1371/journal.pone.0134419
Nandi, S., Subudhi, P. K., Senadhira, D., Manigbas, N. L., Sen-Mandi, S., and Huang, N. (1997). Mapping QTLs for submergence tolerance in rice by AFLP analysis and selective genotyping. Mol. Gen. Genet. 255, 1–8. doi: 10.1007/s004380050468
Neto, E. C., Keller, M. P., Broman, A. F., Attie, A. D., Jansen, R. C., Broman, K. W., et al. (2012). Quantile-based permutation thresholds for quantitative trait loci hotspots. Genetics 191, 1355–1365. doi: 10.1534/genetics.112.139451
Pages, H., Aboyoun, P., Gentleman, R., and Debroy, S. (2016). Biostrings: String Objects Representing Biological Sequences, and Matching Algorithms. R Package Version 2.42.1. Available at: https://rdrr.io/bioc/Biostrings/.
Potokina, E., Druka, A., Luo, Z., Wise, R., Waugh, R., and Kearsey, M. (2008). Gene expression quantitative trait locus analysis of 16 000 barley genes reveals a complex pattern of genome-wide transcriptional regulation. Plant J. 53, 90–101. doi: 10.1111/j.1365-313X.2007.03315.x
Price, A. H., Cairns, J. E., Horton, P., Jones, H. G., and Griffiths, H. (2002). Linking drought-resistance mechanisms to drought avoidance in upland rice using a QTL approach: progress and new opportunities to integrate stomatal and mesophyll responses. J. Exp. Bot. 53, 989–1004. doi: 10.1093/jexbot/53.371.989
Quarrie, S. A., Laurie, D. A., Zhu, J., Lebreton, C., Semikhodskii, A., Steed, A., et al. (1997). QTL analysis to study the association between leaf size and abscisic acid accumulation in droughted rice leaves and comparisons across cereals. Plant Mol. Biol. 35, 155–165. doi: 10.1023/A:1005864202924
Ranjan, A., Budke, J., Rowland, S. D., Chitwood, D. H., Kumar, R., Carriedo, L. G., et al. (2016). eQTL regulating transcript levels associated with diverse biological processes in tomato. Plant Physiol. 172, 328–340. doi: 10.1104/pp.16.00289
Sakai, H., Lee, S. S., Tanaka, T., Numa, H., Kim, J., Kawahara, Y., et al. (2013). Rice Annotation Project Database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol. 54, e6. doi: 10.1093/pcp/pcs183
Salvi, S., and Tuberosa, R. (2015). The crop QTLome comes of age. Curr. Opin. Biotechnol. 32, 179–185. doi: 10.1016/j.copbio.2015.01.001
Sasidharan, R., Mustroph, A., Boonman, A., Akman, M., Ammerlaan, A. M., Breit, T., et al. (2013). Root transcript profiling of two Rorippa species reveals gene clusters associated with extreme submergence tolerance. Plant Physiol. 163, 1277–1292. doi: 10.1104/pp.113.222588
Schadt, E. E., Monks, S. A., Drake, T. A., Lusis, A. J., Che, N., Colinayo, V., et al. (2003). Genetics of gene expression surveyed in maize, mouse and man. Nature 422, 297–302. doi: 10.1038/nature01434
Siangliw, M., Toojinda, T., Tragoonrung, S., and Vanavichit, A. (2003). Thai jasmine rice carrying QTLch9 (SubQTL) is submergence tolerant. Ann. Bot. 91, 255–261. doi: 10.1093/aob/mcf123
Smyth, G. K. (2005). “Limma: linear models for microarray data,” in Bioinformatics and Computational Biology Solutions Using R and Bioconductor, eds R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, and W. Huber (New York, NY: Springer), 397–420. doi: 10.1007/0-387-29362-0_23
Swanson-Wagner, R. A., Decook, R., Jia, Y., Bancroft, T., Ji, T., Zhao, X., et al. (2009). Paternal dominance of trans-eQTL influences gene expression patterns in maize hybrids. Science 326, 1118–1120. doi: 10.1126/science.1178294
Tang, D. Q., Kasai, Y., Miyamoto, N., Ukai, Y., and Nemoto, K. (2005). Comparison of QTLs for early elongation ability between two floating rice cultivars with a different phylogenetic origin. Breed. Sci. 55, 1–5. doi: 10.1270/jsbbs.55.1
The Gene Ontology Consortium (2014). Gene ontology consortium: going forward. Nucleic Acids Res. 43, D1049–D1056. doi: 10.1093/nar/gku1179
Van Veen, H., Mustroph, A., Barding, G. A., Vergeer-Van Eijk, M., Welschen-Evertman, R. A., Pedersen, O., et al. (2013). Two Rumex species from contrasting hydrological niches regulate flooding tolerance through distinct mechanisms. Plant Cell 25, 4691–4707. doi: 10.1105/tpc.113.119016
Voesenek, L. A., and Bailey-Serres, J. (2015). Flood adaptive traits and processes: an overview. New Phytol. 206, 57–73. doi: 10.1111/nph.13209
Wang, J., Yu, H., Weng, X., Xie, W., Xu, C., Li, X., et al. (2014). An expression quantitative trait loci-guided co-expression analysis for constructing regulatory network using a rice recombinant inbred line population. J. Exp. Bot. 65, 1069–1079. doi: 10.1093/jxb/ert464
Wang, J., Yu, H., Xie, W., Xing, Y., Yu, S., Xu, C., et al. (2010). A global analysis of QTLs for expression variations in rice shoots at the early seedling stage. Plant J. 63, 1063–1074. doi: 10.1111/j.1365-313X.2010.04303.x
Wang, Y., Han, Y., Teng, W., Zhao, X., Li, Y., Wu, L., et al. (2014). Expression quantitative trait loci infer the regulation of isoflavone accumulation in soybean (Glycine max L. Merr.) seed. BMC Genomics 15:680. doi: 10.1186/1471-2164-15-680
West, M. A., Kim, K., Kliebenstein, D. J., Van Leeuwen, H., Michelmore, R. W., Doerge, R. W., et al. (2007). Global eQTL mapping reveals the complex genetic architecture of transcript-level variation in Arabidopsis. Genetics 175, 1441–1450. doi: 10.1534/genetics.106.064972
Xu, K., Xu, X., Fukao, T., Canlas, P., Maghirang-Rodriguez, R., Heuer, S., et al. (2006). Sub1A is an ethylene-response-factor-like gene that confers submergence tolerance to rice. Nature 442, 705–708. doi: 10.1038/nature04920
Yano, K., Takashi, T., Nagamatsu, S., Kojima, M., Sakakibara, H., Kitano, H., et al. (2012). Efficacy of microarray profiling data combined with QTL mapping for the identification of a QTL gene controlling the initial growth rate in rice. Plant Cell Physiol. 53, 729–739. doi: 10.1093/pcp/pcs027
Yilmaz, A., Mejia-Guerra, M. K., Kurz, K., Liang, X., Welch, L., and Grotewold, E. (2011). AGRIS: the Arabidopsis gene regulatory information server, an update. Nucleic Acids Res. 39, D1118–D1122. doi: 10.1093/nar/gkq1120
Keywords: quantitative trait locus, expression QTL, eQTL hotspots, ethylene response factor, submergence, abiotic stress, Oryza sativa
Citation: Kuroha T, Nagai K, Kurokawa Y, Nagamura Y, Kusano M, Yasui H, Ashikari M and Fukushima A (2017) eQTLs Regulating Transcript Variations Associated with Rapid Internode Elongation in Deepwater Rice. Front. Plant Sci. 8:1753. doi: 10.3389/fpls.2017.01753
Received: 18 July 2017; Accepted: 25 September 2017;
Published: 13 October 2017.
Edited by:
Avinash Mishra, Central Salt and Marine Chemicals Research Institute (CSIR), IndiaReviewed by:
Xusheng Wang, St. Jude Children’s Research Hospital, United StatesHanwei Mei, Shanghai Agrobiological Gene Center, China
Hiroki Saito, Kyoto University, Japan
Copyright © 2017 Kuroha, Nagai, Kurokawa, Nagamura, Kusano, Yasui, Ashikari and Fukushima. 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) or licensor 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: Takeshi Kuroha, tkuroha@m.tohoku.ac.jp Atsushi Fukushima, atsushi.fukushima@riken.jp