System analysis of differentially expressed miRNAs in hexaploid 1 wheat display tissue-specific regulatory role during Fe deficiency 2 response 3

28 Background 29 Iron (Fe) is an essential mineral element, and its deficiency in soil largely affects crop 30 productivity. In plants, the molecular mechanisms underlying the genetic regulation of Fe 31 deficiency responses have yet to be well understood. Specifically, microRNA (miRNA) 32 mediated regulation of Fe deficiency response and its regulatory network is largely elusive. 33 In the current work, we utilized a whole genome transcriptomic approach to identify the Fe 34 deficiency-responsive miRNAs to understand the molecular mechanisms of Fe deficiency 35 response in wheat seedlings. The study also identifies nine novel miRNAs putatively 36 involved in Fe deficiency response. Further, the identified miRNAs showed tissue 37 preferences relating them to differential mechanisms against Fe deficiency. 38 Results 39 In the present study, we performed small RNA-targeted whole genome transcriptome 40 analysis to identify the involvement of sRNAs in Fe deficiency response. The analysis 41 identified 105 differentially expressed miRNAs corresponding to Fe deficiency response, 42 among them, 9 miRNAs were found to be novel in this study. Interestingly, tissue-specific 43 regulation of Fe deficiency response also participates through miRNA-mediated regulation. 44 We identified 17 shoot specific miRNAs and 18 root-specific miRNAs with altered 45 expression. We validated the tissue specificity of these miRNAs by stem-loop quantitative 46 RT-PCR. Further, an attempt was made to predict their targets to speculate their participation 47 in Fe deficiency response. This miRNA target prediction analysis suggested a few major 48 targets of the identified miRNAs, such as multicopper oxidases, E3 ubiquitin ligases, GRAS 49 family, and WRKY transcription factors previously known to play key roles in Fe 50 homeostasis. Our analysis of selected miRNAs also confirmed a temporal regulation of the 51 response. 52 Conclusion 53 The first information generated here will classify the repository of wheat miRNAs (with few 54 novel miRNAs) for their role in Fe deficiency response. Our work provides insights into

16.03 million reads were obtained for +Fe and -Fe root libraries. After removing low-quality 129 reads, 14.31 million and 13.72 million clean reads yielded for +Fe and -Fe shoot, whereas 130 14.77 million and 15.77 million for +Fe and -Fe root, respectively. After further refining the 131 sequence reads, we ended with 6.4 million to 10.01 million total sRNA reads with 1.74 132 million to 2.28 million unique reads encompassing different types of small RNAs (Table S1). 133 Chromosomal distribution of the mapped sRNA reads on the wheat reference genome 134 showed a predominant contribution from A and B genomes compared to D genomes ( Figure  135 S2). Interestingly, all the genomes of chromosome 7 contributed to the total sRNA reads. 136 Interestingly, we observed that the sense strand of the genome is more responsive 137 towards Fe deficiency as more than 61% of the sRNA reads observed were coded from the 138 sense strand, while the antisense strand was only able to contribute for around 5-10% of the 139 fraction, irrespective of the tissue (Table S1; Figure S2). Further, we observed a decrease in 140 sRNA reads coded by both sense and antisense strands in the shoot in response to Fe 141 deficiency (4.46% and 1.03% decrease for sense and antisense strands, respectively). In 142 contrast, to shoot, there was an increase in sRNA reads in roots in response to Fe deficiency, 143 preferably over-representing sRNAs coded by antisense strand (0.2% and 1.48% for sense 144 and an antisense strand, respectively). 145 To gain an insight into sRNA distribution in the genome, we analyzed whether the 146 coding region or the non-coding region is responsive to sRNA-mediated regulation of Fe 147 deficiency response. Further, the positional mapping of generated sRNA reads from the root 148 and shoot accounted for 11.15 % and 6.24 % of reads mapping to exon and 2.54 % and 3.90 149 % to intron, respectively (Table S2). 150 151

