Mechanisms of wheat (Triticum aestivum) grain storage proteins in response to nitrogen application and its impacts on processing quality

Basis for the effects of nitrogen (N) on wheat grain storage proteins (GSPs) and on the establishment of processing quality are far from clear. The response of GSPs and processing quality parameters to four N levels of four common wheat cultivars were investigated at two sites over two growing seasons. Except gluten index (GI), processing quality parameters as well as GSPs quantities were remarkably improved by increasing N level. N level explained 4.2~59.2% and 10.4~80.0% variability in GSPs fractions and processing quality parameters, respectively. The amount of N remobilized from vegetative organs except spike was significantly increased when enhancing N application. GSPs fractions and processing quality parameters except GI were only highly and positively correlated with the amount of N remobilized from stem with sheath. N reassimilation in grain was remarkably strengthened by the elevated activity and expression level of glutamine synthetase. Transcriptome analysis showed the molecular mechanism of seeds in response to N levels during 10~35 days post anthesis. Collectively, we provided comprehensive understanding of N-responding mechanisms with respect to wheat processing quality from N source to GSPs biosynthesis at the agronomic, physiological and molecular levels, and screened candidate genes for quality breeding.

Wheat is one of the major cereal crops in the world, and provides an important source of vegetable protein for human beings. Compared to rice and corn, wheat is the source of numerous foods depending on the unique properties of its grain storage proteins (GSPs) 1 . The quantity and composition of GSPs are associated with processing and nutritional quality of wheat, and thus affect wheat price in essential market.
GSPs in wheat, mainly referring to gliadins and glutenins, together accounting for 60~80% of the total seed proteins, are important determinants of processing quality, as they are responsible for the elasticity and extensibility of dough that determine the processing qualities of various end-products 2 . Gliadins are monomeric, and have been divided into three classes, i.e. α/β-, γ-, and ω-gliadins. Glutenins are polymeric fractions composed of high-molecular-weight (HMW-GSs) and low-molecular-weight glutenin subunits (LMW-GSs) linked with intermolecular disulfide bonds 3 .
In wheat, 60~95% of the grain nitrogen (N), which is principally used for GSPs synthesis, derives from the remobilization of N stored in vegetative organs before anthesis, particularly in leaf blades, stem and sheaths 4,5 . During senescence after anthesis, proteinaceous components of vegetative organ cells are degraded into amino acids, amides and ammonium 6 . The major part of ammonium is incorporated into amino acids for export by glutamine synthetase (GS) and glutamate synthase (GOGAT) cycle 7 . Subsequently, amino acids as the major N-transport compound are exported via the phloem into developing grains to synthesize GSPs. N remobilization is closely correlated to the synthesis of GSPs 8 , while there is little imformation on the relationship between N remobilization and GSPs synthesis in wheat. In addition, a large number of studies focus on GS enzyme in leaf due to its major importance for the re-assimilation of ammonium into exported amino acids during senescence 7,9,10 , whereas its function in developing seeds is usually ignored.
N is one of the curial macronutrients for plant growth. It is generally accepted that N application, as an important agronomic input, plays a vital role in the well balance of yield and processing quality in wheat production, by delaying vegetative senescence and increasing the amount of N reserve at anthesis 11,12 . N remobilization is greatly strengthened by the elevated GS activity under abundant N level 13,14 . Adequate N supply exerts a significant increase in protein content and dough quality 15,16 . To date, researchers have investigated the response of protein accumulation during grain development of wheat to varied N by transcriptome, proteomics and metabolite profiling methods 15,17,18 . However, very little information covers the entire grain filling stage and exclusively illustrates physiological, biochemical and molecular changes in GSPs biosynthesis under sufficient N supply. Our knowledge of the underlying basis for the establishment of processing quality after N application remains fragmentary.
Here, we performed field trials with four common wheat cultivars under four N treatments at two sites over two wheat growing seasons. Transcriptome profiling in developing grain was combined with the phenotypic and physiological characters, such as the accumulation of GSPs fractions, processing quality parameters, the status of N remobilization in wheat plants, and dynamic enzyme activities related to N reassimilation. This research would deepen our understanding on the mechanisms of the establishment of wheat processing quality in response to N application. In addition, the identified genes would be important candidates for wheat quality breeding in the future.

