Identifying candidate genes and patterns of heat-stress response in rice using a genome-wide association study and transcriptome analyses

Yingxue Yng,Cho Zhng,De Zhu,Huiying He,Zhorn Wei,Qioling Yun,Xioxi Li,Xu Go,Bin Zhng,Hongsheng Go,b,Bo Wng,Shuimin Co,Tinyi Wng,Yuhu Li,Xiomn Yu,Longbio Guo,Gunjing Hu,*,Qin Qin,,*,Lingung Shng,*

a Shenzhen Branch,Guangdong Laboratory of Lingnan Modern Agriculture,Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs,Agricultural Genomics Institute at Shenzhen,Chinese Academy of Agricultural Sciences,Shenzhen 518120,Guangdong,China

b College of Agriculture,Ludong University,Yantai 264025,Shandong,China

c China National Rice Research Institute(CNRRI),Chinese Academy of Agricultural Sciences,Hangzhou 311401,Zhejiang,China

Keywords:Heat stress Rice GWAS Alternative splicing Gene pyramiding

ABSTRACT Because high temperatures impair rice production,it is desirable to elucidate the regulatory mechanisms involved in rice response to heat stress.The objectives of this study were to identify candidate genes and characterize their response patterns during rice adaptation to high temperatures at the seedling stage.Ten heat-associated quantitative-trait loci were identified in a genome-wide association study.Comparison of transcript abundances in heat-sensitive and heat-tolerant rice pools under heat stress revealed approximately 400 differentially expressed genes.The expression of genes from heatsensitive accessions changed more than those from heat-tolerant accessions under heat stress.Alternative splicing(AS)events responded to heat stress in rice.The types of AS variants significant different between the heat-sensitive and heat-tolerant accessions.Expression patterns differing between the heat-sensitive and heat-tolerant accessions were identified for genes known to be involved in heat stress.We identified eleven genes associated with rice heat stress response.These genes could be pyramided to breed heat-tolerant rice accessions.

1.Introduction

Rice is the main food source for more than half of the world’s population.With global warming,rice grain yield reduction by high-temperature stress will become increasingly severe[1].A temperature increases of only a few degrees Celsius can damage the reproductive organs of crops,reducing pollen viability and pistil fertility,and reduce fruit setting rate,resulting in partial or total crop loss[2-4].Seedling growth too influences plant morphology.If rice seedlings encounter temperatures above 35°C,the types and contents of leaf proteins can be greatly altered,and the plant height,biomass,and other characteristics may be impaired[5].Thus,the development of high temperature-tolerant rice accessions contributes to safe and sustainable food production and is urgently needed to deal with global warming[6].

In response to high-temperature stress,plants use multiple mechanisms and regulatory networks to control the expression of genes involved in physiological and biochemical adaptation.The heat-shock transcription factor Hsf is a key factor in activating the high-temperature transcriptional regulation network.The heat-shock protein Hsp is regulated by Hsf and is induced in response to high-temperature stress.Silent mutations in the HsfA1 genes of Arabidopsis and potato can activate high temperature-responsive genes.OsHSP70CP1,a gene that encodes a protein that is localized to the rice chloroplast,is up-regulated and expressed under high temperature stress,and is necessary for the development of chloroplasts at high temperatures[7].Li et al.[8]identified in African rice(Oryza glaberrima)the quantitative trait locus(QTL)TT1,which encodes a proteasomeα2 subunit involved in the degradation of ubiquitinated proteins.Expression of TT1 protects cells by degrading cytotoxic denatured proteins,increasing heat tolerance.Shen et al.[9]found that overexpression in rice of the receptor-like kinase ERECTA(ER)increased heat tolerance,with plants carrying a silent mutation in the ER gene showing lower heat tolerance.Wang et al.[10]showed that the DEAD-box RNA helicase TOGR1 is a molecular chaperone that can increase rice growth by maintaining the steady state of rRNA under high-temperature stress conditions.In addition to the genes mentioned above,many alleles of heat-responsive genes in rice remain to be discovered,and the molecular mechanism of high-temperature stress regulation has not been studied in depth.

