New transcriptomic insights in two RNAi wheat lines with the gliadins strongly down-regulated by two endosperm specific promoters

Mirim Mrín-Snz, Julio C.Msru Iehis, Frncisco Brro,*

a Department of Plant Breeding, Institute of Sustainable Agriculture (IAS), Spanish Council for Scientific Research (CSIC), 14004 Córdoba, Spain

b Department of Biotechnology, Faculty of Chemical Sciences, National University of Asunción, 111421 San Lorenzo, Paraguay

Keywords:Celiac disease Gliadin Protein compensation RNAi lines Triticum aestivum

ABSTRACT The gluten proteins of wheat grain are responsible for visco-elastic properties of flour,but they also trigger the immune-response of celiac disease.In this work,two low-gliadin RNA interference(RNAi)wheat lines that differ for the promoter driving the silencing(D-hordein and γ-gliadin promoters for D783 and D793 lines, respectively), were characterized at transcriptomic, and protein fraction levels in the grain.Silencing of gliadins provides a readjustment in the grain protein fractions that also affects to the nongluten proteins(NGP),which were increased in both RNAi lines.Determination of wheat gluten by means of the R5 monoclonal antibody also showed a strong reduction in the content of gluten in both RNAi lines.Moreover, fructans, an oligosaccharide linked with the development of non-celiac wheat sensitivity(NCWS) were also significantly decreased in RNAi lines.The down-regulation of gliadins fractions also impacts to other metabolic processes,particularly on carbohydrate metabolism,enzyme regulator activity and response to stress.Genes and transcription factors regulated by ABA were up-regulated, which could suggest the implication of this phytohormone on the stress response observed in the RNAi lines.

1.Introduction

The grain storage proteins of wheat are formed by monomeric and polymeric proteins that largely determine the viscoelastic properties of wheat flour.Polymeric proteins are denoted as glutenins, they are insoluble in alcohol as they form large polymers linked by interchain disulfide bonds,and are divided into two fractions; the high molecular weight (HMW) glutenins and the low molecular weight(LMW)glutenins.In contrast,gliadins are mainly monomeric proteins with alcoholic solubility,and they are further divided into three structural groups; ω-, α/β-, and γ-gliadins [1].Despite their role in breadmaking quality,gluten proteins also contain epitopes responsible for triggering the immuno-response in celiac disease (CD) [2] after ingestion of gluten-containing foods.CD is a chronic enteropathy characterized for lesions in the small intestine, flattering of the microvilli, hyperplasia of crypt cells and infiltration of leukocytes.The response is triggered by human leukocyte antigen (HLA) HLA-DQ2 and HLA-DQ8-presented peptides that are recognized by CD4 T cells (cluster of differentiation 4 T cells) [3].In addition to CD, other pathologies are also associated with wheat intake; allergies, which comprises the wheatdependent exercise-induced anaphylaxis (WDEIA) and baker’s asthma, and non-celiac wheat sensitivity (NCWS) [4].Wheat gluten is poorly digested in the human intestine and,when an individual eats food containing gluten, non-digested peptides cross the mucosa in the small intestine where the human transglutaminase 2(tTG2)deaminates glutamine residues present in gluten peptides into glutamic acid,increasing the binding affinity of these peptides to HLA-DQ2/8 molecules.After activation of T-lymphocytes by interaction with antigen-presenting cells,a spectrum of proinflammatory cytokines is released, damaging the enterocytes and producing the intestinal lesions typical of CD.The prevalence of CD is about 1% with a higher prevalence (1.5%) in Northern European countries [5].The prevalence of NCWS is still unknown due to cases of self-diagnoses and the lack of standardized diagnostic criteria or biomarkers, but globally it can affect a large population ranging from 0.6% to 13% [4].