Results
Processing quality parameters. Processing quality parameters of the four cultivars shown a similar trend after increasing N application at two sites in two consecutive growing seasons (Table 1). GPC, SV, GC, GI, DT and ST were significantly influenced by N level and genotype at both sites in two growing seasons. The percentages of variability of GPC, SV, GC, DT and ST explained by N level accounted for 57.8~80.0%, 23.3~72.6%, 37.0~63.7%, 16.2~41.9% and 10.4~31.6%, respectively, while that of GI was only 3.8~16.0%, indicating that GI was mainly controlled by genotype. Only DT and ST showed significant N level × genotype interaction at both sites in two growing seasons, and the percentages of variability of DT and ST explained by N level × genotype interaction were much lower than that explained by N level or genotype. To summarize, GPC, SV, GC, DT and ST of four cultivars were statistically improved when increasing N level, whereas GI was decreased by 5.3~14.0%, especially in 2015-2016 growing season.

Content and composition of glutenin and gliadin fractions in flour.
The content and composition of glutenin and gliadin fractions were significantly affected by N level and genotype but not the N level × genotype interaction (Table 2). Genotype contributed more variance to the content and composition of glutenin and gliadin fractions than N level, indicating the diversity of cultivars used in this experiment. The application of N significantly promoted the accumulation of glutenin and gliadin fractions, and the percentage was highest for HMW-GSs from N 0 to N 225 treatment (86.6% at Chongzhou and 50.7% at Renshou), followed by LMW-GSs, α/β-gliadin, γ-gliadin and ω-gliadin. Notably, Glu/Gli was decreased when increasing N level.
RP-HPLC data were analyzed by principal component analysis. The first four principal components (PCs) accounted for 94.1% of the variation in the distance matrix formed from the similarities between samples (Supplementary Table S1). PC1 and PC2 covered 43.7% and 32.7% of the total variation, respectively, indicating that the first two PCs could explain the variation of protein fractions under differed N levels. PC1 axis mainly depended on total gliadins, HMW-GSs, γ-gliadin, ω-gliadin and α/β-gliadin. LMW-GSs, ω-gliadin%, and total glutenins were important for PC2 axis. The four cultivars were gathered in two groups. Meanwhile, the four N levels were clearly separated from each other. The higher N levels had a trend of migration from negative axis to positive axis, and the migrated distance was positively associated with N level (Fig. 1). N remobilization from vegetative organs. The difference of N concentration (DNC) between maturity and anthesis was highest in flag leaf blade, and then in the lower leaf blades, followed in spike and stem with sheath (Fig. 2). Significant differences were identified in the amount of N remobilized from vegetative organs, among four N levels for four cultivars. The amount of remobilized N at Renshou was higher than that at Chongzhou. DNC in flag leaf blade, the lower leaf blades and stem with sheath were significantly increased with the elevated N application at both sites, while DNC in spike differed between both sites, and even decreased at Renshou.
Correlations among processing quality parameters, content and composition of GSPs, and N remobilization. Significant and varied correlations were broadly observed between content and composition of GSPs and processing quality parameters (Table 3). For processing quality parameters, GI and GC had a stronger correlation with the composition than the content of GSPs. For the composition of glutenin and gliadin fractions, γ-gliadin% and Glu/Gli were significantly associated with all processing quality parameters.
DNC in stem with sheath displayed the highest correlation with processing quality parameters and with the content and composition of GSPs (Table 4). DNC in stem with sheath was significantly and positively correlated with GPC, LMW-GSs, ω-gliadin and GC, while significantly and negatively correlated with GI (Table 4). It suggested that N remobilized from stem with sheath was critical for the accumulation of LMW-GSs and ω-gliadin, resulting in increased GPC in grain and altered processing quality. In addition, GI was the only one which had high correlation with DNC in vegetative organ in all above processing quality parameters. Extremely significant and positive correlation was found between GI and DNC in spike, whereas negative correlation was detected between GI and DNC in leaf, stem and sheath.
Enzyme activity of GS. GS activities in grains exhibited rapidly decreasing trends during 10~35 DAA, especially during 10~20 DAA (Fig. 3). GS activities of SM482 and MM51 under N 150 and N 225 treatments were significantly higher than those under N 0 and N 75 treatments at most time points. It suggested that increased N application was conducive to improve GS activities in grains, which was supported by RNA-seq and qRT-PCR (Fig. 4). Transcriptome analysis (N 0 VS N 225 ) revealed that three GS transcripts, one encoding GSr1 and two encoding GSr2 isoforms, were always up-regulated (Log2 > 0.7) with high expression level during grain filling stage.
KEGG and GO enrichment analysis. There were 2293 DEGs (N 0 VS N 225 ) screened during 10~35 DAA (Supplementary Table S2), with 158 DEGs being detected at more than one time point (Supplementary  Table S2B). To validate the expression profiles obtained by RNA-seq, qRT-PCR was performed to measure the expression of 10 DEGs (with high fold changes, with high expression level (FPKM > 100), and in starch or protein biosynthesis) during 10~35 DAA. As a result, a high correlation (R 2 = 0.818) was identified between qRT-PCR and RNA-seq data, confirming the validity of RNA-seq data ( Supplementary Fig. S1).
DEGs involved in amino acid metabolism. DEGs related to grain amino acid metabolism were mainly observed during 25~35 DAA. Genes associated with alanine, aspartate and glutamate metabolism exhibited the most dominant changes under N 225 treatment (Supplementary Table S5). Alanine and glutamate metabolism were strengthened by the up-regulation of alanine aminotransferase 2 and GSr1, GSr2 and GSr3, respectively, under N 225 treatment. However, aspartate metabolism was diminished by the three down-regulated DEGs encoding asparaginase. Five DEGs involved in phenylalanine, tyrosine and tryptophan biosynthesis were upregulated as well. No DEG involved in the metabolism of other amino acids was identified.
DEGs encoding GSPs. The expression of GSPs genes exhibited an increasing trend from 10 DAA, and reached the maximum during 30~35 DAA (Table 5). Ten DEGs encoding globulins were detected, six of which   Table 3. Correlation coefficients between content and composition of glutenin and gliadin fractions and wheat processing quality parameters. Abbreviations are used as in Tables 1 and 2. "*" and "**" represent significance at P < 0.05 and P < 0.01, respectively. were down regulated (FDR < 0.05). Twelve DEGs encoding glutenins and gliadins were detected, and most of them were up-regulated by 0.26~1.91 (log2; FDR < 0.05), especially the highly expressed DEGs encoding y-type HMW-GSs, ω-gliadin and α-gliadin. These results were basically consistent with the result of RP-HPLC analysis.
DEGs encoding transcription factor. Twenty six DEGs encoding transcription factor were found with at least one fold change (log2, FDR < 0.05), including 11 up-regulated and 15 down-regulated (Table 6). Four, two and two of up-regulated DEGs were responsible for encoding NAC, MYB, MADS-box transcription factors, respectively. The four DEGs encoding NAC appeared during 25~35 DAA. Three and four of down-regulated DEGs were in charge of encoding bHLH and WRKY transcription factors, respectively, and their expression were remarkably altered during 10~25 and 15~35 DAA, respectively.