Genome-wide association study(GWAS)is a powerful tool for characterizing the genetic architecture of complex traits[11].Using diverse germplasm collections,GWAS has been shown to be useful for the high-throughput identification of QTL and for characterizing the genetic architecture of complex traits including yield[12-14],stress tolerance[15,16],and tissue development status[17-19].The method involves searching the genomes from many different accessions of plants and sorting out genetic markers that can be used to predict the presence of a phenotype.

Alternative splicing(AS)occurs during eukaryotic precursor mRNA(pre-mRNA)processing and generates multiple mature mRNA isoforms from the same gene,potentially giving rise to functionally different proteins.AS adds a regulatory layer by which plants can control gene expression patterns by changing the proportions of productive and unproductive variants[20-23].Types of AS events include cassette exon(Cassette),intron retention(IR),alternative 5′splice site(A5SS),and alternative 3′splice site(A3SS).High-throughput RNA sequencing(RNAseq)has advanced AS research in plant science by allowing the identification of previously unknown transcript isoforms and the assessment of dynamic changes in the full complement of transcript isoforms during development or in response to environmental cues[24-26].AS is sensitive to stress and may represent a rapid response to such conditions[27-30].Many genes are alternatively spliced at high temperature,including HSFA2d in rice,HSFA2 and HSFA7d in Arabidopsis,and HSFA2 in grape[31-33].

The aim of the present study was to investigate rice response to heat stress.So we conducted a genome-wide association study in a panel of accessions and comparison between heat-tolerant and heat-sensitive accessions by RNA-seq data.Pattern changes of AS events were analyzed in response to heat stress.The promising candidate genes were identified in our study and they could further be used for genetic breeding of beneficial rice accessions adapted to high temperature environments.And these findings will shed light on our understanding of the underlying genetic bases for the mechanisms that regulate the heat stress response in rice.

2.Materials and methods

2.1.Plant materials

Genomic DNA sequencing information from a panel of 221 rice accessions was collected for GWAS.The descriptions of the 210 accessions were from the 3000 Rice Genomes Project[34]and the remaining 11 were from the study of Wang et al.[35].These lines encompassed domesticated lines of Oryza sativa(japonica,indica,and aus)(Table S1)which were subjected to phenotyping,genotyping and subsequent GWAS analysis.Three heat-tolerant rice accessions;Baxiang(BX),Malaihong(MLH),and Honggenghangu(HGHG),and three heat-sensitive rice accessions;Laohongdao(LHD),Ailuyu(ALY),and Weiguo(WG)from the collection,together with Nagina 22(N22),93-11,and Nipponbare(NIP),were selected for heat stress-response phenotyping and transcriptomic analysis.

2.2.Heat stress treatment and phenotyping

Seeds were immersed in distilled water for three days until germination and then grown in seed trays for nine days under the control conditions(16 h light at 30 °C and 8 h darkness at 28 °C).The seedlings were then subjected to heat-stress treatment of 45 °C for 45 h with light,and were then returned to the control conditions for recovery.For phenotyping,survival scores were recorded after a five-day recovery period from the heat-stress treatment.Numerical survival scores ranging from 0 to 5 were assigned based on the percentage of surviving plants.The survival percentages corresponding to scores 0-5 are presented in Table S2.For heat stress,a score of 0 indicates the lowest tolerance and a score of 5 indicates the highest tolerance.More than 30 seedlings per accession were used in heat-stress treatments.Experiments were performed with at least three biological replicates for each of the seedlings.

2.3.Genome-wide association study(GWAS)

To identify loci associated with heat-stress survival in rice at the seedling stage,we performed GWAS analysis using the 221 accessions,which SNPs information were retrieved from the 3000 Rice Genomes Project and Wang et al.[35].The Efficient Mixed-Model Association eXpedited(EMMAx)[39]was used for GWAS analysis.SNPs with MAF＜0.05 and a missing-data rate＞0.2 were discarded,and the remaining 3,515,687 SNPs were used for GWAS.To present clearer association analyses and subsequent results,we calculated the number of independent SNPs using the Genetic Type I error method[40],and used 0.05/n(independently inherited SNP)as the threshold.The R(version 3.5.0)software package qqman was used to generate Manhattan plots[41].Candidate genes were detected in 200-kb regions upstream and downstream of SNP[42].