The down-regulation of the immunogenic gliadin genes by RNA interference(RNAi)[6]provided wheat lines with low toxicity that could be used to prepare foodstuff with excellent sensory properties and health benefits for NCWS[7]and CD patients.These RNAi lines showed pronounced shifts in the different grain protein fractions that affect the quality,agronomic,and immunogenic properties [7-9].Lines D793 and D783 showed a strong reduction in the gliadin content of the grain, offering a reduced response to T-cell[6,10].These lines contain the same RNAi silencing fragment but driven by two different endosperm specific promoters.Interestingly, in response to the increment in nitrogen supply in both of them, silencing fragment was highly effective in the downregulation of α- and γ-gliadins, at whatever nitrogen supplied[11,12].Despite the strong down-regulation of gliadins, the total grain protein content was not significantly affected in comparison to the wild type,indicating a protein compensation with other protein fractions [11] and specifically, increments in non-gluten proteins (NGPs) like serpins, triticins and globulins were observed[13].However,the effects of RNAi gliadin silencing on other metabolic and physiological processes in the grain, which in the end,could regulate the compensation mechanisms is still unknown.

In this study, we report new insights in the regulation, transcriptome and trafficking of grain storage and non-storage protein genes in two wheat lines with all gliadin genes down-regulated by RNAi.Results reported are important for understanding the processes of protein regulation in the grain and for undertaking new strategies to redesign the wheat grain proteins by RNAi or CRISPR/Cas towards obtaining new varieties suitable for people suffering any gluten intolerance.

2.Material and methods

2.1.Plant material

Two RNAi lines, D793 and D783, derived from the wild type BW208 were used in this work.Lines D793 and D783 contain the plasmids pGhp-ω/α and pDhp-ω/α, respectively, harboring inverted repeat sequences of the most conserved regions from ω-, α/β- and γ-gliadin genes, driven by endosperm-specific promoters;γ-gliadin promoter for line D793 and D-hordein promoter for D783[14,15].The RNAi lines and the wild type were grown in a greenhouse with supplementary lights providing a day/night regime of 12 h/12 h at 24 °C/16 °C.

2.2.Prolamins determination

The gliadin and glutenin fractions were quantified by Reverse Phase High-Performance Liquid Chromatography (RP-HPLC)(1200 Series Quaternary LC System liquid chromatography from Agilent Technologies) with the RP8 column (5 μm in diameter and 250 mm in length,Merck).Three biological replicates for each genotype were analyzed following the protocol established in Pistón et al.[16].

2.3.Total protein, NGPs, gluten content by ELISA, fructans and starch determination

The protein content of whole flour was determined using the Kjeldahl (%N × 5.7) method according to the standard ICC no.105/2 (ICC, 1994), while the starch content was determined according to the standard ICC method no.123/1.(ICC, 1994).Both components were expressed on a 14% moisture basis.The NGPs content was calculated as the difference of total protein and total prolamin content for each line as described in previous work[12].Gluten content was determined by competitive ELISA assays using the monoclonal antibody R5 as described [17].

The fructan content per dry weight was calculated with Megazyme commercial kit following the fructan assay procedure (KFRUC, www.megazyme.com) using 150 mg of two biological samples and two duplicates per biological sample.After removing polisaccharides, the hydrolysis of fructan was performed to measure D-fructose content by spectrophotometric at 410 nm (iMark MicroplateAbsorbance Reader,Bio-Rad).To determine fructan content per dry weight, the moisture content was determined.

2.4.RNA extraction

For RNA sequencing (RNA-seq) data analysis, total RNA was extracted from seeds collected at 23 and 26 days after anthesis(DAA)from three different plants and bulked.For Real-Time quantitative Polymerase Chain Reaction(qPCR)analysis,seeds were collected at 8, 14, 18, 23 and 29 DAA from three different plants as biological replicates.After seed grinding, samples were placed in a tube containing extraction buffer (100 mmol L-1Tris-HCl, pH 9.5,150 mmol L-1NaCl,1.0%sarkosyl,5 mmol L-1DTT(Dithiothreitol), vortex mixed and centrifuged.Acid phenol (pH 4.3)-chloroform extraction was performed using the supernatant.RNA was precipitated with 1/10 vol of 3 mol L-1sodium acetate pH 5.2 and one volume of isopropanol.Resuspended RNA was extracted again using TRIzol reagent (Thermo Fisher Scientific,USA) and chloroform, and then with acid phenol-chloroform mixed with 1/3 vol of cold 3 mol L-1potassium acetate pH 5.5.The supernatant was mixed with salt precipitation solution(1.2 mol L-1NaCl and 0.8 mol L-1sodium citrate tribasic) and 0.25 vol of ethanol for RNA precipitation.

2.5.cDNA synthesis and qPCR experiments

cDNA was synthesized using 1 μg of total RNA and qScript cDNA SuperMix (Quantabio, USA) for RNA-seq and qPCR analyses.qPCR was carried out using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, USA) on CFX Connect Real-time PCR Detection System (Bio-Rad).PCR efficiency of each primer pair was determined for all samples with LinRegPCR version 2012.0[18].Normalization factor was calculated for each sample with geNorm [19] based on the expression levels ofRLI,CDC,andADP-RFgenes.The list of primers used is presented in Table S1.