Discussion
N application is considered as the most important agronomic input due to its contribution to processing quality of wheat, whereas the mechanisms on the contribution of N remain unclear. In the current study, a combination of physiological and transcriptomic analysis was used to identify the basis to partition N into GSPs in response to increased N inputs during grain development, and the putative mechanism based on our research data was summarized in Fig. 6.  Table 4. Correlation coefficients for relation between DNC (difference of N concentration between maturity and anthesis) in vegetative tissues and RP-HPLC data & processing quality parameters. Abbreviations were used as Tables 1 and 2. "*" represents significance at P < 0.05. N application has a beneficial effect on GSPs and its fraction in flour 15,16 . Similar result was observed in the present study. As for GSPs fraction, HMW-GSs and ω-gliadin, especially ω-gliadin, greatly responded to increased N by increasing their accumulation in grain 19,20 . In our experimental conditions, the elevated percentage was the lowest for ω-gliadin, and that was the highest for HMW-GSs at both sites. Meanwhile, principal components analysis indicated that the contribution of HMW-GSs was the highest among GSPs fractions on PC1 axis, consistent with previous result 17 . Taken together, HMW-GSs accounted for the largest effect of N application on GSPs factions, which was possibly responsible for the dominant increase in the ratio of HMW-GSs to LMW-GSs, thereby resulting in pronounced variation in wheat processing quality.
It is well accepted that a large number of grain N is derived from the remobilization of pre-anthesis absorbed N in vegetative organ. Although approximately half of the mobilized N originated from leaf blades, comparable difference happens in elsewhere, in stem or sheath 8,13 . In our experimental conditions, GSPs fractions and processing quality showed strong correlations with the mobilized N from stem with sheath, rather leaf blades or spike, suggesting that the altered GSPs accumulation and processing quality after N application were largely explained by the mobilized N from stem and sheath. Therefore, it would be valuable to pay more attention to the mobilized N from stem and sheath in future programs for wheat breeding. Only GI among these processing quality parameters exhibited clear relationships with the mobilized N from all vegetative organs (Table 4). GI was positively associated with the mobilized N from spike, while negatively associated with that from leaves, stem and sheath. Hence, it was reasonable to speculate that GI was a combined result of the mobilized N from all vegetative organs. The capacity of N remobilization can be inhibited by excess N 21 . Here, only N remobilization from spike  was inhibited by increased N application to some extent (Fig. 2), thereby breaking the homeostasis of composition in GSPs fractions, and eventually leading to the decrease of GI. This observation indicates that the capacity of N remobilization in spike prior to other vegetative organs is repressed by N.
GS, consisting of the cytosolic (GS1) and chloroplastic (GS2), plays an essential role in plant nitrogen metabolism for its responsibility to assimilate ammonium and transform it into glutamine at the first step in GS and GOGAT cycle 22,23 . GS1, including GS1 (a, b and c), GSe (1 and 2), and GSr (1 and 2), contributed more than GS2 to the reassimilation of N 10 . Numerous studies focused on the role of GS enzyme in leaf, whereas its role in grain remained unclear. QTL analysis indicated that GS1 was associated with grain size 24 , and GS2 was co-localized with the locus for GPC 22 . Furthermore, GPC was positively relevant to the activity of GS in leaf 10 , and GS activity in leaf was significantly improved by N application 13,14 . GS activity in developing grain and GPC were positively correlated with N level as well in this research, which was supported by RNA-Seq and qRT-PCR data (Table S5 and Fig. 4). We demonstrated that ammonium reassimilation remained active in grain, thereby providing additional supply of amino acids after N application.
GSPs synthesis needs abundant supply of amino acids, which are principally transported from vegetable organ via phloem and tightly associated with signal transduction between nutrient availability and GSPs accumulation 25 . Increasing N level induced coordinated up-regulation of genes involved in the metabolism of alanine (Ala), glutamate (Glu), phenylalanine (Phe), tyrosine (Tyr) and tryptophan (Trp), but it down-regulated some genes responsible for aspartate metabolism in grain during filling stage. Importantly, DEGs encoding Tyrosyl-tRNA synthetase, Alanyl-tRNA synthetase, Valyl-tRNA synthetase, Glutaminyl-tRNA synthetase exhibited up-regulation during early and late filling stage (Supplementary Table S6). These consistent evidences manifested that the transport and synthesis of Glu, Val, Ala and Tyr responded to N application, resulting in increased accumulation of GSPs. Meanwhile, it was worthy to note that the diminished DEGs encoding asparaginase at high N level might lead to the decrease of asparagine that generate the carcinogen acrylamide during the baking process in mature grain 26 . It suggested that nitrogen application was beneficial to produce more healthy wheat flour, and the adjustment in amino acid metabolism was necessary to adapt GSPs accumulation in response to N application.
GSPs are synthesized on the rough endoplasmic reticulum (rER) (Fig. 6). Firstly, mRNA specially binds to ribosomes on the rERs, thereby forming polypeptide chains, which are the precursors of GSPs. Subsequently, these precursors of GSPs are transferred into the lumen of the rERs where lots of post-translational processes happened, such as cleavage of signal peptide, formation of intra-and inter-chain disulphide bonds, folding, assembly and aggregation of polypeptides. Finally, gliadins and glutenins from rERs are deposited into protein bodies (PBs) in the vacuole of endosperm cell via the Golgi apparatus and autophagy, respectively 27 . These processes are associated with rER-localised chaperones and foldases, which are divided into three groups, protein disulphide isomerases (PDI), peptidyl-prolyl cis-trans isomerases (PPIase), and the binding protein (BiP)/heat shock protein (Hsp70) chaperone [27][28][29] . According to KO enrichment analysis, we found that protein processing at rER and ribosomes were very active. PDI and HSP70 genes were down-regulated by high N supply during   Table S6). Meanwhile, the genes coding for the 15~30 kDa small heat shock proteins (sHSPs) and the molecular chaperones that block the aggregation of unfolded proteins 30 , expressed as PDI and HSP70. Therefore, it was conceivable that rER responded to N starvation by up-regulation of PDI, BiP and sHSPs family during mid-filling stage, thereby enhancing the processing and accumulation of GSPs, and that might lead to high N use efficiency. On the other hand, PPIase gene family including PPIase Pin1 and FK506-binding protein 2 (FKBP) 31 , showed up-regulation during early and late filling stage. Moreover, BAG2, and DnaJ8, encoding co-chaperone that assist HSP70 to carry out chaperone activity 32,33 , as well as calreticulin that possesses molecular chaperone activity in protein folding 34 , were up-regulated during late filling stage. Although PPIases and co-chaperone of HSP70 were often ignored during protein folding in previous studies, they together with calreticulin were crucial for facilitating GSPs synthesis to cope with excess amino acids at rER under high N condition, particularly during early and late filling stage. Ribosomal proteins (RPs) are the major component for ribosome and play crucial roles in protein synthesis. Our study showed that most of RPs genes were repressed at high or low N level in certain period, especially RPs 40s S8 and 60S L12, which exclusively expressed at low N level during 10~20 DAA and at high N level during 25~35 DAA. Consequently, we speculated that ribosomal structure located at rER might change with N supply in grains, and RPs 40s and 60S were good candidate genes for respecting N status. In this study, the expression level of DEGs encoding GSPs was basically consistent with its accumulation except γ-gliadin, which has minor effect on properties of wheat flour 35 . The expression of GSPs genes are modulated by specific transcription factors (TFs) in response to N supply. SPA (Storage protein activator) and DOF (DNA binding with one finger) was a generally accepted TF for GSPs 27 , and they was up-regulated under sufficient N condition during early and late filling stage. In addition, several remarkably upregulated TFs (NAM, MADS-box 29 and 58, MYB 44 and 86) during mid and late filling stage, were also possibly involved in the transcriptional regulation of GSPs genes subjected to high N treatment. NAC, one of super TF family in plant, was highlighted for direct relation with leaf senescence and N allocation to grain during grain filling 36,37 , due to consisting of A to E domains for DNA-binding and protein-protein interactions 38 . However, little is known about its role in developing grain. MADS-box genes, with a highly conserved 180-bp-long motif, were emphasized curial roles in flower development and organ differentiation 39 . A MADS-box gene SlFYFL was proved to delay senescence, fruit ripening and abscission in tomato 40 . MYB TFs, with the presence of 1~4 or more imperfect MYB repeats domain located near the N-terminus, were vital regulators and modulate diverse biological processes, including plant response and tolerance to the nutrient deficiency 41 . Collectively, these TFs could serve as promising candidate genes for wheat quality breeding and deserve further research.