2.4.RNA isolation,library preparation,and sequencing

RNA was extracted from leaf samples collected before and after three hours of heat-stress treatment.RNA extraction using TRIzol reagent(Invitrogen,Carlsbad,CA,USA)was performed according to the manufacturer’s recommendations.Sequencing libraries were prepared using the NEBNext Ultra RNA Library Prep Kit for Illumina(New England Biolabs,Ipswich,MA,USA)following the manufacturer’s instructions.The Illumina NovaSeq(Illumina,San Diego,CA,USA)platform was used to sequence the libraries,which produced 150-base paired-end reads.Three replicates were sequenced for each sample.Clean reads were aligned to the(MSU v7.0)reference genome using hisat 2.0[43].Fragments per kilobase of exon per million mapped reads(FPKM)values were calculated by Stringtie[44](version 2.2.0)to determine gene expression levels,and a gene was considered to be expressed if FPKM≥1 in each library.

2.5.Differential gene expression analysis

To identify candidate genes whose expression responded to heat stress,differential gene expression analysis was performed with the R package DESeq2(1.18.0)using a model based on the negative binomial distribution.Three biological replicates per treatment were used in the analysis.The false discovery rate(FDR)was controlled by the adjusted P-value using the Benjamini and Hochberg approach[45].Genes with an adjusted P-value＜0.05 and|log2(fold change)|＞1 were considered to be differentially expressed genes(DEGs).

2.6.Functional enrichment analysis

GO(Gene Ontology)enrichment analysis was performed using the R package GOseq,with a corrected P-value of 0.05 as the threshold for significantly enriched GO terms.KEGG(Kyoto Encyclopedia of Genes and Genomes)pathway annotation and enrichment analysis of the DEGs were performed using KOBAS software[46].GO analysis for the GWAS was performed using the website https://systemsbiology.cau.edu.cn/agriGOv2/.

2.7.Alternative splicing(AS)analysis

The detection and comparison of AS events were conducted using the comprehensive alternative splicing hunting method(CASH v2.2.1)[47].For AS detection,mapped bam files after hisat2 alignment from replicated samples of the same accession were merged as the input for CASH.To compare AS events between groups of heat-sensitive and heat-tolerant accessions,an isoform was considered present in a group if it was detected in at least two accessions.It was identified as a differentially spliced gene(DSG)if it was present in at least two accessions of a group and absent in the other group with a false discovery rate(FDR)of＜0.05.

2.8.Quantitative PCR(qPCR)

qPCR analysis was performed as described previously[48].Briefly,cDNA was synthesized in vitro from 2μg of total RNA using SuperScriptII reverse transcriptase(Thermo Fisher Scientific,Waltham,MA,USA),and qPCR amplifications were performed using Brilliant III SYBR Green QPCR Master Mix(Agilent).Using OsACTIN as the internal control,gene expression levels were quantified by the comparative 2-△△CTmethod[48].The gene-specific primers used are presented in Table S3.

2.9.Candidate gene identification