2.6.RNA-seq data and bioinformatic analysis

cDNA for mRNA-seq data analysis was high-throughput sequenced by Fundación Parque Científico de Madrid (FPCM) (S/Faraday 7, Madrid, Spain) with Illumina Genome Analyzer (GAxII)platform.The quality control of total RNA (Agilent’s 2100 BioAnalyzer) and the library preparation of poly-A tails-mediated selection (Illumina TruSeq RNA kit) were performed by FPCM.The library type and sequencing details are in Table S2.Sequencing data can be found in the NCBI Sequence Read Archive under the following accession number:PRJNA593367.The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.Bioinformatic analysis of the RNA-seq was performed on a server with 64 cores and 128 Gb of RAM (Random Access Memory) with the operative system GNU/Linux Ubuntu version 18.04 (www.ubuntu.com, United Kingdom).The preprocessing of the RNA-seq data was carried out by removing the adapters with cutadapt (cutadapt v.1.16)[20] and trimming low quality reads, i.e.(1).The common dispersion(d=0.05096259)of a subset of housekeeping genes[27]was calculated.The final value used for DGE with edgeR was increased to 0.1,established for genetically identical model organisms [28] to reduce false positives in the data analysis.

2.7.Functional annotation

Functional annotation of genes was performed for known genes and denovoassembled genes.The last ones were annotated by blastx (blast-2.7.1+) against the non-redunant protein database of the National Center for Biotechnology Information (NCBI) and UniProtKB/Swiss-Prot knowledgebase [29].Genes that have been identified as uncharacterized genes by blastx method were subjected to a blastp search against UniProtKB/Swiss-Prot database,by obtaining protein sequence from transcripts using TransDecoder (TransDecoder 5.3.0) [30].In addition, gene ontology (GO)enrichment analysis was performed to study differential expressed(DE) genes in cell processes.Custom Python v2.7 scripts were implemented in RNA-seq workflow to facilitate data compatibility between software and annotation, and expression genes matrix SQLite database management.

3.Results

3.1.Grain composition, quality and agronomic features of RNAi lines

In this work,RNAi lines D783 and D793 with a strong reduction in α/β-, ω-, and γ-gliadins were characterized in comparison to their BW208 wild type.Both lines harbor the same RNAi silencing fragment but driven by different endosperm specific promoters; a D-hordein promoter for D783 and a γ-gliadin promoter for D793.A general overview of grain components and other quality and agronomic features between RNAi lines and the wild type is showed in Fig.1.All three gliadins fractions were decreased in both RNAi lines, particularly the γ- and α/β-gliadins.In contrast, the HMW glutenin subunits were increased in both RNAi lines but no the LMW, which was also decreased for D793.This readjustment in the grain protein fractions also affects the NGP,which were significantly increased for line D793.Overall,the total protein and starch contents, the major components of wheat grain, were slightly decreased in both RNAi lines in comparison to the wild type,although it was significant only for starch content in line D793.Determination of wheat gluten by means of the R5 monoclonal antibody also showed a strong reduction in the content of gluten in both RNAi lines.Moreover, fructans, an oligosaccharide linked with the development NCWS[31]were also significantly decreased in RNAi lines.Kernel weight, test weight and anthesis were not affected by the down-regulation of gliadins.In contrast, sodium dodecyl sulphate sedimentation (SDSS) test, a parameter used as predictor for wheat quality, was significantly affected in one RNAi line (D793) (Fig.1).

3.2.Transcriptomic analysis of grains from RNAi lines

