Genome-Wide Identification of MicroRNAs in Leaves and the Developing Head of Four Durum Genotypes during Water Deficit Stress

MicroRNAs (miRNAs) are small non-coding RNAs that play critical roles in plant development and abiotic stress responses. The miRNA transcriptome (miRNAome) under water deficit stress has been investigated in many plant species, but is poorly characterised in durum wheat (Triticum turgidum L. ssp. durum). Water stress during early reproductive stages can result in significant yield loss in durum wheat and this study describes genotypic differences in the miRNAome between water deficit tolerant and sensitive durum genotypes. Small RNA libraries (96 in total) were constructed from flag leaf and developing head tissues of four durum genotypes, with or without water stress to identify differentially abundant miRNAs. Illumina sequencing detected 110 conserved miRNAs and 159 novel candidate miRNA hairpins with 66 conserved miRNAs and five novel miRNA hairpins differentially abundant under water deficit stress. Ten miRNAs (seven conserved, three novel) were validated through qPCR. Several conserved and novel miRNAs showed unambiguous inverted regulatory profiles between the durum genotypes. Several miRNAs also showed differential abundance between two tissue types regardless of treatment. Predicted mRNA targets (130) of four novel durum miRNAs were characterised using Gene Ontology (GO) which revealed functions common to stress responses and plant development. Negative correlation was observed between several target genes and the corresponding miRNA under water stress. For the first time, we present a comprehensive study of the durum miRNAome under water deficit stress. The identification of differentially abundant miRNAs provides molecular evidence that miRNAs are potential determinants of water stress tolerance in durum wheat. GO analysis of predicted targets contributes to the understanding of genotypic physiological responses leading to stress tolerance capacity. Further functional analysis of specific stress responsive miRNAs and their interaction with targets is ongoing and will assist in developing future durum wheat varieties with enhanced water deficit stress tolerance.

MicroRNAs (miRNAs) are small non-coding RNAs that play critical roles in plant development and abiotic stress responses. The miRNA transcriptome (miRNAome) under water deficit stress has been investigated in many plant species, but is poorly characterised in durum wheat (Triticum turgidum L. ssp. durum). Water stress during early reproductive stages can result in significant yield loss in durum wheat and this study describes genotypic differences in the miRNAome between water deficit tolerant and sensitive durum genotypes. Small RNA libraries (96 in total) were constructed from flag leaf and developing head tissues of four durum genotypes, with or without water stress to identify differentially abundant miRNAs. Illumina sequencing detected 110 conserved miRNAs and 159 novel candidate miRNA hairpins with 66 conserved miRNAs and five novel miRNA hairpins differentially abundant under water deficit stress. Ten miRNAs (seven conserved, three novel) were validated through qPCR. Several conserved and novel miRNAs showed unambiguous inverted regulatory profiles between the durum genotypes. Several miRNAs also showed differential abundance between two tissue types regardless of treatment. Predicted mRNA targets (130) of four novel durum miRNAs were characterised using Gene Ontology (GO) which revealed functions common to stress responses and plant development. Negative correlation was observed between several target genes and the corresponding miRNA under water stress. For the first time, we present a comprehensive study of the durum miRNAome under water deficit stress. The identification of differentially abundant miRNAs provides molecular evidence that miRNAs are potential determinants of water stress tolerance in durum wheat. GO analysis of predicted targets contributes to the understanding of genotypic physiological responses leading to stress tolerance capacity. Further functional analysis of specific stress responsive miRNAs and their interaction with targets is ongoing and will assist in developing future durum wheat varieties with enhanced water deficit stress tolerance.