Identification of Fe deficiency-induced miRNAs 152
To characterize the sRNA reads into different subfamilies, we annotated all the sRNAs with 153 Rfam database into rRNA, snoRNA, snRNA, TAS, and miRNA classes (Figure 1A-D). This 154 analysis extended our observation of inverse relations in sRNA reads in shoot and root tissues 155 in response to Fe deficiency. This categorised data showed a significant decrease in miRNA 156 representation in response to Fe deficiency in shoot tissues. In contrast, we observed 157 increased miRNA reads in response to Fe deficiency in root tissues (Figure 1A-D). This data 158 suggested that miRNAs might be acting in a tissue manner in regulating Fe deficiency 159 response in wheat. 160 While analysing the length-based classification of unique sRNA reads, we found that the 20-161 24 nucleotide sRNAs were the most abundant classes in our datasets, representing around 5-162 30% of our data irrespective of the treatment conditions ( Figure 1E). To our interest, the 163 abundance pattern of sRNAs observed in shoot and root tissues w.r.t. to Fe deficiency was 164 also found to be length biased. We observed that 24 nucleotide-long sRNAs follow decreased 165 abundance in shoot while showing an increased abundance in root tissue in response to Fe 166 deficiency ( Figure 1E). 167 As we observed that the 20-24 nucleotide long sRNAs are over-represented in our 168 dataset, which typically lies in the range of miRNAs, we analysed these reads for 169 characteristics of miRNAs. Our analysis identified 105 miRNAs in our transcriptome 170 analysis, out of which nine were novel ( Table S3). As evident from previous reports, most 171 active miRNAs prefer U at the first nucleotide at 5' end, which in addition to high A+U 172 content, provides them low internal stability, promoting them to be processed into mature 173 miRNA through RISC complex (RNA-induced silencing complex). Additionally, A or U at 174 the 10 th position is over-represented in natural plant miRNAs, further contributing to their 175 processibility [39,40]. In agreement of these previous reports, we observed that around 27% 176 of the identified miRNAs showed first nucleotide preference for U and around 45% of 177 miRNAs have nucleotide preference for A/U at 10 th position (Table S3). Hairpin analysis of 178 these identified miRNAs classified them into 36 miRNA encoding families. Interestingly, 179 members for 35 out of 36 hairpin families are represented in the wheat genome, which was 180 the maximum variability observed among all the 66 plant species analysed (Table S4). To 181 further validate our report, we predicted the secondary structure of the identified miRNAs 182 with an RNAfold web server with a minimum free energy index (MFEI) algorithm 183 (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi ). Characteristic stem-loop 184 hairpin formation in all 9 novel miRNAs validated their secondary structure and strengthened 185 their putative function as miRNAs ( Figure S3). 186 Next, to pinpoint miRNAs responsive to Fe deficiency, miRNAs were profiled for 187 differential expression patterns using DEGseq in all four libraries. Our analysis revealed that 188 23 miRNAs were upregulated in the shoot; 17 were shoot specific. Further, 20 miRNAs 189 showed downregulation, with 10 miRNAs downregulated in a shoot-specific manner ( Figure  190 2). We extended a similar analysis for root tissues where we observed that 15 miRNAs 191 showed upregulation and 19 miRNAs showed downregulation with 8 and 10 miRNAs with 192 root-specific behaviour (Figure 2). This relative analysis further allowed us to identify the 193 miRNAs with similar regulation irrespective of the tissues. We found 3 miRNAs (tae-194 miR1122b-3p, -1138 and -9652-5p) commonly upregulated while 6 miRNAs (tae-195 miRNA5049-3p, novel_1, 395b, 395a, 408 and 9666b-3p) commonly downregulated in either 196 tissue in response to Fe deficiency. Further, 7 miRNAs (tae-miRNA9679-5p, -9782, -1118, -197 1117, -9674a-5p, -9654a-3p and novel_10) showed inverse behaviour in terms of their 198 transcriptional induction during Fe deficiency when shoot and root tissues are compared 199 (Figure 2). 200 201