Fig.1.Fold change of grain components and other quality and agronomic features for RNAi lines in comparison to that of the wild type.HMW,high molecular weight glutenin subunits;LMW,low molecular weight glutenin subunits;Gli/Glu,gliadins to glutenins ratio; SDSS, sodium dodecyl sulphate sedimentation; R5, gluten content(mg kg-1)determined by R5 monoclonal antibody.ANOVA was performed to determine statistically differences between genotypes for each component and feature.*, P < 0.05; **, P < 0.01; ***, P < 0.001.NS, Non-significant.

An RNA-seq data analysis was carried out for the two RNAi lines using a pool of grains samples at 23-26 DAA.Read libraries from the three samples were mapped to the wheat reference genome and DGE analysis was performed to identify DE genes between the three lines.In the pipeline, 17.4, 19.2, and 18.8 Mb of reads from BW208,D783,and D793,respectively,were uniquely mapped to the reference genome (Table S2).A total of 129,963 genes were considered:110,790 known genes and 19,173de novoassembly genes (identified as MSTRG name), of which 55,069, 75,773, and 52,144 were considered for DGE analysis in BW208 vs.D783,BW208 vs.D793 and D783 vs.D793 pair-wise comparisons,respectively (Table S2).There were a higher number of DE genes for the BW208 vs.D783 comparison than for the BW208 vs.D793 one(Fig.S2;Table S2).The GO enrichment analysis revealed that carbohydrate metabolic, response to stress or nutrient reservoir activity, among others, concentrate the higher number of DE genes,particularly for the BW208 vs.D783 comparison.In fact,carbohydrate metabolism term makes a remarkable difference between the RNAi lines in the number of DE genes, being this GO more enriched in D783 compared to the wild type (Fig.2).In contrast, nutrient reservoir activity was more enriched in D793 compared to the wild type (Fig.2).We also analyzed which GOs were enriched between up- and down-regulated genes in both RNAi lines compared to the wild type, to find out common responses between RNAi lines:(i) Nutrient reservoir activity was enriched in the set of up-regulated and down-regulated genes.(ii) However, response to stress was enriched in both RNAi lines among up-regulated genes while(iii)carbohydrate metabolic process was enriched among down-regulated genes (Fig.S3A, C).In the same way, we studied DE genes between D783 and D793 which were also DE in the BW208 vs.D793 comparison; nutrient reservoir activity was enriched among down-regulated genes because of the stronger silencing of gliadin genes in D793 line in comparison to that of D783 (Fig.S3B, D).A summary of DE genes in both RNAi lines with a high expression in almost one of the genotypes is provided in Table S3.Among them, serine-type endopeptidase inhibitor (Serpin) genes were over-represented among up-regulated genes, and some non-specific lipid transfer proteins (ns-LTPs) genes were DE in both RNAi lines (Table S3).According to the wheat reference genome, the major differences for gene expression (up or down-regulation) between RNAi lines and the wild type were observed for genes mapping in group 1,6,7 and Unassigned(Un)chromosomes,where gliadins and glutenins, among other genes, are located (Fig.S4).In contrast, the expression of serpins and ns-LTPs,mapping in group 4 and 5 chromosomes,was always higher in both RNAi lines(Fig.S4;Table S4).

Fig.2.GO enrichment analysis of differentially expressed (DE) genes in RNAi lines.All domains (biological process, cellular component and molecular function) are represented in three pair-wise comparisons:BW208 vs. D783, BW208 vs.D793, and D783 vs.D793.First column shows the P-value of pair-wise comparisons based on hypergeometric distribution.Green line represents the threshold for statistical significance (P < 0.001) of GO terms.Numbers on each column represent the number of DE genes in the comparison.

Next,hierarchical clustering of all DE genes allowed their gathering into four subclusters (Fig.S5), and the comparison of the expression patterns between the RNAi lines and the wild type.As shown, the wild type had clear differences for the expression pattern of DE genes to the RNAi lines.However, there are groups of genes in all subclusters with expression patterns notably different between the two RNAi lines,which could be due to different silencing pattern as they have different promoter driving the same RNAi fragment.Expression pattern in subcluster 1 showed an increased gene expression in D783 and D793,while a decreased gene expression was shown in subcluster 2 for both RNAi lines.Interestingly,GO enrichment analysis of genes grouped in subcluster 1 showed that they were mostly related to stress response,hydrolase activity and enzyme regulator activity, while in subcluster 2, nutrient reservoir activity was enriched.