Conclusions
The establishment of wheat processing quality was significantly affected by N fertilizer. Sufficient N supply promoted the remobilization of pre-anthesis absorbed N in vegetative organ and even N reassimilation governed by GS enzyme in grains, resulting in increased accumulation of GSPs. Meanwhile, the high-throughput RNA-seq analysis allowed us to investigate an overall survey of genes and processes potentially related to the response of GSPs to stimuli induced by high N level, which were mainly involved in adaptation in amino acid metabolism, protein processing at ER, transcriptional regulation for GSPs genes. Therefore, these comprehensive observations deepen our understanding in the mechanisms of GSPs in response to N nutrient and optimization of N strategy for wheat cultivar with superior quality.

Materials and Methods
Wheat cultivars and field experiments.  Supplementary Table S7. The soils at Chongzhou and Renshou experimental site were dark brown clay and purple sandy loam, respectively. The basic soil properties in 0~20 cm layer were given in Supplementary  Table S8. The four cultivars were evaluated at four levels of N supplies (urea), i.e. 0, 75, 150 and 225 kg ha −1 (referred to as N 0 , N 75 , N 150 and N 225 , respectively). The field experiment was laid out as a split-block design with three replicates. Each individual plot size was 2 × 5 m at Chongzhou and 1.8 × 4.5 m at Renshou. Seeds were sown at 250 seeds m −2 with 20 cm row spaces in both field spots. The application of N as urea was spilt into 60% before sowing and 40% at stem elongation with immediate irrigation. A prophylactic programme for disease, weed, and pest management was used to maintain undisturbed healthy crop growth. Grains were harvested during mid-May at Chongzhou and early May at Renshou.
The flowering spikes, blooming on the same day, were marked with colored line, and the grains from middle ear at 10, 15, 20, 25, 30 and 35 days after anthesis (DAA) were collected at 2:00-3:00 pm. The collected samples were immediately frozen in liquid N 2 and then stored at −80 °C.
Evaluation of processing quality. Total plot was harvested and threshed using a mini-Vogel machine for recording yield. Two kilos of seeds were randomly picked from each grain sample to evaluate grain processing quality. Grain samples were stored for three months before milling. Milling was done on a Brabender Quadromat Juniors ® (Brabender GmbH & Co. KG, Germany). Samples were conditioned to 14.0% moisture before milling. The grain protein content (GPC), zeleny sedimentation value (SV), wet gluten content (GC) and gluten index (GI) were determined following GB/T 17320-2013, by an automatic azotometer (Kjelec 8400; FOSS, Denmark), a zeleny analysis system (CAU-B, China) and a glutomatic 2200 system (Perten, Sweden), respectively. Dough rheological properties were tested by Farinograph ® -AT (Brabender GmbH & Co. KG, Germany) with Mixer*50.