Analysis of miRNA for their temporal expression responses 202
In order to validate our transcriptomic results, we checked the expression of miRNAs in both 203 the tissues. We randomly selected eight miRNAs from our transcriptomic data and performed 204 stem-loop qRT-PCR analysis with gene-specific primer sets (Table S5). Following this 205 analysis, we found that the expression pattern of all the selected miRNAs was very similar to 206 the one observed in transcriptomic analysis with Pearson's correlation coefficient 207 (R 2 =0.9803) (Figure 3). 208 To further understand miRNA-mediated Fe deficiency response temporally, we 209 studied expression patterns of multiple miRNAs in either of the tissues after subjecting them 210 to Fe deficiency stress for 6, 9, 12 and 15 days. We selected miRNAs found differentially 211 regulated in the wheat shoot (tae-miR6201, -5050, -9774, -1122a, -1137b-5p and -9671-5p) 212 and root (tae-miR1138, -167c-5p, -444a, -9652-5p, -9654a-3p and -397-5p) tissues. Although 213 we found a perfect correlation between our transcriptome and stem-loop qRT-PCR data, 214 varying temporal transcriptional responses of miRNAs were observed during early and 215 prolonged Fe deficiency stress, irrespective of the tissue considered (Figure 4). This 216 approach helped us to characterize the tissue-specific expression responses into three 217 categories. The first includes early responsive miRNAs, where most of the miRNAs showed 218 early induction in the shoot, though, in the root, few miRNAs showed a low early response 219 (tae-miR444a, -9654a-3p and -397-5p). The second category comprises late responsive 220 miRNAs mainly accumulated in root tissues, including (tae-miR444a, -167c-5p, -9654a-3p 221 and -397-5p). Additionally, in shoot tissues, most of them showed weak late responsive 222 nature except two miRNAs (tae-miR1137b-5p and tae-miR9671-5p) which showed 223 constantly induced and late responsive behaviour, respectively, in response to Fe deficiency 224 (Figure 4). The third category highlights the miRNAs with mixed temporal expression 225 during Fe deficiency, which included tae-miR9774 in shoot and tae-miR1138 and tae-miR- partially increased the number of miRNAs expressed ( Figure S4B). Therefore, though we 239 did not find any negative impact of the DD genome, the AABB genome contributed most to 240 the expression of the selected miRNAs. Although, when we analysed the expression levels of 241 selected miRNAs in different genomes, the AABBDD genome contributed the most to the 242 expression levels of miRNAs (48%), while the AABB genome was the least contributor 243 (21%) while the DD genome contributed around 31% for the expression levels of the 244 miRNAs ( Figure S4C). Therefore, in conclusion, we observed that the DD genome has the 245 lowest penetrance in miRNA expression levels and the highest expressivity.