To narrow the set of promising candidate genes associated with heat stress tolerance,we performed gene haplotype analysis.If a SNP identified in the GWAS was found in a gene coding region,genotype analysis was performed to identify genes with significant differences in phenotypes(survival scores)between genotypes,which could imply that they are candidate genes associated with heat-stress response.The gene expression profiles in the transcriptomic data for these genes were then examined.Genes that showed significant differences in gene expression in response to heat stress among most rice accessions were selected.Genes with amino acid changes were also assigned as candidates.If a SNP identified in the GWAS was found in the non-coding region of a gene,the four genes surrounding the SNP position were examined by haplotype analysis of SNPs to identify genes with significant differences in relevant phenotypes between the variation categories.A gene was identified as a candidate if it displayed significant expression differences among most rice accessions in response to heat-stress treatment and significant expression differences between heatsensitive and heat-tolerant accessions.The website RiceVarMap v2.0(https://ricevarmap.ncpgr.cn/)was used to identify SNP variations in candidate genes in the 221 accessions.

3.Results

3.1.Rice heat stress-response phenotyping under controlled growth conditions

Heat tolerance of the 221 accessions is summarized in Table S2.The five heat-tolerant and four heat-sensitive accessions selected are shown in Fig.1A.The five heat-tolerant accessions,93-11,BX,HGHG,MLH and N22 that are labeled in red showed high survival scores ranging from 4.33 to 4.67,and there were no significant changes in survival before and after heat stress treatment,whereas the four heat-sensitive accessions ALY,WG,NIP and LHD that are labeled in white showed lower survival scores ranging from 0.33 to 0.67(Fig.1B).These two groups of accessions were used for the subsequent RNA-seq experiments.To determine the suitable time duration for heat stress treatment in the RNA-seq experiments,qPCR analysis was performed using the heattolerant accession 93-11 and the heat-sensitive accession NIP.For this purpose,five known heat-responsive genes;TOGR1(LOC_Os03g46610),Hsp70cp1(LOC_Os05g23740),TT1(LOC_Os02 g02880),SNAC3(LOC_Os01g09550),and HsfB2c(LOC_Os09g35790)were selected for qPCR analysis using leaf samples collected at six time points(0,0.5,1,3,8,and 24 h)after heat-stress treatment.As shown in Fig.S1A and S1B,the relative expression patterns showed peaks at 1,3,and 8 h after treatment,consistent with previous reports.We selected a time point of 3 h post-treatment for the RNA-seq experiments.

3.2.Association mapping identified potential heat stress-responsive loci

We set the range of 200 kb upstream and downstream of the associated SNP as a QTL.As a result,41 significantly associated SNPs were identified in 10 QTL(HT1 to HT10)for heat stress tolerance(Fig.1C;Table S4).In genomic intervals 200 kb on either side of these SNP sites,670 candidate genes(including 432 transposable elements associated genes)were detected in the QTL regions,ranging from 44 in HT10 to 108 in HT5(Table S5).The known heatresponsive genes OsHPS17.7(LOC_Os03g16040),located in HT4,and OsSIZ2(LOC_Os03g50980),located in HT5,were detected among the candidate genes.GO of candidate genes is described in Table S6.Photosynthesis(GO:0015979)and membrane(GO:0016020)were associated with heat stress[49,50],which showed the reliability of the GWAS analysis.

3.3.Differentially expressed genes(DEGs)in response to heat stress

Using the NIP genome as a reference,1.37 billion reads were aligned.A total of 40,793 transcripts from the nine accessions were identified.RNA-seq data was validated for gene expression by qPCR(Fig.S2).The relative expression of 10 selected genes in heat-tolerant 93-11 and heat-sensitive NIP showed the same pattern whether measured by qPCR or RNA-seq.

To identify genes that were differentially expressed in response to heat stress treatment,gene expression was compared in leaves before and after heat treatment in each of the nine rice accessions.The up-regulated and down-regulated genes in each accession are summarized in Fig.2A.A total of 1891 common DEGs were detected in the five heat-tolerant accessions(Fig.2B)and 3239 common DEGs were detected in all four heat-sensitive accessions(Fig.2C).Comparison of DEGs between the heat-sensitive and heat-tolerant accessions revealed 1990 DEGs in the sensitive and 642 in the tolerant accessions and 1249 in both(Fig.2D).To better screen candidate genes related to high temperature tolerance,we further filtered the above DEG data by considering the effects of the heat-stress treatment per se.DEGs from the first step were sorted into three categories;(1)DEGs that were uniquely expressed in the four sensitive accessions,(2)DEGs that were uniquely expressed in the five heat-tolerant accessions,and(3)their intersections(DEGs expressed in both the heat-sensitive and heat-tolerant accessions).This filter revealed DEGs induced by heat treatment rather than by the accession itself.As expected,fewer DEGs were identified in the heat-sensitive and heat-tolerant accessions:1691 in the sensitive,604 in the tolerant accessions and 1246 in both.A heat map of the DEGs that were expressed mainly in the four sensitive accessions and those that were expressed mainly in the five tolerant accessions is shown in Fig.2E.The two pools showed differing gene expression patterns.Approximately 400 DEGs in the heat map are listed in Table S7.These genes represent accession-specific DEGs associated with heat stress in rice.This analysis showed that the rice leaf transcriptome undergoes dynamic changes in response to heat stress.

Fig.1.The stress-responsive phenotypes of nine rice accessions.(A)The upper and lower panel show the phenotypes respectively before and after heat-stress treatment.Five heat-tolerant accessions:93-11,BX,HGHG,MLH and N22 are labeled in red.Four heat-sensitive accessions:ALY,WG,NIP and LHD,are labeled in white.The white line at the bottom right corner represents the scale of 3 cm.(B)Bar plot of survival scores.Significant differences between control and heat treatment were calculated by Student’s t-test,*,P＜0.05;**,P＜0.01;***,P＜0.001;ns indicates no significant difference.Values are presented as means±SD.(C)Manhattan plot identifying candidate loci involved in heat stress regulation.

Fig.2.Transcriptome analysis of differentially expressed genes(DEGs)in nine rice accessions.(A)Summary of up-regulated and down-regulated gene numbers in nine rice accessions.(B)Venn diagram of DEGs in the heat-tolerant accessions.(C)Venn diagram of DEGs in the heat-sensitive accessions.(D)Comparison of DEGs numbers in tolerant accessions(left circle)and sensitive accessions(right circle).(E)Heat map of gene expression analysis in nine rice accessions with or without heat treatment.DEGs were chosen that were expressed in heat-tolerant accessions or in heat-sensitive accessions.The gene list is presented in Table S7.

3.4.GO and KEGG pathway annotation analysis of the rice DEGs

The GO classification was analyzed based on the above filtered DEGs,and all genes were classified into the three main GO categories of molecular function(MF),cellular component(CC),and biological process(BP).For the GO annotation of the heat-stress treatment DEGs,regardless of whether the DEGs were from the heat-sensitive or heat-tolerant accessions,oxidoreductase activity(GO:0016491),apoplast(GO:0048046),and oxidation-reduction process(GO:0055114)were the primary subcategories identified in MF,CC,and BP,respectively(Fig.S3A).In the GO annotation of DEGs sorted by heat-stress treatment,the dominant terms in the heat-sensitive accessions were oxidoreductase activity(GO:0016491),cell periphery(GO:0071944),and protein folding(GO:0006457)in MF,CC,and BP,respectively(Fig.S3B).The dominant terms in the heat-tolerant accessions were cyclin-dependent protein serine/threonine kinase regulator activity(GO:0016538),cyclin-dependent protein(GO:0000307),and regulation of cyclin-dependent protein serine/threonine kinase activity(GO:0000079)in MF,CC,and BP,respectively(Fig.S3C).These GO results suggest that the DEGs from the heat-sensitive and heattolerant accessions in the different GO categories may result from the differences in heat tolerance between the two groups.In the intersection of DEGs from the heat-sensitive and heat-tolerant rice lines,the GO annotation indicated that oxidoreductase activity(GO:0016491),cell wall(GO:0005618),and response to heat(GO:0009408)were the primary terms in the MF,CC,and BP categories,respectively(Fig.S3D).The GO analysis results suggest that one of the possible mechanisms by which rice perceives and responds to heat stress is its effect on oxidoreductase activity.For the corresponding KEGG pathway annotation of the DEGs from the heat-stress treatment in all accessions,the most enriched term was biosynthesis of secondary metabolites(Fig.S4A).In the above KEGG pathway annotation of the DEGs sorted by heat-stress treatment,the main enriched terms for the heat-sensitive accessions,heat-tolerant accessions,and their intersection were metabolic pathways,biosynthesis of secondary metabolites and metabolic pathways,respectively(Fig.S4B,C and D).Cai et al.[51]reported that a class of secondary metabolites,the flavonoids,and their biosynthetic pathways function in heat-induced defense responses,suggesting that the thermotolerance mechanism involves metabolic pathways.

3.5.Profiling the AS landscape in the response to heat stress

AS events were detected by RNA-seq analysis in nine rice accessions and categorized into eight classes:Cassette,multiple adjacent cassette exons(Cassette_multi),IR,A5SS,A3SS,alternative start exon(AltStart),alternative end exon(AltEnd)and mutually exclusive exons(MXE)(Fig.3A).The most abundant AS type was IR,consistent with previous reports in plants[52](Fig.3B,all events in green).The most abundant AS type differentially affected by heat stress was Cassette(Fig.3B,significant events in black).Significantly differentially affected events accounted for 39.5% of all detected AS events differed significantly between the heatsensitive and heat-tolerant accessions(defined as differentially spliced genes;DSGs),suggesting that alternative splicing is involved in heat-stress response.Nearly all of the detected DSGs also showed differential expression changes in response to heat stress(Table S8).A total of 39 known heat-responsive genes including Hsfs and Hsps were identified as DSGs[53].Their expression profiles presented in a heat map showed expression patterns differing before and after heat stress treatment in most rice accessions(Fig.3C).

GO analysis was performed for genes from the overlap between the DSGs and DEGs in the response to heat stress.The most significant GO terms included RNA binding(GO:0003723)and hydrolase activity(GO:0016787)in MF,intracellular(GO:0005622)and intracellular part(GO:0044424)in CC,and organelle organization(GO:0006996)in BP(Fig.S5A,B,and C).KEGG pathway annotation indicated that spliceosome,mRNA surveillance,and aminoacyltRNA biosynthesis were the most highly enriched pathways(Fig.S5D).These results suggest that regulation of alternative splicing could act together with gene expression to modulate biological functions during the response to heat stress.

3.6.Analysis of selected candidate genes by combining GWAS and transcriptome data

An analysis that combined GWAS and transcriptome data was performed to further reduce the number of promising candidate genes.A set of 11 candidate genes were identified by examining SNP variation and gene expression,and all of them were promising targets in the regulation of heat stress.Among these genes,nine(located in HT1,HT2,HT3,HT4,and HT5)showed significantly different phenotypes in the variation categories and gene expression as determined by the RNA-seq data in the response to heat stress(Figs.S6,S7,and S8).The individual annotations can be found in the candidate gene list in Table S5.We focused our attention on the other two genes,LOC_Os03g16460 and LOC_Os05g07050,owing to their functional annotation and high expression(Fig.4).LOC_Os03g16460,located in HT6,encodes a putative uncharacterized protein.LOC_Os05g07050,located in HT4,was annotated as encoding a putative pre-mRNA processing splicing factor 8.Investigating the SNPs in these two genes in the collection of 221 accessions(Fig.4A,D)revealed seven mutation sites classified into five categories for LOC_Os03g16460 and three mutation sites classified into two categories for LOC_Os05g07050 in comparison to the reference genotype(Ref)(Fig.4C,F).This finding also showed that there were predicted amino acid changes in the encoded proteins due to the mutations.For LOC_Os03g16460,plants carrying a single mutation in the H1 or H7 categories displayed relatively low heat tolerance(low survival scores),whereas plants in the category with all seven mutated sites showed the highest heat tolerance(Fig.4B).Similarly,for LOC_Os03g16460,plants in the category with all three mutation sites displayed higher heat tolerance than plants in category H3 with a single mutation(Fig.4E).Categories with more mutation sites tended to have higher heat tolerance,mostly in indica rice accessions,whereas categories with fewer mutation sites tended to have reduced heat tolerance,mostly in japonica accessions.The indica or japonica preferences agree with the observation that indica rice generally has higher heat tolerance than does japonica rice.The expression of the LOC_Os03g16460 and LOC_Os05g07050 genes from the RNA-seq data is shown in Fig.4G and H.The two genes showed significant differences in gene expression due not only to heat stress per se,but to accession.The expression in heat-sensitive accessions differed from that in heat-tolerant accessions,suggesting their potential role in perceiving high temperature signals.

4.Discussion

4.1.Strategies for detecting heat stress-associated QTL

Several strategies can be used for the detection of heat stress-associated QTL.Many QTL have been identified by linkage mapping;however,fine-mapping requires the time-consuming establishment of large segregating populations[54].In addition,QTL mapping is based on bi-parental populations that lack genetic diversity;this deficiency can be compensated by GWAS analyses even though this method also has the drawbacks of false associations and rare alleles.Transcriptome analysis has facilitated the identification of candidate genes that function the response to heat stress,and can limit the range of candidate genes detected by QTL mapping or GWAS.From the high-throughput RNA-seq datasets,we could get information from analysis of alternative splicing,since AS provides a post-transcriptional mechanism in response to stress for various functions[23].Alternative splicing provides another layer for heat-stress regulation[55].AS analysis may identify promising candidate loci.Previously reported genes encoding heat-shock proteins and heat-stress transcription factors,such as HSP82B and HsfA3,were identified in the present study as DSGs in response to heat stress treatment.

Fig.3.Analysis of alternative splicing(AS)after heat stress treatment.(A)Categorization of AS events.Blue boxes represent exons.Yellow boxes represent partial exons to be excluded or included.Dotted lines connect exons.(B)Detection of AS events in nine accessions.All AS events(green)and significant splicing events between heat-sensitive and heat-tolerant accessions(black)are presented.(C)Heat map of expression profiles for 39 known heat-responsive genes.

Studies[56,57]have reported using GWAS for the identification of candidate genes associated with heat stress tolerance.These studies used multiple genetic populations in their analyses.In the present study,we collected gene expression data from 221 accessions of Oryza sativa that included accessions from the japonica,indica and aus subgroups that exhibited diverse traits such as tolerance to different abiotic stresses.Most of these studies focused their analyses on the flowering stage in rice,because decreased spikelet fertility can eventually affect crop yield.However,in the present study we investigated the effects of heat stress on rice during the seedling stage.For this reason,our candidate QTL identified by GWAS mapping are different from those reported previously.The known heat-responsive genes OsHPS17.7(LOC_Os03g16040)and OsSIZ2(LOC_Os03g50980)were also detected in our analysis,indicating the reliability of the present study.

4.2.Promising candidate genes in heat stress

Our study has demonstrated that the combined results from GWAS and transcriptome analyses can reduce the number of target genes.We were able by this means to identify promising candidate genes involved in tolerance to heat stress.LOC_Os03g16460 is predicted to encode an uncharacterized protein.It was previously named OsFes1A because it shares a relatively low degree of sequence identity with OsFes1C,a gene that encodes a potential nucleotide exchange factor for OsBiP1 that responds to endoplasmic reticulum and salt stress[58].In Arabidopsis,AtFes1B(AT3G53800),together with AtFes1A(At3g09350)and AtFes1C(At5g02150),are three AtFes1 genes belonging to the ARM(armadillo repeat)superfamily[59].Expression of these genes is strongly induced by heat stress and is repressed in the athsfa1a athsfa1b double mutant[60].In the present analysis,LOC_Os03g16460 also displayed strong expression by RNA-seq in response to heat stress,suggesting its possible involvement in heat-stress regulation.The finding from analyzing the variable sites in LOC_Os03g16460 indicates that the nucleotide change at position 9,076,778 may influence high-temperature tolerance in rice.LOC_Os05g07050 was annotated as encoding a putative premRNA-processing-splicing factor 8.The orthologous gene AT1G80070 in Arabidopsis encodes a protein that affects premRNA splicing and is essential for embryonic development[61]because it enables RNA binding and protein binding.The orthologous genes in Brachypodium,maize,and sorghum also act as RNA splicing factors.Thus,LOC_Os05g07050 may have a splicing factorassociated function in rice.Our study showed that alternative mRNA splicing patterns undergo dynamic changes in response to heat stress,and because splicing factors are involved in these processes,it is possible that LOC_Os05g07050 could be associated with the heat stress response via regulation of AS by splicing factors.By analyzing the variable sites in LOC_Os05g07050 and the related heat-response phenotypes in the variation categories,we found that the A to C transversion at position 3,711,495 results in an amino acid change from Ile to Leu,and increases heat tolerance compared to the control(Ref),suggesting that this amino acid change may influence high-temperature tolerance in rice.The direct roles of the genes identified in the seedling response to heat stress in this study have not been reported to date.Further research may reveal how these genes are associated with heatstress regulation.

Fig.4.Variation analysis of candidate genes LOC_Os03g16460 and LOC_Os05g07050 in 221 rice accessions.(A)Gene schematic of LOC_Os03g16460 showing SNP variation sites.(B)Heat-responsive phenotypes in the variation categories.Ref represents control of genotype.H1 to H7 represent different haplotypes.(C)Variation categories and relevant amino acid changes of LOC_Os03g16460 and the distribution between indica and japonica.(D)Gene schematic of LOC_Os05g07050 showing SNP variation sites.(E)Heat-responsive phenotypes in the variation categories of LOC_Os05g07050.(F)Variation categories and relevant amino acid changes of LOC_Os05g07050 and the distribution between indica and japonica.(G)Gene expression of LOC_Os03g16460 from RNA-seq data.(H)Gene expression of LOC_Os05g07050 from RNA-seq data.Significant differences of phenotypes between variation categories were calculated by Mann Whitney U test.Significant differences between control and heat treatment in G and H were calculated by Student’s t-test.No significant difference was labeled with ns.*,P＜0.05;**,P＜0.01;***,P＜0.001;****,P＜0.0001.Values are presented as means±SD.

4.3.Effects of candidate gene pyramiding in heat stress

Gene pyramiding,in which multiple elite genes are combined in a single genotype,is a breeding technique used to obtain a desired phenotype.One way to improve rice tolerance to heat stress is by incorporating multiple heat tolerance genes into existing elite accessions.In the present study,we identified 11 highly promising candidate genes associated with heat stress.We examined the stacking effect of these genes to select for the desired genotypes(Fig.5).Four gene pyramiding combinations were found to increase heat tolerance relative to a heat-sensitive line.The highest tolerance to heat stress was achieved by pyramiding all 11 of the identified candidate genes.Thus,the gene pyramiding scheme indicated in the present study provides a valuable direction for rice breeding.The results of our study will contribute to the development of elite rice accessions with increased tolerance to heat stress.

Fig.5.Analysis of group pyramiding combination for evaluating heat tolerance.Ref represents control of the genotype.Different groups represent different QTL combinations.Group 1,HT1+HT2+HT6;group 2,HT1+HT2+HT3+HT6;group 3,HT1+HT2+HT3+HT4+HT6;group 4,HT1+HT2+HT3+HT4+HT5+HT6.HT refers to the QTL name shown in Table S4.The percentage at the top represents the number of accessions corresponding to the QTL in the group accounting for the total number of accessions.Significant differences between Ref and groups were calculated by Student’s t-test.*,P＜0.05;**,P＜0.01;***,P＜0.001;****,P＜0.0001.

5.Conclusions

The present study was an investigation of heat-tolerancerelated genes in rice at the seedling stage using an approach integrating GWAS,transcriptomics,and AS analysis.Eleven promising candidate genes were identified,and may prove to be useful targets for dissecting via genetic and biochemical analyses the molecular mechanisms that underlie heat stress.

CRediT authorship contribution statement

Lianguang Shang:Supervision.Qian Qian:Supervision.Guanjing Hu:Supervision.Yingxue Yang:Conceptualization,Data curation,Formal analysis,Writing-review & editing,Validation.Chao Zhang:Conceptualization,Data curation,Formal analysis,Writing-review & editing,Validation.De Zhu:Conceptualization,Data curation,Formal analysis,Writing-review & editing,Validation.Huiying He:Validation.Zhaoran Wei:Validation.Qiaoling Yuan:Validation.Xiaoxia Li:Validation.Xu Gao:Validation.Bin Zhang:Validation.Hongsheng Gao:Validation.Bo Wang:Validation.Shuaimin Cao:Validation.Tianyi Wang:Validation.Yuhua Li:Validation.Xiaoman Yu:Validation.Longbiao Guo:Validation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Agricultural Science and Technology Innovation Program,Shenzhen Science and Technology Program(KQTD2016113010482651),Projects Subsidized by Special Funds for Science Technology Innovation and Industrial Development of Shenzhen Dapeng New District(RC201901-05),Guangdong Basic and Applied Basic Research Foundation(2019A1515110557),the National Natural Science Foundation of China(31771887),and Genome-wide Association Study of High Nitrogen Use Efficiency Related Traits in Rice and Breeding Application(2020-KYYWF-0237).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.02.011.