3.3.Expression of prolamins and other grain protein genes

The expression of three different HMW genes (1Bx7,1Dx5, and1Ax2), encompassing the three bread wheat genomes, and three NGP genes were performed throughout grain development by qPCR analysis(Fig.3).The expression of all three HMW genes were higher in both RNAi lines from 14 to 26 DAA, particularly for line D793 with the maximum peak at 18 DAA.Data from RNA-seq analysis at 23-26 DAA also showed higher expression of the HMW genes for the RNAi lines (Fig.S6A, C).In addition, DE analysis of RNA-seq data confirmed the down-regulation of ω-, α/β-, and γgliadin genes in the RNAi lines (Fig.S6 C).

On the other hand, the expression of triticin and serpin genes was higher in both RNAi lines from 14 to 26 DAA while the globulin gene expression was similar for all three lines (Fig.3D).Interestingly, the expression pattern of all three genes was different throughout grain development.This increase in the expression of serpin and triticin genes observed by qPCR was confirmed in the RNA-seq data analysis in the case of serpin (Table S3; Fig.S6A, C).

There are other genes encoding proteins of interest involved in triggering wheat-related human diseases as α-amylase/trypsin inhibitors (ATIs).Although there were no DE ATI genes in D783 and D793 lines compared to BW208 (Fig.S6C), it can be appreciated that there are differences in gene expression between genotypes (Fig.S6C).

Fig.3.qPCR normalized expression values for three HMW glutenin subunit genes (1Bx7, 1Dx5, and 1Ax2), Glo3 (TraesCS4A01G296000, TraesCS4A01G296100), triticins(TraesCS1D01G067100) and serpin (MSTRG.33959) genes for BW208, D783, and D793 genotypes through grain development.Mean values are represented, and standard deviation was calculated with three biological replicates.DAA,days after anthesis;HMWGS,high molecular weight glutenin subunit;BW,BW208.ANOVA was performed for each DAA between each RNAi line and the wild type.*, P < 0.05; **, P < 0.01; ***, P < 0.001.

Next,we examined the expression of transcription factors(TFs)related to prolamin genes.In Table S5 a comprehensive list of TFs related to seed storage proteins found in the RNA-seq and their homologs/homeologs in barley, wheat or rice is provided.There was only one TF(TraesCS2D01G117000),homeolog to OsGZF1 from rice, that was significantly up-regulated in D783 in comparison to the wild type.In rice, this OsGZF1 TF participate in the regulation of the seed storage protein GluB-1 [32].

3.4.Expression of genes related to stress response in RNAi lines

As showed previously, many genes related to stress response were enriched in the RNAi lines (Fig.2).For example, dehydrins and Rab genes (response to abiotic stimulus) were highly expressed in both RNAi lines (Table S3).Among TFs related with defense response, or abiotic and biotic stress response with high fold-change (FC) in the RNAi lines, 12 were associated with response to abiotic, biotic and/or abscisic acid (ABA) (Table 1),and some of them were significantly up-regulated in both RNAi lines but, particularly in line D783.

Table 1 Putative transcription factors (TFs) related with defense response, and abiotic and biotic stress response.

Next, we set up qPCR analysis for some of the DE genes related to stress throughout grain development (Fig.4).The expression of dehydrin (Dhdn4) gene was significantly higher in D793 for some of the last days of grain development and, from 23 DAA, D783 showed no significantly higher expression than that of the wild type (Fig.4 A).The expression of the ABA INSENSITIVE 5 (ABI5)TF was also analyzed because of its relation to ABA-dependent signaling pathway and its role in abiotic stress response [40,43].The expression of this gene was higher in both RNAi lines compared to the BW208 line,but only significant for D793(Fig.4 B).Finally,the basic Leucine Zipper (bZIP) TF (homolog to rice OsbZIP20/RITA-1/RISBZ3) was reported as implicated in the ABA-dependent stress response [36,37] and in the transcriptional activation of prolamin genes and other storage protein genes as globulin in rice [44].The RNA-seq and qPCR data analysis confirmed that this TF was up-regulated in both RNAi lines.This bZIP TF present higher expression levels in the RNAi lines than the wild type for most of the grain development stages, significantly higher for both RNAi lines at 18 and 26 DAA (Fig.4 C).