Target prediction of wheat miRNAs revealed an adaptive response against Fe deficiency 282
To identify the regulatory network involved in Fe-homeostasis by these distinct miRNAs, we 283 predicted their targets and analysed them for their involvement in Fe deficiency response. 284 Among the most interesting targets for these Fe deficiency responsive miRNAs revealed a 285 complex network involving multicopper oxidases (MCOs), transcription factors like GRAS 286 and MADS-box, major facilitator superfamily (MFS) transporters, E3 ubiquitin-protein 287 ligases, oxidoreductases, protein kinases etc.. In prediction analysis, we observed that 288 miRNA397-5p could target transcripts encoding for MCOs. Tae-miRNA171b was predicted 289 to target transcripts encoding for GRAS TF. On similar lines, sulphate transporters are 290 predicted to be targeted by Fe deficiency-induced miRNA395a and 395b. As we proposed to 291 understand the miRNAs mediated regulatory network against Fe deficiency, we predicted the 292 differential response of shoot and root-expressed miRNAs in terms of their targets. 293 Interestingly, we observed that among the targets for root-induced miRNAs, GRAS family 294 TFs, NAC and MADS TFs, SCR-like genes, serine/threonine phosphatases, and sugar 295 transporters. 296 While among the targets of root down-regulated miRNAs included potassium 297 transporter, CNX (molecular chaperon), bHLH TFs, MFS transporters, Zn finger TFs, F-box 298 related, sulphate transporters, kinases, cell wall-related genes. Our analysis supported our 299 previous report suggesting over-accumulation of bHLH TFs during Fe deficiency [36]. 300 miRNA-regulated expression of CNX type of molecular chaperones suggested an adaptive 301 response for root against the detrimental effect of ROS accumulation during Fe deficiency. 302 Apart from it, we observed that a significantly downregulated miRNA in the root (tae-miR 303 5049-3p) was found to target S-adenosyl-L-methionine-dependent methyltransferase while 304 the S-adenosyl-methionine is the precursor for mugineic acid (MAs) family of PS [41]. This 305 indicates that PS biosynthesis in wheat might be regulated through tae-miR 5049-3p. In the 306 shoot, however, we observed that targets associated with redox enzymes, kinases and 307 phosphatases were over-represented. Among the TFs families, NAC and myb were among 308 the targets of Fe-responsive miRNAs and phytohormone (like auxin and JA) associated genes 309 were also targeted through these miRNAs. Therefore, the target prediction analysis, on the 310 one hand, strengthens the involvement of miRNAs in Fe deficiency response in wheat. At the 311 same time, it also suggests a tissue-specific regulation of their targets. to different abiotic stresses. However, an attempt has yet to be made to unravel the miRNA-319 mediated control of Fe homeostasis in hexaploid wheat. Although genes involved in Fe 320 homeostasis were reported earlier, the miRNA-mediated targets were not addressed in wheat 321 [36]. Noteworthy, such molecular responses largely depend on the genotype and specific 322 stress condition. This study attempts to identify wheat miRNA that could be linked to the 323 target pathway functions to identify the critical regulatory miRNA-mRNA interaction 324 involved in Fe deficiency conditions. A total of 105 miRNAs were identified in shoots and 325 roots, respectively. Our work identified 9 novel miRNAs with distinct expression responses 326 during Fe deficiency with a typical stem-loop structure ( Figure S2). 327 In this study, sRNA libraries were generated from the roots and shoots of wheat 328 seedlings subjected to Fe deficiency for different time points. This was done to generate the 329 inventory of miRNA that could assist in collating miRNA that may differentially express at 330 any time. Our analysis resulted in the identification of multiple root and shoot specific 331 miRNAs in response to Fe deficiency ( Figure 2; Table S3). Among these differentially 332 expressed miRNAs, a subset of randomly selected miRNAs from root and shoot were 333 employed to validate the RNA-Seq data. qRT-PCR of these tissue specific miRNAs drew a 334 strong correlation with respect to RNA-Seq data. . Furthermore, the characterization of 335 differentially expressed miRNAs from root (tae-miR1138, -167c-5p, -444a, -9652-5p, -336 9654a-3p and -397-5p) and shoot (tae-miR6201, -5050, -9774, -1122a, -1137b-5p and -9671-337 5p) revealed the spatio-temporal expression of these Fe-responsive miRNAs. 338 Fe deficiency is a major nutritional disorder that limits crop productivity. In plants, 339 multiple miRNA gene families are known to be involved in the Fe deficiency responses 340 [27,47]. In addition, correlation studies were done where the expression of specific miRNA 341 was observed in high and low Fe genotypes of wheat and rice [30,48]. These studies support 342 that miRNA-mediated control could occur concerning the Fe flux in a tissue-specific manner. 343 In our analysis size of the majority of the filtered reads was 21-24 nt ( Figure 1E) analysis tae-miR397-5p shows potential targets for multiple wheat MCOs (Table S6) Altogether, miRNA profiling suggests their involvement in regulating Fe deficiency 412 responses. Our work provides evidence that miRNA perturbed due to Fe deficiency targets a 413 subset of previously reported Fe-responsive genes (Figure 7), This reflects that miRNA-414 mediated control of Fe-responsive genes contributes to such regulatory mechanism in 415 hexaploid wheat. Specifically, the miRNA can target the genes primarily involved in 416 changing the Fe redox status and its uptake. Similarly, the regulation of multiple TFs was 417 also targeted by miRNA (Figure 7). Overall, the generated datasets will serve as an

Bioinformatics analysis and miRNA identification 453
The data obtained from high throughput sequencing was processed into raw sequencing reads 454 by CASAVA base recognition. Low-quality raw reads containing adaptors, 5' primer 455 contamination, and polynucleotide tails, and reads with >50% bases having a Qphred less 456 than or equal to 5, and the ones in which >10% base information were indeterminable were 457 discarded. The clean reads ranging from 18 to 30 nucleotides were mapped to the reference 458        Table S1. Small RNAs obtained in Control (Fe-EDTA 80µM) and -Fe treated (Fe-EDTA 844 2µM) Wheat C-306 small RNA libraries of root and shoot 845