RNAseq analysis reveals pathways and candidate genes associated with salinity tolerance in a spaceflight-induced wheat mutant

Salinity stress has become an increasing threat to food security worldwide and elucidation of the mechanism for salinity tolerance is of great significance. Induced mutation, especially spaceflight mutagenesis, is one important method for crop breeding. In this study, we show that a spaceflight-induced wheat mutant, named salinity tolerance 1 (st1), is a salinity-tolerant line. We report the characteristics of transcriptomic sequence variation induced by spaceflight, and show that mutations in genes associated with sodium ion transport may directly contribute to salinity tolerance in st1. Furthermore, GO and KEGG enrichment analysis of differentially expressed genes (DEGs) between salinity-treated st1 and wild type suggested that the homeostasis of oxidation-reduction process is important for salt tolerance in st1. Through KEGG pathway analysis, “Butanoate metabolism” was identified as a new pathway for salinity responses. Additionally, key genes for salinity tolerance, such as genes encoding arginine decarboxylase, polyamine oxidase, hormones-related, were not only salt-induced in st1 but also showed higher expression in salt-treated st1 compared with salt-treated WT, indicating that these genes may play important roles in salinity tolerance in st1. This study presents valuable genetic resources for studies on transcriptome variation caused by induced mutation and the identification of salt tolerance genes in crops.

ity-tolerant mutant st1 in a screen of a large number of induced wheat mutants grown in hydroponics with high salinity. The germination rate of st1 was significantly higher than that of WT when assayed on 250 mM NaCl (Fig. 1a,b). Moreover, st1 seedlings clearly grew better than WT when plants were exposed to high concentrations of NaCl (Fig. 1c). Mutant plants treated with 200 or 300 mM NaCl had significantly increased shoot weight compared to WT (Fig. 1d). With 200 mM NaCl treatment, the sodium concentration in the shoot was notably decreased in st1 compared to WT (Fig. 1e). Additionally, the malondialdehyde (MDA) content, which is typically held to reflect the intensity of damage to the plasma membrane under stress, was significantly lower in the leaves of st1 than in the leaves of WT (Fig. 1f). All of these results indicate that the spaceflight-induced wheat mutant st1 plants were more tolerant to salt stress than were WT plants.