Fig.4.qPCR normalized expression values for Dhdn4 (dehydrin 4) (MSTRG.41119), ABI5 (TraesCS3A01G371900, TraesCS3A01G372200) and bZIP (TraesCS7B01G391800),ortholog to rice OsbZIP20/RITA-1/RISBZ3 through grain development.Mean values are represented, and standard deviation was calculated with three biological replicates.DAA, days after anthesis; Dhdn4, dehydrin 4; ABI5, ABA INSENSITIVE 5; BW, BW208.ANOVA was performed for each DAA between each RNAi line and the wild type.*,P < 0.05; **, P < 0.01; ***, P < 0.001.

Fig.5.qPCR normalized expression values for three NAC TFs through grain development,orthologs to NAC92(XM_020326243.1 and XM_020341606.1)from Aegilops tauschii,to ANAC092 (AT5G39610.1) from Arabidopsis thaliana, and to ONAC20 (Os01g0104500) from rice.Mean values are represented, and standard deviation was calculated with three biological replicates.DAA,days after anthesis;BW,BW208.ANOVA was performed for each DAA between each RNAi line and the wild type.*,P<0.05;**,P<0.01;***,P < 0.001.

3.5.Expression of genes related to carbohydrate metabolism and other metabolic process in RNAi lines

Both RNAi lines show a slight decrease in the grain starch content, only significant for line D793 compared to the wild type(Fig.1).Carbohydrate metabolism process was enriched for both RNAi lines, especially in D783 with up to 111 DE genes (Fig.2).Some of the genes implicated in starch synthesis were up- or down-regulated in both RNAi lines, as starch synthase, betaamylase, starch branching enzyme and sucrose synthase among others (Table S3).

We also monitored de expression of three NAC TF genes for the RNAi and the wild type lines throughout seed development.In addition to response to stress, NAC TFs are also related to many important metabolic processes as seed development, starch synthesis and grain nitrogen concentration [45-48].We designed the primers based onTraesCS7B01G094000, the homeologsTraesCS7D01G154200andTraesCS7B01G056300, andTraesCS7A01G569100.These genes are orthologs toNAC92(XM_020326243.1 and XM_020341606.1) fromAegilops tauschii,toANAC092(AT5G39610.1) fromArabidopsis thaliana, and toONAC20(Os01g0104500) from rice.ONAC20is related to seed development and seed size/weight [47] andNAC92participates in many processes as regulation of seed germination,anther development, age-related resistance and stress response among others[49,50].The expression of these NAC TFs is shown in Fig.5.For all of them, higher expression in the RNAi lines in comparison to that of the wild type was found during grain maturation, significant in many of the stages for both RNAi lines.In the RNA-seq data,the expression of three of them (TraesCS7D01G154200,TraesCS7B01G056300, andTraesCS7A01G569100) was also higher in the RNAi lines, although it was only significant forTraesCS7D01G154200in D783 compared to the wild type.

4.Discussion

Obtaining wheat lines lacking the immunogenic component of gluten is a very appealing goal as an alternative to produce foods for people that need to follow a gluten-free diet.The RNAi technology has demonstrated to be highly effective in the downregulation of wheat gliadin genes [6].In addition to RNAi,CRISPR/Cas9 technology was also used for the reduction of wheat gliadins [51].In this work, two RNAi lines, D783 and D793, with the same silencing fragment but different endosperm specific promoter, were evaluated for protein fractions, and transcriptomic features.Results reported here, confirmed the strong downregulation of gliadin proteins in both RNAi lines [6,8].However,this down-regulation was lower for ω-gliadins, indicating that for this gliadin fraction the silencing construct is not as effective as for the α/β-and γ-gliadins.In fact,in a nitrogen fertilization study using these two RNAi lines, α/β- and γ-gliadins were silenced in both lines, even at high nitrogen concentrations.However, the ωgliadins were increased when increasing the nitrogen, confirming the lower effectiveness of silencing [11,12].