Introduction
Durum wheat (Triticum turgidum L. ssp. durum) is the only tetraploid wheat species (2n = 4x = 28, genomes AABB) grown commercially throughout the world. Water deficit stress is one of the main abiotic factors that cause durum yield loss in Mediterranean environments. Water deficit stress in early reproductive stages has been shown to adversely affect grain yield and biomass through reduced grain number in durum [1]. Nonetheless, Liu et al. also demonstrated that genotypic variation in morphological and physiological responses exists in durum wheat when grown in water limited conditions [1]. Investigating water deficit stress tolerance mechanisms and genotypic differences within a plant species is an important strategy for understanding the basis of stress response and for selection of genotypes with improved water stress tolerance. The genetic mechanism(s) associated with tolerance against abiotic stresses is not well documented in durum wheat, partly because the full genome sequence is still unavailable. Understanding gene regulatory pathways underlying stress responses may lead to new strategies to enhance stress tolerance in durum wheat.
In plants, small non-coding RNAs of 20-24 nucleotides (nts) have been identified as important regulators of genome integrity, virus and pathogen defence, development and importantly, abiotic stress response pathways [2][3][4]. Small RNAs are broadly divided into microRNAs (miR-NAs) and small interfering RNAs (siRNAs). MicroRNAs are global regulators of gene expression mainly through post-transcriptional repression or translational inhibition [5][6][7]. The general molecular networks related to their complex biogenesis and silencing have now been widely characterised [8][9][10][11]. Plant miRNAs control the expression of their targets by binding to imperfect reverse complementary sequences, resulting in degradation and/or translational repression of the cognate target mRNAs [5,11].
Functional analyses of miRNAs and their targets in plants have demonstrated that miRNAs are associated with diverse biological processes including reproductive development and abiotic stress tolerance [12][13][14]. A large number of studies with different plant models have revealed the up-or down-regulation of certain responsive miRNAs when subjected to various abiotic stresses such as water deficit, salinity, heat and cold stress (Table 1). Stress-responsive miRNAs have displayed different regulation patterns between species. However, some stress responsive miRNAs might also exhibit different expression patterns when comparing genotypes of the same plant species; as shown in cowpea exposed to drought stress [15], wheat exposed to dehydration stress [16] and maize exposed to salt stress [17]. Such genotype-specific responses of miRNA help explain the genetic basis of the phenotypic and physiological differences between genotypes of the same species under stress conditions [15,18]. Furthermore, miRNAs have been shown to display spatio-temporal patterns specific to certain plant tissues, suggesting the involvement of tissue-specific miRNAs in various developmental processes [16][17][18]. These tissue-specific patterns have been studied in bread wheat [19,20], but not specifically in durum wheat.
As indicated in Table 1, although numerous miRNAs have been identified in many plant species, including cereals like barley (Hordeum vulgare) [35][36][37], rice (Oryza sativa) [38,39], Brachypodium distachyon [40,41], and bread wheat (Triticum aestivum) [16,19]; only one mature miRNA sequence from Triticum turgidum is recorded in the current miRBase v21. A holistic evaluation of cereal miRNA-mediated response mechanisms under stress conditions is far from complete [42], with very little known about miRNAs and their regulatory functions in relation to water deficit stress across multiple durum genotypes.
This study provides insight into miRNA-mediated water deficit stress regulatory pathways, using four Australian durum genotypes with different water deficit sensitivity [1]. Using Illumina sequencing, we identified 110 conserved miRNAs and 159 novel miRNA hairpin candidates in durum. Statistical analysis has revealed 66 conserved water deficit stress responsive miRNA as well as a number of conserved tissue-and genotype-specific miRNAs. In addition, 16 conserved and five novel miRNA hairpins showed contrasting regulatory patterns under water deficit stress between stress tolerant and sensitive genotypes. To our knowledge, this is the first report of water deficit stress responsive miRNAs identified through direct small RNA cloning and sequencing in durum wheat. Furthermore, target prediction and Gene Ontology (GO) analysis suggests that miRNA targets function in a broad range of biological processes such as metabolic process, response to stimuli, reproduction and development. Comparisons of miRNA profiles in different genotypes under stress in combination with the investigation of target functions and their gene ontologies is a promising approach in predicting miRNA-mediated stress signalling mechanisms in durum wheat, which may have the potential for improving abiotic stress tolerance in breeding programs [42,43].
For novel miRNA identification, a customised bioinformatics approach (Approach #2) was developed. Putative miRNA hairpins were identified using the latest International Wheat Genome Sequencing Consortium's (IWGSC) Chromosomal Survey Sequences (CSS) of bread wheat [44], due to the limited availability of durum wheat sequence. This process resulted in the identification of an initial set of 6,643 loci representing 3,421 non-redundant sequences. Of these non-redundant sequences, 2,710 sequences passed checks by RNAFold and miRcheck, which satisfied in silico requirements of the biogenesis pathway of miRNAs in plants. Of these 2,710 candidate miRNA hairpin sequences, 237 matched the expectations for a true miRNA in terms of their read coverage profile (Category A) using three Boolean metrics as described in the Materials and Methods. Of these, 78 contained an exact match to at least one known mature miRNA from miRBase (Table 2), while the remaining 159 putative novel miRNAs had no match to any known mature miRNAs in the miRBase (S3 Table). Some conserved durum miRNAs are genotype-or tissue-specific regardless of water-deficit stress (Approach #1) Differential miRNA expression profiles were observed between the water deficit stress sensitive (EGA Bellaroi and Tjilkuri) and tolerant (Tamaroi and Yawa) genotypes across both treatments. Comparisons were made between Tamaroi versus EGA Bellaroi, and Yawa versus Tjilkuri, separately, based on their breeding history and genetic background. A total of 70 miRNAs were differentially expressed between different durum genotypes (Fig 1). Among these miR-NAs, four groups displayed interesting expression patterns between the water deficit stress tolerant and the sensitive genotypes (Table 3): I) miRNAs predominantly expressed in water deficit tolerant genotypes under both treatments (7 miRNAs); II) miRNAs predominantly expressed in the water deficit sensitive genotypes under both treatments (5 miRNAs); III) miR-NAs predominantly expressed in the water deficit sensitive genotypes under water deficit stress treatment, but predominantly expressed in the water deficit tolerant genotypes under the control treatment (9 miRNAs); IV) miRNAs predominantly expressed in the water deficit tolerant genotypes under the water deficit stress treatment, but predominantly expressed in the water deficit sensitive genotypes under the control treatment (1 miRNA  (Table 3). In group III, Osa-miR5071 was more abundant in Yawa compared to Tjilkuri in the control treatment libraries (1.78 fold in the flag leaf and 2.10 fold in the developing head, respectively); but was more abundant in EGA Bellaroi compared to Tamaroi in the water deficit treatments (1.78 fold in the flag leaf and 1.56 fold in the developing head, respectively) ( Table 3).
A comparison between all flag leaf and developing head samples identified miRNAs displaying differential abundance between different tissues, irrespective of genotype and treatment. While a total of 110 conserved miRNAs were identified in all sRNA libraries, 86 miRNAs were differentially abundant between flag leaf tissue and the developing head tissue (Fig 2). A total of nine miRNAs were predominantly expressed in the developing head tissue in all four durum genotypes across both treatments while 37 miRNAs were predominantly expressed in the flag leaf tissue (Table 4). For example, Bdi-miR171d was more abundant (from 2.99 to 9.35 fold greater) in the developing head libraries compared to the flag leaf libraries in the four durum genotypes irrespective of the treatment. In contrast, Tae-miR156 was more abundant (from 4.60 to 8.66 fold greater) in the flag leaf libraries compared to the developing head libraries in the four durum genotypes irrespective of the treatment (Table 4).
Water deficit stress-responsive conserved miRNAs in durum (Approach #1) Differential expression of conserved miRNAs were found between water deficit stressed and corresponding control libraries in both the flag leaf and developing head tissues of each durum genotype. Using the criteria described in the Materials and Methods, 66 conserved mature miRNAs were determined to be water deficit stress-responsive miRNAs (Fig 3 and S4 Table).
Hierarchical clustering of the water deficit stress-responsive miRNAs illustrated that several miRNAs showed different regulation patterns under water deficit stress between stress tolerant and sensitive genotypes (Fig 3), whereas certain miRNAs showed the same regulation patterns (e.g. Gma-miR408d was up-regulated under stress of all four durum genotypes in the flag leaf tissues). More interestingly, a small number of stress responsive miRNAs showed up-regulation in water deficit stress sensitive genotypes while those same miRNAs were down-regulated in the tolerant genotypes. For example, in the developing head libraries, Bdi-miR7757 was upregulated in the sensitive genotypes (EGA Bellaroi and Tjilkuri), but was down-regulated in the tolerant genotypes (Tamaroi and Yawa) (Fig 3 and S4 Table). Moreover, some miRNAs responded to water deficit stress only in stress tolerant or sensitive genotypes. In the head libraries, there were 26 miRNAs that were only down-regulated in the stress tolerant genotype Yawa, but not in the stress sensitive genotypes EGA Bellaroi or Tjilkuri (Fig 3). In summary, through further analysing the differentially expressed miRNAs identified through Approach #1, 57 conserved miRNAs were identified as being responsive to water deficit stress, as well as being differentially abundant across different genotypes and tissue types (Fig 4).
Conserved and novel miRNA hairpins showed inverted expression profiles in response to water deficit stress across genotypes (Approach #2) Using the Limma Bioconductor package [45,46], 23 of the 237 putative miRNA hairpins in Category A were found to have a significant tolerance × treatment interaction term. On manual inspection of the miRNA hairpin read-coverage profiles in Category A, 21 of these 23 miRNA hairpins represent strong candidates as they have good read-coverage signatures (Fig 5). Of these 21 candidates, we determined that 16 perfectly matched at least one known mature miRNA in the miRBase, with some hairpins matching to the same conserved miRNA ( Fig 6A). The remaining five novel candidate miRNA hairpins, representing four mature novel miRNAs, do not contain a perfect alignment to any known mature miRNAs ( Fig 6B and S2 Fig). For example, miRNA hairpin Ttu pre-miR008 representing Ttu-miR008 was down-regulated in both flag leaf and developing head tissues under water deficit stress in the stress tolerant genotypes (Tamaroi and Yawa), but was up-regulated in the stress sensitive genotypes (EGA Bellaroi and Tjilkuri).

Validation of differentially expressed miRNAs in durum wheat by quantitative real-time PCR (qPCR)
To validate differentially expressed durum miRNAs predicted by high-throughput sequencing, miRNA was quantified using qPCR. Ten selected stress responsive durum miRNA candidates including seven conserved miRNAs (identical to Ath-miR167d, Gma-miR408d, Bdi-miR5054, Osa-miR5071, Bdi-miR5200, Bdi-miR528 and Zma-miR528a) and three novel miRNAs (Ttu-miR007, Ttu-miR038 and Ttu-miR109) were screened using flag leaf and developing head tissues of four durum genotypes simultaneously. Comparative fold changes of expression levels of miRNA are shown in Fig 7. The expression level changes of conserved miRNAs detected by qPCR were compared with those determined by Illumina sequencing (S5 Table). Most miR-NAs showed similar trends in their expression profile across Illumina sequencing results and qPCR results. For example, in the Illumina sequencing analysis, Zma-miR528a was determined to be down-regulated under stress in the head libraries of Tjilkuri, Tamaroi and Yawa (1.2, 5.1, and 2.7 fold reduction), and up-regulated in EGA Bellaroi (3.6 fold increase). When tested by qPCR, the same miRNA was up/down-regulated in the same libraries and varieties (2.6, 4.9, 1.6 fold reduction, and 2.2 fold increase, respectively). While the expression values between the two platforms are not exactly the same, this has been reported previously and is expected based on the two different quantification methods used [47].   Putative targets of novel water deficit stress responsive durum miRNAs, GO analysis and qPCR To infer the biological functions of the novel water deficit stress responsive miRNAs in durum, putative mRNA target genes were predicted using the psRNAtarget program (http://plantgrn. noble.org/psRNATarget/) with the wheat DFCI gene index (TAGI) version 12 as a reference. A total of 130 targets were identified for four novel stress responsive durum miRNAs (S6 Table). Ttu-miR008 had the highest number of putative target genes (101) while Ttu-miR109 had the lowest (5). On the basis of sequence complementarity between miRNAs and putative target genes, the possible inhibition type between miRNA and their targets was predicted [48,49]. Out of 130 predicted mRNA targets, the inhibition of 109 mRNA targets (83.8%) is caused by cleavage activity, while 21 targets (16.2%) are inhibited through translational repression (S6 Table).
All of the predicted targets were analysed through Gene Ontology (GO) using the Blast2GO server (https://www.blast2go.com/) to further evaluate their putative functions (S6 Table). The BLASTX search obtained the most significant BLAST hits for each target across different species (S3 Fig). According to the ontological definitions of their GO terms, all targets were grouped into three GO categories (S7 Table). At the cellular level (Fig 8A), predicted targets are primarily associated with the nucleus (28.4%), followed by either the mitochondrion or plastid (17.9% each). In evaluating molecular functions, the majority of the targets are potentially involved in either organic or heterocyclic compound binding (16.8% each), ion binding (13.6%), or small molecule binding (10.7%) (Fig 8B). Biologically, nearly half of the targets were classified as being involved in metabolic processes (41.4%) (which includes catabolic, cellular, nitrogen compound, organic substance, primary, and wax metabolic processes) (Fig 8C). The remaining targets were involved in a broad range of biological processes including cellular processes (16.2%), regulation (10.1%), localisation (10.1%), response to stimuli (8.1%), and most significantly, response to stress (5.1%) (Fig 8C). Many of the predicted targets are annotated to be transcription factors, elongation factors, protein phosphatases, and osmotic stress receptors that are associated with multiple stress response processes (S6 Table).
Seven selected targets of Ttu-miR008 were quantified using qPCR (S8 Table). For example, TC438017 (non-specific lipid-transfer protein) and CV779294 (non-specific lipid-transfer protein a-like). In the flag leaf under water stress, TC438017 was up-regulated in the stress tolerant genotypes (4.26 fold in Tamaroi and 2.79 fold in Yawa), whereas it was down-regulated in the stress sensitive genotypes (2.72 fold in EGA Bellaroi and 1.11 fold in Tjilkuri). Similarly, CV779294 was up-regulated in the stress tolerant genotypes (1.34 fold in Tamaroi and 1.40 fold in Yawa), while being down-regulated in the sensitive genotypes (2.37 fold in EGA Bellaroi and 2.41 fold in Tjilkuri). In addition, TC447684 (Glossy 1 protein-GL1) was shown to be upregulated in the developing head of the stress tolerant genotypes (1.22 fold in Tamaroi and 1.13 fold in Yawa), while being down-regulated in the developing head of the sensitive genotypes (1.17 fold in EGA Bellaroi and 1.52 fold in Tjilkuri). Overall, of the seven targets quantified several were negatively correlated with Ttu-miR008, which was down-regulated in the stress tolerant genotypes but up-regulated in the stress sensitive genotypes (Fig 6B).

Discussion
The miRNAome in durum wheat under water deficit stress Water deficit is a major abiotic stress that limits the production of many crops in rain-fed environments. Plant responses to water deficit stress are regulated by complex genetic and epigenetic networks. Interactions between miRNAs and their target mRNAs through sequence-specific binding offer an inheritable and accurate regulation pathway for plants to respond to environmental stimuli at both the translational and post-transcriptional level. To date, extensive efforts have been made to discover water deficit stress-associated miRNAs in many plants including Arabidopsis [24], rice [22], maize [50], soybean [51], barley [52] and bread wheat [16,53]. However, there has rarely been any study on water deficit-stress responsive miRNAs in Triticum turgidum, with only the ssp. dicoccoides being investigated but under shock drought conditions [23].
As an important cereal, mostly grown in rain-fed Mediterranean environments under stressful and variable conditions, durum wheat offers an attractive alternative to studying the much more complex bread wheat genome. With climate change models predicting increased rising temperatures and decreased rainfall, understanding the water deficit stress response pathway(s) in durum wheat has become an important research objective for breeding programs. Using deep sequencing of small RNA libraries in this study, we discovered significant changes that occur with the miRNAome in four durum genotypes under water deficit stress and across two tissue types. Illumina sequencing yielded approximately 623 million reads which were subsequently trimmed and processed to remove inherent redundancy, obtaining a total of 301 million unique sRNA sequences. The highest proportion of the sequenced RNAs was 24 nt in length, which is in agreement with previous studies where 24 nt sRNA fragments constituted the majority of small RNA populations, thereby implicating the function of Dicer proteins during the formation of miRNAs [25,29,54]. Since durum wheat (2n = 4x = 28, genomes AABB) is an ancestral source of the A and B genomes of bread wheat (2n = 6x = 42, genomes AABBDD) and only a partial genome sequence for Triticum turgidum ssp. durum is available, the International Wheat Genome Sequencing Consortium's (IWGSC) Chromosomal Survey Sequences (CSS) of bread wheat was used to identify novel putative miRNA hairpins in durum sRNA libraries [44]. From the 110 conserved miRNAs and 159 novel miRNA hairpins identified, 66 conserved miRNAs and four novel miRNAs were water deficit stress responsive. Further experimental validation including Poly (A)-qPCR and miRNA Ã examination will assist in confirming novel durum miRNA hairpins and their precise excision of the miRNA/miRNA Ã duplex [20,55]. In this study, ten representative stress responsive miRNAs (seven conserved   and three novel) were validated by Poly (A)-qPCR. Poly (A)-qPCR has been shown in bread wheat to provide more accurate and consistent quantification of miRNA expression than stemloop qPCR [56]. Among water deficit stress responsive miRNAs identified in this study, some miRNAs have been found to be associated with abiotic stress response in previous studies; including miR156, miR159, miR167, miR319, miR393, miR398, and miR408. The expression patterns of some of these water deficit stress responsive miRNA were similar to results previously reported. For example, miR159 was up-regulated 1.75 times under water deficit stress in Tjilkuri. Similarly in maize, the expression level of miR159 was significantly increased during drought stress [21]. The up-regulation of miR162, miR167, miR393 under water deficit stress has been commonly observed in different plants ( Table 1), indicating that some miRNA stress-responsive pathways are more than likely to be conserved across different plant species including durum wheat. In contrast, some conserved miRNAs, as well as novel durum miRNAs, were found to be water deficit stress responsive for the first time, including miR1136, miR1432, miR5048, miR5054, miR5071, miR5200 and miR6300. Their regulation pattern indicates that these miRNAs are possibly involved in species-specific response pathways. Most interestingly, the expression profiles of 16 conserved and five novel miRNA hairpins showed inverted regulatory patterns between water deficit stress tolerant and sensitive genotypes, suggesting the regulatory roles of miRNAs in some stress response pathways are genotype-specific (Fig 6). The four durum wheat genotypes used in this study have different levels of water deficit tolerance, which is reflected through their genotypic physiological responses [1]. The distinct genotype differences in miRNA expression profiles could lead to inverted regulation of their functional target genes, which might activate different physiological responses between genotypes [16]. In a recent study of dehydration associated miRNA in wheat, contrasting expression patterns of 13 conserved miRNA (including Tae-miR160a, Tae-miR166h, Tae-miR172a, and Tae-miR393) were also observed between stress tolerant and sensitive genotypes [16]. In the current study, several conserved miRNAs were found to be predominantly expressed in specific genotypes, with or without water deficit stress treatments. For example, miR5200 was consistently more abundant in the water deficit stress sensitive genotypes (EGA Bellaroi and Tjilkuri) than the stress tolerant genotypes (Tamaroi and Yawa) in both the control and water deficit stress libraries. Based on the prediction and further analysis of miRNA targets, we can infer that different capacities for water deficit stress tolerance between durum wheat genotypes may arise from the differential physiological regulation triggered by target genes, which are regulated by genotypic stress responsive miRNAs.

Regulation of miRNA and their targets may contribute to genotypic variation in stress tolerance capacity in different durum genotypes
In the present study, in silico target gene predictions and GO analysis were carried out for four novel water deficit stress responsive miRNAs. This bioinformatics strategy has been applied previously in bread wheat to successfully predict and construct possible miRNA/mRNA target stress regulatory pathways, which were further experimentally validated [16,19,53,[57][58][59]. A total of 130 target genes for four novel durum miRNAs were predicted to encode proteins of diverse functions. GO analysis indicated that these targets are involved in a broad range of biological processes and varied physiological responses in durum wheat, such as biosynthetic activity, binding activities with proteins and nucleic acids, protein transport, abscisic acid  (ABA) metabolic processes, photosynthetic activity and leaf senescence. Significantly, stress responsive expression of seven predicted target genes were validated by qPCR. The negative correlation of several targets with their corresponding miRNA implies the involvement of miRNA-mRNA target regulation in stress response pathways in durum.
A significant number of targets are predicted to possess nucleic acid binding activities and encode transcription factors involved in signalling and defence, which contribute to stress tolerance in different durum genotypes. For example, auxin response factor (ARF) 18-like is a target of Ttu-miR008. ARFs bind to auxin response elements to usually negatively regulate expression of auxin-inducible genes such as GH3 (Gretchen Hagen3), Aux/IAA (auxin/indole-3-acetic acid) and SAUR (small auxin-up RNA) [60]. Several auxin-responsive genes have been identified to respond to various abiotic stress conditions such as drought, salinity and cold in Arabidopsis, rice and sorghum, indicating the cross-talk between auxin signalling and abiotic stress responses [61][62][63]. In durum, Ttu-miR008 is down-regulated under stress in the tolerant genotypes suggesting that ARF18-like protein increases thereby repressing auxin-inducible genes enhancing auxin signalling. This might affect processes which require a lower auxin:cytokinin ratio, such as lateral root development [64]. In maize and wheat, the development of lateral roots in the stress tolerant genotype is enhanced from the accumulation of auxinresponsive factors [16,17]. However, the role of miRNA and ARF in lateral root development in durum needs to be confirmed with further experimental validation.
Other targets also contribute to water stress tolerance in durum as signalling factors including protein kinases and protein phosphatases. For example, a target of Ttu-miR008 (TC451175) was annotated as a probable protein phosphatase 2C (PP2C). Studies in Arabidopsis and rice have shown that PP2C genes were induced by diverse environmental stimuli and acted as positive regulators in ABA-mediated signalling pathways well known to be involved in stress responses [65,66].
However, there are also other targets of Ttu-miR008 which could contribute to water deficit stress tolerance in different ways such as maintaining osmotic pressure of the plant or homeostasis of the cell. For example, the target CV769573 identified in this study as an ABA 8'hydroxylase, is a key enzyme in ABA degradation [67]. ABA is crucial for various stress responses, including regulation of stress-responsive genes, stomatal closure, and metabolic changes [68]. ABA is rapidly increased in response to environmental stress [67], suggesting a role for removing ABA 8'-hydroxylation to ensure increased ABA levels. Equally rapid elimination of stress induced ABA when stresses are relieved is essential [69]. Indeed, dehydration stress can cause steady increases in ABA degradation in Arabidopsis over time [70]. Although requiring confirmation, ABA 8'-hydroxylase may therefore decrease to a lesser extent in tolerant genotypes suggesting they have a lower ABA requirement during water deficit stress.
Also identified and quantified in this study was the Glossy 1 (GL1) protein, which is yet another target of Ttu-miR008 (TC447684). GL1 functions in the biosynthesis pathway of cuticular wax, which provides protection against environmental stress. In rice, OsGL1 over-expression plants showed increased cuticular wax accumulation on the leaf surface and were more tolerant to drought stress at reproductive stages compared to the wild type [71]. The inhibition of GL1 is reduced through the down-regulation of stress responsive miRNA, leading to enhanced wax production, thus preventing water loss. This helps to explain the genotypic difference in the reduction of relative water content in leaves, in response to water deficit stress between stress tolerant and sensitive durum genotypes [1].
Two other quantified functional targets, TC438017 (non-specific lipid-transfer protein) and CV779294 (non-specific lipid-transfer protein a-like), examined by qPCR may also assist to explain the genotypic difference in maintaining osmotic pressure. Lipid transfer proteins (LTPs) help to repair stress-induced damage in membranes or alter the lipid composition of membranes. In pepper, the accumulation of LTP transcripts induced by environmental stresses is associated with cuticle formation, which contributes to the avoidance/tolerance of low tissue water potential and water content [72,73]. In this study, TC438017 and CV779294 were negatively correlated with their corresponding miRNA showing genotypic expression patterns in response to water deficit. The up-regulated accumulation of LTPs observed only in stress tolerant durum genotypes helps to explain the genotypic differences in the maintenance of leaf water potential and relative water content [1], suggesting the participation of miRNA/target interaction in genotypic physiological response pathways in durum. Experimental examination of these miRNA-regulated targets also helps demonstrate the validity of prediction analysis using bioinformatics.

Conclusion
The present study provides a comparative description of the miRNAome in durum wheat between water deficit tolerant and sensitive genotypes in response to water deficit stress, suggesting that there are multiple miRNA regulation patterns which might contribute to, and partly explain, the distinct water deficit stress sensitivities between different durum genotypes. The first comprehensive durum small RNA dataset generated provides a good foundation for future characterisation of the molecular mechanisms underlying water deficit stress tolerance in durum. This was achieved through Illumina sequencing, which enabled profiling of the miR-NAome in water deficit stress tolerant and sensitive durum wheat genotypes across different tissues and treatments. We have identified 110 conserved miRNAs and 159 novel miRNA hairpins in durum wheat, including 66 conserved miRNAs and five novel miRNA hairpins (representing four novel miRNAs) that are water deficit stress responsive. A total of 16 conserved miRNA hairpins (representing 11 conserved miRNAs) and five novel miRNA hairpins (representing four novel miRNAs) showed distinct down-regulation profiles in the water deficit stress tolerant genotypes while the same miRNAs were up-regulated in sensitive genotypes. This demonstrates that regulation patterns of the same miRNAs may vary extensively across genotypes of the same species, in response to environmental stimuli. Target prediction and GO analysis of four novel genotype-specific regulated miRNAs provide evidence for the potential involvement of miRNAs in a broad range of biological processes, including stress response pathways. Several potentially valuable target genes have been identified and are now undergoing further experimental validation, which will be reported elsewhere.

Plant material and growth conditions
Four durum wheat genotypes (EGA Bellaroi, Tamaroi, Tjilkuri and Yawa) were used in this study. Seeds were obtained from Durum Breeding Australia's (DBA) southern node breeding program (The University of Adelaide). Tamaroi and Yawa are water deficit stress tolerant genotypes; while EGA Bellaroi and Tjilkuri are water deficit stress sensitive genotypes [1]. Plants were grown at 22°C/12°C day/night temperature with a 12 h photoperiod with watering to field capacity (12% soil water content (SWC)) from germination to booting stage when the water limiting stress treatment was imposed for 15 d (6% SWC or 50% field capacity; water deficit stress group, WG) or field capacity maintained (control, CG), as per Liu et al. [1].

Sampling and total RNA extraction
After 15 d of water deficit stress, the flag leaf and developing head were collected with sterile razor blades and frozen immediately in liquid nitrogen. Frozen tissues were ground to a fine powder in liquid nitrogen using a sterile mortar and pestle, pre-chilled to -80°C. Total RNA was isolated using the TriPure isolation reagent kit (Roche Diagnostics, Australia) and treated with RQ1 RNase-Free DNase I (Promega, Australia) following the manufacturer's instructions.
The concentration and quality of extracted RNA samples were measured by spectrophotometric analysis at 260 nm and 280 nm using a NanoDrop Lite spectrophotometer (Thermo Scientific, USA). RNA integrity was assessed by agarose gel electrophoresis. A total of 96 RNA samples (4 durum genotypes × 2 tissue types × 2 treatment groups × 6 biological replicates = 96) were extracted and stored at -80°C for downstream applications.

Small RNA library construction and deep sequencing
For small RNA library construction, 5 μg of total RNA was size-fractionated on a 15% denaturing TBE urea polyacrylamide gel and small RNAs (15 to 40 nt) were excised using an NEB miRNA marker (New England Biolabs, UK) as a guide. Small RNAs was eluted in 0.3 M NaCl by rotating the tube overnight at 4°C. Eluted RNA was passed through a Spin-X column and then precipitated using glycoblue (Ambion, USA) and isopropanol. The sRNA pellets were washed and air-dried at room temperature, then re-suspended in DEPC-treated water [74]. A total of 96 small RNA libraries were constructed from flag leaf and developing head of durum wheat plants that had been treated or not treated with water deficit stress (4 durum genotypes × 2 tissue types × 2 treatment groups × 6 biological replicates = 96) using NEB Next 1 Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, UK) following the manufacturer's instructions. For each flag leaf sRNA library and head sRNA library, a unique index primer was used for multiplexing purposes using the NEBNext 1 Index Primer Set (New England Biolabs, UK). The final cDNA product was purified using Pippin Prep™ System (Sage Science, USA). Prior to sequencing, quality and quantity of the amplified small RNA cDNA libraries was evaluated on an Agilent 2100 Bioanalyzer system (Agilent Technologies, USA) and Qubit fluorometer (Invitrogen, USA). All 96 small RNA libraries were sequenced using Illumina sequencing technology on a HiSeq2500 machine after cluster generation. All sequencing reads were submitted to the NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/), and are accessible under the accession number GSE69339.

Identification of conserved miRNAs (Approach #1)
In this study, Approach #1 was developed to identify conserved miRNAs in durum wheat using CLC Genomics Workbench v7.0 (CLC Bio, Denmark). Briefly, raw sequencing reads were first processed by trimming adaptor sequences and removing low-quality reads. Sequences shorter than 15 nt and larger than 50 nt were excluded from further analysis. Trimmed reads were generated for each small RNA library and then annotated to determine the presence of known plant miRNAs. Durum small RNA sequences were aligned with known miRNAs in miRBase using CLC Genomics Workbench v7.0 based on their sequence homology, allowing up to two mismatches in alignment [15]. Conserved miRNAs in common monocot and dicot species (Triticum aestivum, Triticum turgidum, Brachypodium distachyon, Zea mays, Oryza sativa, Hordeum vulgare, Sorghum bicolor, Arabidopsis thaliana, and Glycine max) deposited at miRBase v20 (June 2013) were used as references for annotation.
Normalisation of miRNA abundance in each library was carried out using a value referred to as RPM (reads per million). The RPM value was obtained by dividing the reads number of a miRNA with the total number of putative sRNA reads in each library and multiplying by a million. Matched sequences with no more than two mismatches and with an abundance of over two RPM in at least 50% of the 96 libraries were considered as candidate conserved miRNAs.

Identification of differentially expressed conserved miRNAs (Approach #1)
Differentially expressed conserved miRNAs were identified based on the RPM. To identify differentially expressed miRNAs, the following criteria were used: 1) number of miRNA reads was set as 0.01 by default when the sequencing read was 0; 2) normalised reads (RPM) was at least 10 in one of the libraries in comparison; and 3) the fold-change of normalised reads of libraries in comparison was greater than 1.5 [16,75]. For expression analysis, reads of unique mature miRNAs deposited in miRBase were used as they are an active and functional form of mature miRNAs [29]. Tissue-specific conserved miRNAs were identified by comparing flag leaf libraries with head libraries. Genotype-specific conserved miRNAs were identified by comparing water deficit sensitive varieties and water deficit tolerant varieties. Comparisons were made only between EGA Bellaroi and Tamaroi, or Tjilkuri and Yawa due to their breeding background. Water deficit stress-responsive miRNAs were identified by comparing control treatment libraries with water deficit stress treatment libraries. Heat maps of differentially expressed miRNAs were generated in R (version 3.1.2) (http://www.r-project.org/). Where the fold change of some conserved miRNA candidates were not analysed due to their low reads in the sequencing results (RPM were less than 10 in both libraries for differential expression comparison), their log2 fold change under stress was recorded as zero in the clustering analysis.
Small RNA-Seq data pre-processing for novel miRNA identification (Approach #2) To identify novel miRNAs in durum wheat, a customised bioinformatics approach (Approach #2) was developed. Small RNA-Seq raw reads were 5' and 3' adapter trimmed and the output partitioned into two sets of reads: 1) those that had been trimmed and were 19-26 bp long, and 2) those that did not contain any adapter sequence. The first set represents non-redundant (NR) 3' adapter trimmed reads, which were used to identify putative pre-miRNA hairpin. In order to remove reads which are derived from the breakdown products of longer mRNA's rather than true sRNA molecules, the second set of reads and the NR sRNA reads were de novo assembled to generate a reference against which sRNA reads would be filtered. This was done using Velvet (v 1.2.09, https://www.ebi.ac.uk/~zerbino/velvet/) with a kmer length of 17 and read tracking enabled. The NR set of 3' adapter trimmed reads (19-26 bp) were filtered to remove those which were either: a) low abundance (< = 5 reads in all samples); b) mapped to known wheat rRNAs; c) mapped to the wheat chloroplast or mitochondrial genomes; d) mapped to the 50bp+ long de novo assembled contigs; e) mapped to UniVec (build 7.1, http:// www.ncbi.nlm.nih.gov/VecScreen/UniVec.html) data set; or f) mapped to the Triticeae Repeat Sequence Database (TREP) database of grass repeat sequences [76]. The mappings in steps b-f above was performed using Bowtie2 (v 2.2.3; http://bowtie-bio.sourceforge.net/bowtie2/index. shtml) with parameters which allowed up to two mismatches and no indels. The NR reads which passed all the above filters were used to identify candidate pre-miRNA hairpins.

Identification of miRNA precursors and novel miRNA candidates in durum wheat (Approach #2)
Since only a partial genome sequence for Triticum turgidum ssp. durum is available, the International Wheat Genome Sequencing Consortium's (IWGSC) Chromosomal Survey Sequences (CSS) [44] was used to identify putative miRNAs. The NR sRNA sequences which passed the filters were mapped to the IWGSC CSS using BioKanga v3.4.3 (http://sourceforge.net/projects/ biokanga/) in order to identify all possible contigs from which the sRNA sequence could have been derived. For each NR 3' adapter trimmed read, all perfect alignment locations in the IWGSC CSS were identified. Using a subset of reads and CSS contigs involved in those perfect alignments, we also identified all imperfect alignments (two-five mismatches). The candidate pre-miRNA hairpins were defined using all pairwise combinations of perfect to imperfect alignments of a given read within a CSS contig. Additional constraints were applied such that the perfect and imperfect alignments were in opposite orientations and separated by 54-1000 bp. A NR set of these regions ±20 bp, were processed by RNAFold (http://rna.tbi.univie.ac.at/ cgi-bin/RNAfold.cgi) and then miRcheck (http://web.wi.mit.edu/bartel/pub/software.html) to ascertain if they could form hairpin structures with characteristics associated with the miRNA biogenesis pathway in plants, indicating the formation of a miRNA/miRNA Ã duplex from stem-loop hairpins based on their read coverage profile [55]. Three primary criteria were applied as follows: 1) A peak of reads in the first or last 50 bp of the hairpin sequence all aligned to the same strand/stem (the miRNA site); 2) a second, smaller peak of complementary reads aligned on the opposite end to the miRNA strand/stem (the miRNA Ã site); 3) a small proportion of reads mapping between the above two defined regions (the loop). All candidate miRNA hairpin sequences were classified into one of eight categories (A-H, where A has a read coverage profile matching the expectations for a true miRNA) using three Boolean metrics based on their read coverage profile: 1) if !95% of the reads mapped to one strand of the hairpin; 2) if !95% of the reads mapped to one of the terminal 50 bp of the hairpin; and; 3) if 5% of the reads mapped to the loop region of the hairpin (Table 2). Putative miRNA hairpins were further characterised by identifying if their sequence contained any perfect matches to the 35,828 mature miRNAs from miRBase v21 (accessed July 2014).

Identification of stress responsive novel miRNA hairpins (Approach #2)
To identify novel water deficit stress-responsive miRNA hairpins, the Limma Bioconductor (v3.18.13) package [45] was used to perform a statistical analysis using linear models based on the RPM data. Different durum varieties were recoded with binary values which indicated water deficit stress tolerance or water deficit stress sensitivity. Of the many possible contrasts that could be made, the tolerance × treatment interaction term was of primary interest in the linear model. This effectively identified hairpins which showed differential expression to water deficit stress and that this response was different between water deficit stress sensitive and water deficit stress tolerant cultivars. Pre-miRNA hairpins from Category A ( Table 2) which had a significant tolerance × treatment interaction were then inspected to ascertain if their read-coverage profiles followed what we expected from a true mature miRNA and miRNA Ã .

Quantitative real-time PCR (qPCR) of miRNA candidates
In order to evaluate the expression of miRNA candidates, poly-A tailing combined with qPCR was performed for a select group of seven conserved and three novel stress responsive miRNAs with the 96 durum total RNA samples which were used for sRNA library construction. For each sample, 1 μg of total RNA was poly-A tailed and reverse-transcribed with the NCode VILO miRNA cDNA synthesis kit (Invitrogen, USA) following the manufacturer's instructions. The final cDNA product was diluted to 100 μL. qPCR was performed using the ViiA™ 7 Real-Time PCR system (Applied Biosystems, USA). In each 10 μL qPCR reaction (six biological replicates for each sample), 1 μL diluted cDNA template and primers (3 pmol of each forward and reverse) were mixed with SYBR 1 Green reagent (iQ TM supermix, BioRad, USA). The forward miRNA primers were designed based on the full mature miRNA sequences (S9 Table). The reverse primer was the universal reverse primer provided in the NCode VILO miRNA cDNA synthesis kit. The qPCR running conditions were: 95°C for 2 min, followed by 40 cycles of 95°C for 15 s, 56/58/60°C for 15 s, and 70°C for 10 s, followed by 72°C for 10 min. Melting curve analysis was used to detect the specificity of the amplified product. The relative expression ratio was calculated using the comparative CT ( ΔΔ C T ) method with GAPDH [Gen-Bank: AF251217] as the reference gene.

Target prediction, functional GO analysis and target qPCR
The putative mRNA targets of stress responsive novel miRNAs were identified using psRNA Target Server (http://plantgrn.noble.org/psRNATarget/) with the following parameters: prediction score cut-off value = 3.0, length for complementarity scoring = 20, and target accessibility = 25. Mature novel miRNA sequences were used as queries and the wheat DFCI gene index (TAGI) version 12 was used as the reference genome dataset [19]. All the predicted targets were evaluated using the functional enrichment analysis tool at Blast2GO (http://www. blast2go.com) [77,78]. BLASTX was employed to perform a homology search against the NR protein databases in NCBI to obtain the most significant BLAST hits for each target using the Blast function with Blast2GO. Default parameters were used in the mapping and annotation steps to obtain GO terms for each target transcript in Blast2GO. The annotation results were further improved by analysing conserved domains/families using the InterProScan function. GO terms for three categories (cellular component, molecular function and biological processes) were determined for each annotated target. All the annotated targets were classified on the basis of their GO term enrichments in each category. Seven selected functional targets were quantified using qPCR with the same cDNA libraries employed in the miRNA qPCR. Target qPCR was performed using the comparative CT ( ΔΔ C T ) method with GAPDH as the reference gene [GenBank: AF251217]. Target primers were designed to include the predicted miRNA/ mRNA binding region in the amplified product ensuring the quantification of uncleaved targets, in order to examine the correlation of miRNA and regulated targets. Target transcript sequences, primer locations and primer sequences are listed in S10 Table. Supporting Information  Table. Fold-change of selected water deficit stress responsive miRNA candidates identified by Illumina sequencing and qPCR. Fold changes have been determined by comparing the RPM in Illumina sequencing or comparing relative expression ratio in qPCR between the control treatment and the water deficit stress treatment in different tissues of four durum wheat genotypes. Bold fold change value indicates that the miRNA was more abundant in the water deficit stress treatment libraries whereas unbolded fold change indicates that the miRNA was more abundant in the control treatment libraries.  Table. Predicted targets of four novel durum stress responsive miRNAs and their GO analysis results. Definitions: Column E (Expectation)-The expectation scoring of the complementarity between miRNAs and their targets. The maximum expectation threshold score was set at 3.0. Column F (Target Accessibility (UPE))-The maximum energy required to open (unpair) the secondary structure around the target site on the target mRNA. Column O (Multiplicity)-Multiplicity of the target site representing the number of target sites within a specific target transcript. (XLS) S7 Table. Combined Gene Ontology classification in GO levels of 130 predicted targets of four novel miRNAs. Definitions: Column A (Level)-The GO level represents the position of a GO term in the GO hierarchy. The level of a GO term is the number of GO terms between that term and the Root Term of the Ontology. Column E (Node score)-The node score is the sum of sequences directly or indirectly associated to a given GO term weighted by the distance of this term to the term of its direct annotation, i.e. the GO term the sequence is originally annotated to. This confluence score takes into account the number of sequences converging at one GO term and at the same time penalises by the distance to the term where each sequence was actually annotated. Column F (%Seq)-The percentage of sequences annotated with a particular GO term among all the sequences annotated within the same GO level. Column G (#Seq)-The number of target sequences annotated with that particular GO term. (XLS) S8 Table. Fold-change of seven selected functional targets of Ttu-miR008 quantified by qPCR. Green values indicate that the targets were up-regulated under water deficit stress, while red values indicate that the targets were down-regulated under water deficit stress. Bold fold change values indicate negative correlation with Ttu-miR008. FL = Flag leaf libraries; H = Head libraries. (XLS) S9 Table. Forward primers used in qPCR validation of seven conserved and three novel miRNAs in durum. Each forward primer was designed based on the full sequence of the mature miRNA. (XLS) S10 Table. Target transcript sequences, primer locations and primer sequences used in qPCR validation of seven selected target genes. (XLS)