RNAseq analysis of WT and mutants shoot under salinity stress. The transcript changes in
WT and the st1 mutant with or without salinity treatment were investigated using RNAseq analysis based on high-throughput sequencing. Sequencing of each genotype/treatment group was done independently twice, and each sample was pooled from five plants; thus, 8 cDNA libraries were constructed. The raw transcriptome datas generated from this study have been submitted to the Sequence Read Archive (SRA), National Centre for Biotechnology Information (NCBI) with the accession number SRR5277542. The number of raw reads from each library ranged from 43 million to 65 million. After removing poor quality reads, adaptor contamination, and low quality regions, between 42 to 63 million clean reads were obtained from each of the libraries. From each library >67% of reads mapped to the release-31 version of wheat reference genome, including >58% with unique mapping and <9% multiple mapping (Supplementary Table S1). Furthermore, >98% and >95% of the clean reads had quality scores at the Q20 or Q30 (an error probability for base calling of 1% or 0.1%) 37 level, respectively. The density of reads was equally distributed across all chromosomes of the wheat reference genome ( Supplementary  Fig. S1) and 70% of the genes in the wheat genome database were detected within these samples, suggesting the RNAseq datas are well covered the wheat transcriptional region. Additionally, there was a clear linear relationship between the gene expression levels of the two biological replicates for each genotype/treatment group (R 2 > 0.97, Supplementary Fig. S2).
Among the homozygous variants, 6211 SNPs and 414 InDels, distributed among 21 chromosomes, were detected (Table 1). This is consistent with the findings of a study of transcriptome sequence variation in a space-induced mutant of Kentucky bluegrass, where more SNPs than InDels were also identified 38 . Of the 6211 SNPs identified in our study, there were more transitions than transversions with a transitions: transversions (Ti: Tv) ratio of 2.0 ( Fig. 2a). High Ti: Tv values have also been reported for other types of mutagenesis like fast-neutron irradiation 39 and EMS mutation 40 . However, the distribution of transitions and transversions is roughly similar in gamma-ray-induced mutation 40 . Among the four possible types of transversions, the number of C/G to G/C transversions was the highest, accounting for 34.4% of point mutations, while the number of T/A to A/T transversions was the lowest, accounting for 16.0% (Fig. 2a). In contrast, it has been reported that EMS-, gamma-ray-, and fast-neutron-induced mutations exhibit lower numbers of C/G to G/C transversions compared with other types of transversions 39, 40 . We next looked at the SNPs located within genes. In six genes SNPs resulted in termination mutations, and 110 genes of SNPs resulted in missense mutations. Among the genes with termination and missense mutations, only seven genes had increased expression and four genes had decreased expression (fold change ≥2) in the st1 compared to WT (Supplementary Table S2). Furthermore, the enriched GO terms and KEGG pathways for genes with termination and missense mutations were determined (Supplementary Table S2 and Fig. 2b,c). REVIGO was used to summarize the enriched GO terms (http://revigo.irb.hr/). The top ten biological GO terms in biological process enriched in genes with termination and missense mutations, defined based on the lowest over-represented p values, were analyzed by REVIGO (Fig. 2b). This program removes redundant GO terms and the similarity between terms is reflected by semantic space 41 . Interestingly, two mutated genes are annotated with the enriched GO term "sodium ion transport" (Fig. 2b). It has been reported that the Na + /H + antiporter SOS1 is important for the exclusion of sodium in Arabidopsis root cells 17 , and the wheat sodium transporters HKT1; 5A 42 and HKT1; 5D 43 function in decreasing Na + concentration in the leaves and thus improve salt tolerance. These mutated genes associated with sodium ion transport offer important clues for the mechanism of salt tolerance in the st1 mutant plants, and their gene functions need to be further analyzed. KEGG pathway analysis indicated that some important pathways, such as "metabolic pathways", "carbon metabolism" and "plant hormone signal transduction", were enriched in mutated genes (Fig. 2c, Supplementary Table S2), suggesting that spaced-induced mutation of metabolism genes led to variation in metabolism and changes the characteristics of the mutant.
Analysis of differentially expressed genes (DEGs) was another way to identify genes that may be responsible for salinity tolerance in the st1 mutant. Therefore, DEGs were identified from the RNA-seq data based on the criteria of fold change ≥2 for the comparisons salinity (S)_st1 vs st1, S_WT vs WT, S_st1 vs S_WT and an FDR < 0.05 (Supplementary Table S3). Cluster analysis of DEGs suggested that the number of genes down-regulated by Student's t-tests were used to assess the significance of differences from WT plants, *P < 0.05, **P < 0.01. salinity in both WT and st1 was generally higher compared with other subclusters (Fig. 3). In the S_st1 vs st1 comparison group, 3560 DEGs, including 1230 with up-regulated expression and 2330 with down-regulated expression, were identified (Fig. 4a). And 2395 transcripts (962 up-regulated and 1433 down-regulated genes) showed significant changes in S_WT vs WT (Fig. 4b). In S_st1 vs S_WT, 1585 genes were up-regulated and 1262  genes were down-regulated (Fig. 4c). Venn diagrams indicate that a large proportion of DEGs (63.3%) in S_WT vs WT overlap with the S_st1 vs st1 DEGs. However, only a small proportion of DEGs (17.4%) in the S_st1 vs S_WT comparison overlap with the salinity responsive genes (i.e. S_st1 vs st1 or S_WT vs WT DEGs). Because st1 is more tolerant to salinity compared with WT, the salt responsive genes of st1 that are also differentially expressed in the S_st1 vs S_WT comparison (in total 369 genes, Fig. 4d) are more likely to play roles in the salinity tolerance of st1.
Functional categorization of salinity-treated WT and st1 mutant DEGs. The top ten GO terms in biological process category enriched in DEGs of the S_st1 vs S_WT group based on the lowest p values were analyzed by REVIGO. For the up-regulated genes, "oxidation-reduction process" was the only significantly enriched GO term (p < 0.05) (Fig. 5a), suggesting the importance of this process for salinity resistance in st1. Consistent with this finding, previous transcriptome profiling experiments have indicated that oxidation-reduction processes are related to salt tolerance in various plants 25,27,28 . In contrast, in the S_st1 vs S_WT comparison, enriched highly significant terms (p < 0.05) for down-regulated DEGs included "response to water", "response to inorganic substance", "response to abiotic stimulus", "lipoprotein metabolic process", "lipid transport", "lipid localization", "oxidoreduction coenzyme metabolic process", "NADPH regeneration" and "glycerophospholipid biosynthetic process" (Fig. 5b). The fact that genes related to these processes, especially those involved in response to water, inorganic substance and abiotic stimulus, were up-regulated in the salt-treated WT compared to the salt-treated mutant, indicates that the WT is more sensitive than st1 to high salt, which is consistent with its salt-sensitive phenotype ( Fig. 1). Sixty-six GO biological process terms were significantly enriched in the S_st1 vs st1 up-regulated genes, whereas only 11 GO terms were enriched in S_WT vs WT up-regulated genes (Supplementary Table S4). Furthermore, the number of up-regulated genes involved in "response to water stimulus", "response to abiotic stimulus" and "response to chemical stimulus" was generally higher for the S_WT vs WT comparison compared with the S_st1 vs st1 comparison (Supplementary Table S4), indirectly suggesting the greater sensitivity of WT to salinity compared with st1. It is possible that the up-regulation of multiple genes involved in various biological processes in st1 under salinity contributes to salinity tolerance. KEGG pathway enrichment analysis showed that, among the top 20 enriched pathways (based on corrected p-value) for the S_st1 vs S_WT up-regulated genes, the "Starch and sucrose metabolism", "Linoleic acid metabolism", "Galactose metabolism" and "alpha-Linolenic acid metabolism" pathways were significantly  (Fig. 5c). These pathways have been well-documented to play roles in salinity tolerance 3,44,45 . For the down-regulated genes in this comparison, no pathway was significantly enriched, but the "Biosynthesis of secondary metabolites" pathway had the largest number of genes (Supplementary Fig. S3 and Supplementary Table S5). The pathways involved in "Arginine and proline metabolism", "Alanine, aspartate and glutamate metabolism", "Linoleic acid metabolism", "Plant hormone signal transduction", "Citrate cycle (TCA cycle)", "Pyruvate metabolism", "Biosynthesis of secondary metabolites", and "Butanoate metabolism" were commonly significantly enriched (p < 0.05) in S_st1 vs st1 and S_WT vs WT up-regulated genes (Supplementary Table S5). This suggests that these pathways are active under salinity stress. Among these pathways, the "Butanoate metabolism" pathway has previously not been shown to function in salinity response. 4-amino-butanoate is one of the intermediate products of butanoate metabolism, and it is closely connected with the PA pathway, which is suggested to be involved in plant resistance to salinity stress [10][11][12][13] . Additionally, "Butanoate metabolism" is also closely associated with the salinity responsive "Alanine, aspartate and glutamate metabolism" and "TCA cycle" pathways 27 . Therefore, the expression of genes involved in butanoate metabolism was affected by salinity stress in st1 and WT.
Key genes responsible for salt tolerance in the st1 mutant. Since the salt responsive genes showing differential expression in the S_st1 vs S_WT comparison might be important for salinity tolerance in st1, the overlap between DEGs from the S_st1 vs st1 and S_st1 vs S_WT comparisons were analyzed in detail ( Table 2,  Supplementary Table S6), especially the up-regulated genes. We also confirmed the RNA-seq results for selected DEGs by qPCR (Fig. 6). Among the up-regulated genes, most genes were involved in arginine and proline metabolism, oxidoreductase, hormone metabolism, transcription factor and signaling regulation (Table 2). Additionally, genes involved in "Alanine, aspartate and glutamate metabolism", "S assimilation and metabolism", "Glutathione metabolism", "Flavonoid metabolism", "disease related" and "energy metabolism", were also up-regulated in both the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Supplementary Table S6). Generally, distinct genes were observed in the down-regulated group (Supplementary Table S6).
Among genes involved in arginine and proline metabolism, transcript levels of five genes encoding polyamine oxidase (PAO) family members were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons ( Table 2). Similar expression patterns for these genes were observed in qPCR analysis (Fig. 6, Q1-Q4). The mRNA of the gene encoding arginine decarboxylase (ADC), an enzyme that catalyzes the synthesis of Put, was also detected as having increased expression in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2). PAO and ADC genes are critical for PA metabolism, and PAs play important roles in plant resistance to environmental stresses 11,14,[46][47][48] . For example, it has been reported that a cotton PAO is required for the mediation of spermine and camalexin signaling to resist a destructive fungal pathogen 16 . Additionally, PAO has also been reported to contribute to improving drought stress tolerance in grapevine 15 . It has been shown that Put plays a positive role in salt tolerance by attenuating oxidative damage in soybean roots 12 . Furthermore, by modulating nutrient acquisition and the PA pool, inoculation with arbuscular mycorrhizal fungi improves adaptation to saline soils 49 . Several studies have suggested that the modulation of PA catabolism is important for salinity tolerance [50][51][52][53] . In this study, the expression levels of key genes in PA metabolism were increased in the S_st1 vs st1 and the S_st1 vs S_WT comparisons, suggesting that the modulation of PAs is important for salt tolerance in st1.    Table 2. The RNA-seq values represent the ratio of the expression level in st1 to the expression level in WT. The expression level for each sample is the average of two replicates. For qPCR, the data are means ± SD from three independent replicates, and different letters indicate significant differences between samples at P < 0.05 based on SAS statistic analysis.
Oxidation-reduction processes are also critical for salinity tolerance in plants. In accordance with the enrichment of the GO term "oxidation-reduction process" in up-regulated genes of the S_st1 vs S_WT comparison, several oxidoreductase genes, including genes encoding cytochrome P450 monooxygenases (CYP), ascorbate oxidoreductase and glutaredoxin, had increased expression in the S_st1 vs st1 and the S_st1 vs S_WT comparisons. Specifically, eleven genes putatively encoding CYP were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2). Similar expression patterns for these genes were observed in qPCR analysis. (Fig. 6, Q5-Q9). CYPs catalyze the biooxidation of various substrates through activation of molecular oxygen and play important roles in metabolic processes and stress responses 54 . It has been suggested that increased expression of CYP94 genes in rice plants alleviates jasmonate responses and enhances salt tolerance 55 . Mutation of CYP709B3 in Arabidopsis conferred sensitivity to ABA and low tolerance to salt stress, indicating that CYP709B3 plays a role in ABA responses and salinity tolerance 56 . In addition, two genes encoding L-ascorbate oxidase-like and ascorbate-dependent oxidoreductase, which are involved in ascorbate recycling, had increased expression in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2). Ascorbate plays important roles in protecting plants against various stresses by regulating cellular H 2 O 2 levels 57 . Two glutaredoxin-C1 genes, the glutathione-dependent enzyme involved in reducing disulfide bridges, were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons ( Table 2). Overexpression of tomato glutaredoxin results in tolerance to oxidative, drought, and salt stresses, suggesting that glutaredoxin plays a crucial role in regulating abiotic tolerance 58 . On the whole, the oxidoreductases are involved in scavenging reactive oxygen species and maintaining oxidation-reduction homeostasis 3 and thus protect the salt-treated plants from damage.
Hormones, such as jasmonate, ethylene, abscisic acid and cytokinins, are also known to be involved in the regulation of salinity tolerance in plants 59 . Increased expression of allene oxide cyclase in wheat and Arabidopsis resulted in increased jasmonate content and in improved resistance to salinity 60 . In our study, allene oxide synthase, a key gene involved in the biosynthesis of jasmonic acid, was highly expressed in salt-treated st1 ( Table 2). qPCR analysis also showed that this gene had significantly higher expression in salt-treated st1 than in untreated st1 or in salt-treated WT (Fig. 6, Q10). Two 12-oxophytodienoate reductase genes and four lipoxygenase family genes, which also participate in the jasmonate synthesis pathway, and had increased expression levels in the S_st1 vs st1 and the S_st1 vs S_WT comparisons ( Table 2). The high expression of these jasmonate biosynthetic genes may contribute to salinity tolerance in st1 under high salinity. Makhlouf et al. reported that an ethylene response factor gene is probably associated with salt tolerance in wheat 61 ; ethylene also plays roles in salt responses in plants 59,62,63 . Five genes involved in ethylene biosynthesis or regulation, were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2), and these expression patterns were confirmed by qPCR (Fig. 6, Q11-Q13). Changes in the expression of ethylene-related genes observed in the S_st1 vs st1 and the S_st1 vs S_WT comparisons indicate that ethylene may also contribute to salinity tolerance in st1. Interestingly, the expression of putative ethylene-responsive transcription factor (Q12) and TaAP2-D (Q13) genes was barely detectable in WT and salinity-treated WT but highly expressed in st1 (Fig. 6). The two novel genes may be involved in salinity tolerance in the st1 mutant plants. In addition, a transcription factor involved in abscisic acid signal transduction and two genes related to cytokinin showed increased expression levels in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2 and Fig. 6, Q14-Q15). Taken together, the high expression of genes related to jasmonate, ethylene, abscisic acid, and cytokinins in the salt-treated st1 mutant suggest that these hormones may explain the improved tolerance to salt stress observed in this mutant.
More than ten transcription factor and signal transduction genes were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons (Table 2). Similar expression patterns were observed for the selected genes in qPCR analysis. (Fig. 6, Q16-Q18). It has been reported that high levels of expression of the Na + -induced TaMYB73 gene in Arabidopsis results in increased tolerance to salt stress 64 . In our study, two salt-induced genes encoding MYB-related proteins were more highly expressed in salt-treated st1 compared with salt-treated WT plants ( Table 2), suggesting that this gene may function in tolerance to salt stress in st1. Three genes encoding CBL-interacting protein kinases (CIPK) were up-regulated in the S_st1 vs st1 and the S_st1 vs S_WT comparisons ( Table 2). The CBL/CIPK network has been shown to regulate the Na + efflux transporter SOS1 in Arabidopsis 18,23,65 . Several CIPK genes from various plants have been reported to function in salt tolerance [21][22][23] , and it is reasonable to speculate that the high levels of expression of CIPK genes in st1 promotes tolerance to salt stress. Additionally, some transcription factors, including a MADS-box gene (Fig. 6, Q16), zinc finger proteins (Fig. 6, Q17-18), and a homeobox-leucine zipper protein, were newly identified as putative salt tolerance genes in this study (Table 2 and Fig. 6).
In conclusion, comparative analysis of the shoot transcriptomes of WT and the space-induced wheat mutant st1 not only indicates that multiple genes and pathways are related to plant responses to salinity stress, but also provides a catalog of transcriptomic sequence variation induced by spaceflight. The data generated in this study will be a valuable resource for the identification of key genes related to salinity tolerance in wheat and also for studies of the sequence variation resulting from induced mutation.