Wheat lines with down-regulation of one or more gliadin fractions were reported to present a compensatory mechanism involving other grain proteins [52,9] providing nitrogen contents comparable to that of the wild types.In this work, the increment of HMW in both RNAi lines is the major responsible of the increasing of the total glutenin content as the LMW had a different behavior in D783 and D793 lines, with a significant decrease for line D793.Similar changes were previously reported for RNAi lines with all the gliadin fractions down-regulated[11]where the readjustment of the protein composition is not the same among the RNAi lines.In the qPCR data, the increase of HMW transcripts abundance occurs at 18 DAA in both RNAi lines, and it could account for the higher HMW protein content observed in the grain for lines D783 and D793.This increment in the HMW glutenin content in both RNAi lines is important as they contribute to the viscoelastic properties of dough, and could be related to the unique bread-making quality features observed in these RNAi lines,particularly the tolerance to over-mixing[8].The synthesis of both LMW and HMW is essentially regulated at the transcriptional level, and therefore,differences for the HMW gene expression between RNAi lines and the wild type could be attributed to the transcriptional regulation.In the stage of development in which we have carried out the RNA-Seq, we have found that gene TraesCS2D01G117000,ortholog toOsGZF1 in rice, and described as regulator of the seed storage protein GluB-1[32]was up-regulated in RNAi lines,particularly D783.In addition, we obtained higher expression of a bZIP wheat genes orthologs toOsbZIP20/RITA-1/RISBZ3gene from rice,which is involved in transcriptional activation of prolamin genes and other storage protein genes as globulins [44].

For NGPs,such as globulins,serpins,and triticins,their increase was reported to compensates the down-regulation of gliadins to maintain the total grain protein content [13,53].An increase in the number and size of the inclusions at the surface of protein bodies was also reported for the RNAi low gluten wheat lines [53].These inclusions contain triticins that raise the 5% content of total grain protein [53,54].García-Molina et al.[9] reported that serpin spots were increased in the RNAi low gluten wheat lines.The qPCR and RNA-seq analysis agreed to protein reported data and confirmed that the higher content in the NGPs is related to a higher gene expression during grain development.Serpins are related to wheat allergy as they contain Immunoglobulin E(IgE)-binding epitopes, and therefore, the increment in this protein fraction in the RNAi lines may indicate the increment in wheat allergy properties.Serpins have also a role in response to abiotic and biotic stresses[55] and they have been described as possible prolamin proteases inhibitors[56]so it could be possible that a regulation mechanism against the decrease of gliadins has been activated.

In addition to grain proteins, there are other wheat grain components also responsible for allergies and intolerances.One of them are fructans (oligo- or polysaccharides with short chain of fructose units and terminal glucose), which are hydrolyzed partially in the intestine causing symptoms as bloating and abdominal pain.These belong to fermentable oligosaccharides, disaccharides,monosaccharides and polyols(FODMAPs)and are related to NCWS[57].Fructans are used as reserve carbohydrates with similar function as starch and sucrose.Although genes encoding fructan degrading enzymes were not DE in both RNAi lines in comparison to that of the wild type,at least between 23 and 26 DAA,the fructan content in mature grains of RNAi lines was lower than the wild type, in the same way as total starch decreased.In the RNA-seq data analysis, carbohydrate metabolic process was enriched for both RNAi lines,and was more significant for D783 compared with the wild type,involving some DE carbohydrate synthesis enzymes genes such as sucrose synthase (SS), starch branching enzyme(SBE)and starch synthase(SSS).Interestingly,D783 and D793 have different behavior in terms of up- or down-regulated status of these genes compared to the wild type.This indicates that the variation in fructan and starch contents during kernel maturation could be part of signaling pathways to regulate carbohydrate metabolism and storage [58].

Other genes related to wheat pathologies are ns-LTPs and ATIs which show high reactivity with IgE related to baker’s asthma[59].In this work, ns-LTPs were up-regulated in both RNAi lines while the expression of ATIs genes at 23-26 DAA was not affected by the silencing of gliadins.However, proteomic studies showed that ATIs were over-accumulated in RNAi lines, particularly in D793 line [9].This clearly indicates that the expression data obtained by RNA-Seq at 23-26 DAA for ATIs are not predictive of the protein readjustments that finally operate in the grain.The increase in ATIs and serpins,both related to allergies,make necessary further studies, particularlyin vivostudies to determine the possible impact of these lines on allergies.