Determination of the contents of gliadin and glutenin fractions.
For each sample, 45 mg flour was used to extract gliadins and glutenins. The gliadin and glutenin fractions were analysed according to previous method 42,43 with minor modifications, respectively. The gliadin and glutenin fractions were both filtered through 0.45 µm nylon filter (Teknokroma, Barcelona, Spain) before performing Reverse-phase high-performance liquid chromatography (RP-HPLC) analysis. Gliadin and glutenin extracts were applied to a ZORBAX 300SB-C18 reverse phase analytical column (4.6 × 250 mm, 5-Micorn; Agilent Technologies, Santa Clara, CA) using a 1100 Series Quaternary LC System liquid chromatograph (Agilent Technologies, Santa Clara, CA) with a DAD UV-V detector. Absorbance was monitored with the DAD UV-V module at 210 nm. The major analytical parameters were set as 60 °C for column temperature, 1.00 ml/min for flow rate, 20 µl for sample volume, eluting gradient and the variable concentrations of acetonitrile with 0.06% trifluoroacetic acid for gliadins and glutenins gradually growing from 21% to 48% (v/v) in 55 min and from 25% to 48% (v/v) in 50 min, respectively. Column washing time for gliadins and glutenins between two adjacent samples were 15 and 10 min, respectively. The integration procedure was handled automatically by the Agilent Technologies Chemistry Station software with minor manual adjustment. The amounts of HMW-GSs, LMW-GSs and gliadins were estimated by integrating the relevant peak areas in the chromatograms.
Measurement of the ratio of glutenins to gliadins. Albumins, globulins, gliadins and glutenins were sequentially extracted using distilled water, 10% NaCl, 70% ethyl alcohol and 0.2% NaOH, respectively. After adding extracting solution, the mixture was shaken for 30 min at 220 rpm, and was then centrifuged for 15 min at 4000 rpm. The supernatant was extracts, and the precipitate was resuspended with next extracting solution. Every extraction step was repeated for three times. Protein content in the supernatant was tested with an automatic azotometer (Kjelec 8400; FOSS, Denmark). The ratio of glutenins to gliadins (Glu/Gli) was calculated by dividing their contents. season. There were two biological replicates for grain sample at each stage. The midsection of ear with grain was freeze-dried in liquid N 2 , and grains were stored at −80 °C for RNA extraction. Grains were ground into find powder in liquid N 2 and total RNA was extracted as described by Plant RNA extraction kit V1.5 (Biofit, China, http://www.biofit.com). Three µg RNA per sample was used as input material for the RNA sample preparation. The mRNA samples were enriched by using oligo (dT)-magnetic beads and then cut into fragments with fragmentation buffer at 80 °C. First-strand cDNA was synthesized with random hexamers using the fragments as templates. Second-strand cDNA was subsequently performed using RNase H, DNA polymerase I, and dNTPs. The cDNA was purified by QiaQuick PCR kit and eluted with EB buffer. After terminal repair, poly A and adaptor sequences were connected to the cDNA end. cDNA libraries were established by PCR amplification after screening fragment size and assessing quality on Agilent 2100 Bioanalyzer system. The qualified cDNA libraries were sequenced on Illumina HiSeq2500 TM system at Gene Denovo Co., Ltd. (Guangzhou, China).
SCientifiC REPORTS | (2018) 8:11928 | DOI:10.1038/s41598-018-30451-4 Quantitative real-time PCR. To validate the reliability of expression profiles in RNA-seq data (N 0 VS N 225 ), 10 genes were selected for quantitative real-time PCR (qRT-PCR) analyses, using a SYBR premix Ex Taq tm RT-PCR kit (Takara, Dalian, China). Primers were designed by Primer Premier 5.0 (Premier Biosoft, Palo Alto, Canada) and listed in Supplementary Table S9. GAPDH (Genbank accession number KU246046.1), TUBA-2A (DQ435659.1), and TEF1 (M90077.1) were used as reference. The relative expression value was calculated by the ∆∆Ct method. Statistical analysis. All the data were calculated by using Excel 2010 (Microsoft, USA). Analysis of variance (ANOVA) was performed using the DPS software version 12.01 45 . Least significant difference (LSD) was used to compare difference among experimental mean values (P < 0.05). A bivariate correlation procedure was performed to analyze the relationships between these traits.