Materials and Methods
Plant materials. The core parent Jimai 20 from Shandong province was used as the WT, and mutation was induced by spaceflight on the Shijian-8 satellite in 2006. The satellite Shijian-8 was in orbit for 15 days; its mission was specific to seed-breeding programs. The average irradiation dose for the seeds in the cabin of the satellite was about 2.894 mGy, and the mean gravity was 1. 3 × 10 −3 g 66 . The M 6 generation was used for salt-tolerance screening. Salt-tolerance screening was conducted in 1.5% NaCl solution for several days, and then the growth of the mutants was monitored. One mutant, named st1, with a high germination rate and superior growth characteristics was selected for further analysis.
Scientific RepoRts | 7: 2731 | DOI:10.1038/s41598-017-03024-0 Salt treatment and physiological analysis. After germination for 2-3 days in deionized water, the WT and mutant were grown in half concentration nutrient solution without NaCl, or with 100, 200, 300 mM NaCl for five days under 200-300 µmol m −2 s −1 light at 21 °C in a growth chamber. The composition of nutrient solution was according to Guo et al. 67 . The shoot weights from five independent replicates of WT and st1 plants treated with different concentrations of NaCl were measured. The Na + content in the shoot was measured from at least three replicates according to a previously-published method 2 . The malondialdehyde (MDA) content was measured as described previously 68 . For germination analysis, the WT and st1 seeds were directly cultured in water containing 250 mM NaCl for several days. The germination rate was calculated as the number of seeds germinated divided by the total number of seed.
RNA extraction and library construction. The shoots from five plants for each treatment were pooled as one sample replicate. Two replicates for each treatment were prepared to isolate total RNA. Total RNA was isolated by using an RNeasy Plant Mini Kit (Qiagen) according to the product instructions. DNase I (Takara) and an RNA purification kit (Tiangen) were used to eliminate DNA contamination. The purity, concentration and integrity of RNA was assessed by using a NanoPhotometer ® spectrophotometer (IMPLEN), a Qubit ® RNA Assay Kit with a Qubit ® 2.0 Flurometer (Life Technologies), and an RNA Nano 6000 Assay Kit with the Bioanalyzer 2100 system (Agilent Technologies). Sequencing libraries were generated using an NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB) following the manufacturer's recommendations. Briefly, mRNA was purified by using poly-T oligo-attached magnetic beads. After fragmentation, the cDNA was synthesized, and NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. PCR products were purified (AMPure XP system), and library quality was assessed on an Agilent Bioanalyzer 2100 system. Finally, eight libraries were successfully constructed.
Transcriptome analysis. The libraries were deep sequenced with the Illumina Hiseq platform. For sequence quality control, clean reads were obtained by removing reads containing adapters, reads containing poly-N homopolymers and other low quality reads from the raw data. The Q20, Q30, and GC content in the clean reads were calculated to further assess quality. Release-31 version of the wheat reference genome and gene model annotation files were downloaded from the genome website (ftp://ftp.ensemblgenomes.org/pub/plants/release-31/ fasta/triticum_aestivum/dna/). An index of the reference genome was then built using Bowtie v2.2.3, and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12. For SNP analysis, Picard-tools v1.96 and samtools v0.1.18 were used to sort, mark duplicated reads, and reorder the bam alignment results of each sample. GATK2 (v3.2) software was used to perform SNP calling.

Quantification of gene expression levels and DEG analysis.
The number of reads mapped to each gene were counted using HTSeq v0.6.1. The expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) of each gene was then calculated based on the length of the gene and the read count mapped to this gene. FPKM, reflecting both the effect of sequencing depth and gene length for the read count, is a commonly used method for estimating gene expression levels 69 . DEG analysis was performed using the DESeq R package (1.18.0). DESeq provides statistical routines for determining DEGs using a model based on a negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Genes with an adjusted P-value determined to be <0.05 (FDR < 0.05) by DESeq and that had a fold change value ≥2 (|Log2 fold change| ≥ 1) between two groups were considered to be differentially expressed.
GO and KEGG enrichment analysis. For gene ontology (GO) mapping, the GO terms of DEGs associated with homologies (GO; http://www.geneontology.org) were extracted. GO enrichment analysis of DEGs was implemented using the GOseq R package, in which gene length bias was corrected. GO terms with corrected P values less than 0.05 were considered to be significantly enriched. REVIGO was used for analysis of the enriched GO terms (http://revigo.irb.hr/); this program removes redundant GO terms and attempts to reflect the similarity of given terms by semantic space 34 . The ten biological process category GO terms with the lowest p values for enrichment in both the up-or down-regulated genes in the S_st1 vs S_WT comparisons were analyzed by REVIGO. KEGG pathways for the DEGs were retrieved (http://www.genome.jp/kegg/), and KOBAS software was used to test the statistical significance of the enrichment of DEGs in KEGG pathways.
Quantitative real-time PCR. The shoots from salt-treated and control wheat plants were sampled. Five plants were mixed as one replicate, and three independent biological replicates with three technical replicates of each biological replicate were analyzed. Total RNA was isolated using an RNeasy Plant Mini Kit (Qiagen). DNA contamination was removed using DNase I (Takara) and an RNA purification kit (Tiangen). First-strand cDNA was synthesized using an iScript cDNA synthesis kit (Bio-Rad). Quantitative real-time PCR was conducted by using SsoFast EvaGreen Supermix Kit (Bio-Rad) on a CFX 96 Real-Time System (Bio-Rad). The primers used for quantitative real-time PCR are detailed in Supplementary Table S7.