The down-regulation of gliadins fractions also impacts to other metabolic processes.Two of the top enriched GOs in our analysis were response to abiotic and biotic stresses.Both of them comprise many genes as dehydrin,rab protein gene,the rRNA N-glycosidase,wheatwin-1 and subtilisin-chymotrypsin inhibitor WSCI among others, related to drought stress and plant defense against pests,and bacteria and fungi pathogens [60-66].These genes were upregulated in both RNAi lines, mainly in D783.Moreover, some TFs that regulate the stress response, were also up-regulated in the RNAi lines, as DREB TFs that regulate dehydrin genes [33].Among these TFS, the ABI5 stands out as having a key role in the stress response in presence of ABA [43] and it has higher gene expression in both RNAi lines, particularly in D793 line as reveled also by qPCR.ABA participates in the stress response by regulating the expression of many genes which have a key role in this response, being part of ABA-dependent stress response process[67].In the present work, some genes regulated by ABA were upregulated, which could suggest the implication of this phytohormone on the enrichment of stress response metabolic processes observed in the RNAi lines.For example, dehydrin genes has ABA-responsive elements (ABRE) motifs in their promoter sequence in wheat [68]; the pathogen defense Wheatwin-1[69,70]a class II chitinase of PR-4 family with a high sequence similarity with ZmPR4 protein fromZea mays,could be up-regulated by ABA, as it was reported that ABA can induce the expression ofZmPR4[71]; ABI5 and NAC TFs have a role in ABA-mediated response to abiotic stresses [67,43]; the NAC TF coding geneTraesCS5A01G468300, up-regulated in RNAi lines, could be also involved in ABA-dependent response as its homolog in rice(OsNAC19) is induced by this phytohormone, providing resistance to blast fungus [38].In addition, other grain components as ns-LTPs and HMW have been studied as genes transcriptionally regulated by this phytohormone in barley and wheat [72,73].In a recent RNA-seq analysis in wheat under heat stress condition, γgliadin and LMW were up-regulated across 13 and 30 DAA, suggesting a relation between these prolamin genes and the heat stress condition [74].However, in our experimental conditions we have not found an increment in the expression of 9-cisepoxycarotenoid dioxygenase and ABA 8′-hydroxylase, main enzymes involved respectively in biosynthesis and catabolism of ABA[67,75].Although activation of response to stress mechanisms could be a consequence of the down-regulation of gliadins by RNAi lines, more studies are required to determine the molecular and physiological basis of this response and the possible implication of ABA.However,the pre-activation of these stress response mechanisms in these lines is a very interesting aspect since it could provide them advantages in certain conditions where a rapid response to stress situations is needed.

5.Conclusions

Wheat lines D783 and D793 have a strong down-regulation in the gliadin fraction, which are compensated by NGPs and HMW glutenins since they are increased in D783 and D793 lines.This compensation with NGPs, particularly serpins and ATIs, could somehow increase allergies.However, these lines also had other transcriptomic changes compared to the wild type, involving enzyme regulator activity,carbohydrate metabolism,and response to stress.Carbohydrate metabolic process showed the downregulation of some key enzymes in starch biosynthesis that could provide lower starch content in RNAi lines.Finally, the silencing of wheat gliadins has a clear impact on the activation of genes and TFs related to the response to stress, leaving a much more complex metabolic restructuring, and that ultimately would also be responsible for the differences in agronomic traits such as kernel weight and/or yield.Many of these metabolic processes are regulated by ABA, which would indicate that this phytohormone may be involved.

Availability of data and materials

Sequence data from this article can be found in the NCBI Sequence Read Archive under the following accession number:PRJNA593367.The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

CRediT authorship contribution statement

Francisco Barro:conceptualization,funding acquisition,supervision, writing - original draft, writing -review & editing.Miriam Marín-Sanz:investigation, methodology, writing - original draft,writing - review & editing.Julio C.Masaru Iehisa:investigation,methodology, writing - review & editing.

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

Technical assistance from Ana García is acknowledged.This research was funded by the Spanish Ministry of Science and Innovation(Project PID2019-110847RB-I00)and the European Regional Development Fund (FEDER).

Appendix A.Supplementary data

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