Transgenerational Effects of Water-Deficit and Heat Stress on Germination and Seedling Vigour—New Insights from Durum Wheat microRNAs

Water deficiency and heat stress can severely limit crop production and quality. Stress imposed on the parents during reproduction could have transgenerational effects on their progeny. Seeds with different origins can vary significantly in their germination and early growth. Here, we investigated how water-deficit and heat stress on parental durum wheat plants affected seedling establishment of the subsequent generation. One stress-tolerant and one stress-sensitive Australian durum genotype were used. Seeds were collected from parents with or without exposure to stress during reproduction. Generally, stress on the previous generation negatively affected seed germination and seedling vigour, but to a lesser extent in the tolerant variety. Small RNA sequencing utilising the new durum genome assembly revealed significant differences in microRNA (miRNA) expression in the two genotypes. A bioinformatics approach was used to identify multiple miRNA targets which have critical molecular functions in stress adaptation and plant development and could therefore contribute to the phenotypic differences observed. Our data provide the first confirmation of the transgenerational effects of reproductive-stage stress on germination and seedling establishment in durum wheat. New insights gained on the epigenetic level indicate that durum miRNAs could be key factors in optimising seed vigour for breeding superior germplasm and/or varieties.


Introduction
Global food demand is constantly increasing with the rise in our world population. However, even with a finite amount of arable land, agricultural production is still expected to provide more food and fibre, as well as to adapt to the environmental challenges (i.e., climate resilience). High quality seeds or grains are therefore essential for sustainable agriculture and food security. Cereal crop production starts with seed germination to establish vigorous healthy seedlings in the field. Successful germination and robust seedling growth are critical early stages of crop life to ensure reliable canopy establishment and production efficiency [1,2]. Seeds with high vigour are expected to germinate rapidly and uniformly, thus warranting field emergence with minimum delay after sowing [3]. Seedlings with high vigour are expected to be more resource-use-efficient and to have competitive advantages over unwanted weed species [1,2].
For cereal crops, successful germination and vigorous seedling emergence relies on high quality seeds. Durum wheat (Triticum turgidum L. ssp. durum) is a vital cereal crop that provides substantial economic output in the Mediterranean region [4] due to its unique use in multiple food products such as pasta and couscous. Durum is better adapted than common hexaploid wheat to semi-arid environments, Plants 2020, 9,189 3 of 20 water-deficit and heat stress group), L6CG (seeds from L6 control group), and L6WH (seeds from L6 water-deficit and heat stress group).
For seed germination measurements (Table 1), genotype had a significant impact (p < 0.05) on the germination potential (Gp), germination rate (Gr), and germination index (Gi) without interaction with the parent treatment factor. Overall, DBA Aurora had better germination performance than L6. In particular, AuCG and AuWH almost had double the germination potential of L6CG and L6WH. Parent treatment only had a significant impact (p < 0.05) on Gp and Gi without interaction with the genotype. Overall, water-deficit and heat stress treatment on the parents negatively affected the germination performance of the next generation. Coleoptile length was not significantly affected by any of the factors. Interestingly, a significant interaction (p = 0.005) was detected between genotype and parent treatment for young root length. For DBA Aurora, the parent stress treatment did not significantly affect the young root of the progenies. However, for L6, stress on the parents significantly reduced the length of the young roots. Table 1. Germination potential (Gp, %), germination rate (Gr, %), germination index (Gi), coleoptile length, and young root length of the four seed groups. AuCG, seeds from DBA Aurora control group. AuWH, seeds from DBA Aurora stress group. L6CG, seeds from L6 control group. L6WH, seeds from L6 stress group. 3.1% n.a 1.6 n.a n.a l.s.d Genotype × Parent treatment n.a n.a n.a n.a 3.38 For the seedling vigour index I (SVI-I) and seedling vigour index II (SVI-II) ( Table 2), both genotype and parent treatment had a significant impact (p < 0.05) without interaction with each other. Overall DBA Aurora outperformed L6. Water-deficit and heat stress (WH) treatment of the parents negatively impacted both genotypes, but to a greater extent in L6. For other seedling measurements (Table 2), genotype had a significant impact (p < 0.05) on shoot fresh weight, shoot dry weight, and seedling height, with no interaction with the parent treatment factor. Overall, DBA Aurora again displayed a better performance for these traits when compared to L6. Parent treatment also had a significant impact (p < 0.05) on shoot fresh weight and seedling height without interaction with the genotype factor. The WH treatment of the parents significantly reduced the shoot fresh weight and seedling height of the next generation (p < 0.05). Interestingly, a significant interaction (p < 0.05) was detected between genotype and parent treatment for all root traits (root fresh weight, root dry weight and root length). For DBA Aurora, the WH treatment of the parents had no impact on the root fresh weight and dry weight of the next generation seedlings. For L6, the WH treatment of the parents significantly reduced the root fresh weight and dry weight. The WH treatment of the parents promoted the root length in both genotypes, but to a greater extent in DBA Aurora.

The microRNAome in the Progenies
Four sRNA libraries (TCG, TWH, SCG and SWH) were constructed and sequenced from the four seedling groups (AuCG, AuWH, L6CG and L6WH), respectively. A total of 83,991,830 raw reads were obtained (Table 3). Approximately 52.04 million clean sRNA reads within the 18-25 nt range remained after data processing. Among these, 26.16 million non-redundant unique sRNA reads were obtained. The size distribution of sRNA reads exhibited a two-peak pattern (a weak peak at 21 nt and a strong peak at 24 nt) (Figure 1), a classic profile expected for cleaved products of plant Dicer-like enzymes [20,22,34].

The microRNAome in the Progenies
Four sRNA libraries (TCG, TWH, SCG and SWH) were constructed and sequenced from the four seedling groups (AuCG, AuWH, L6CG and L6WH), respectively. A total of 83,991,830 raw reads were obtained (Table 3). Approximately 52.04 million clean sRNA reads within the 18-25 nt range remained after data processing. Among these, 26.16 million non-redundant unique sRNA reads were obtained. The size distribution of sRNA reads exhibited a two-peak pattern (a weak peak at 21 nt and a strong peak at 24 nt) (Figure 1), a classic profile expected for cleaved products of plant Dicer-like enzymes [20,22,34]. A total of 2202 MIR-miRNA entries (different combinations of MIR loci origin and mature miRNA products) were identified (Table S1). All MIR-miRNA entries were categorised into five groups (G1-G5), with G1-G4 representing conserved miRNAs and G5 representing novel miRNAs A total of 2202 MIR-miRNA entries (different combinations of MIR loci origin and mature miRNA products) were identified (Table S1). All MIR-miRNA entries were categorised into five groups (G1-G5), with G1-G4 representing conserved miRNAs and G5 representing novel miRNAs ( Table 4). The sequence reads that mapped to the durum genome and ESTs (expressed sequence tags) but their extended genome sequences could not form secondary hairpins (G3) had the highest number of unique miRNAs. The group that had neither aligned reference pre-miRNAs nor sequencing reads that could be further mapped to the genome (G4) had the lowest number of unique miRNAs (Table 4). In Table  S1, MIR-miRNA entries with the same pre-miRNA index number within each group share the same pre-miRNA sequence. As one pre-miRNA can potentially produce multiple unique mature miRNAs, each unique miRNA is labelled with its own miRNA index number within each group. The pre-miRNA and miRNA index numbers in Table S1 can be linked to Table 4. For Group 2 to Group 4, the total number of pre-miRNAs and the total number of unique mature miRNAs correspond to the highest pre-miRNA index number and the highest miRNA index number within each group in Table S1 (e.g., the last MIR-miRNA entry for Group 2 in Table S1 is indexed as 446-483). However, the highest miRNA index number in Group 1 (Table S1) is different to the total number of unique miRNAs in Group 1 ( Table 4) as the same mature miRNA (unique miRNA) can be produced from different pre-miRNAs with different index numbers. Table 4. Categorisation of all MIR-miRNA entries identified in this study. G1-G4 represent conserved miRNAs while G5 represents novel miRNAs. Unique miRNAs represent distinct mature miRNA products. G1, pre-miRNAs can be mapped to the durum genome and ESTs (expressed sequence tags). G2, sequencing reads mapped to the durum genome and ESTs with their extended genome sequences forming secondary hairpins. G3, sequencing reads mapped to the durum genome and ESTs, but their extended genome sequences do not form secondary hairpins. G4, neither aligned reference pre-miRNAs nor sequencing reads can be further mapped to the genome. G5, novel miRNAs where reads cannot be mapped to the miRBase but can be mapped to the genome, and secondary hairpins can be formed from extended genome locations. TCG, constructed from DBA Aurora CG seedlings. TWH, constructed from DBA Aurora WH seedlings. SCG, constructed from L6 CG seedlings. SWH, constructed from L6 WH seedlings.  G1  161  270  158  252  158  256  158  259  158  257  G2  446  483  368  392  372  405  382  410  300  322  G3  838  890  663  694  701  736  722  763  611  644  G4  48  49  34  34  31  33  42  43  34  34  G5  437  485  426  446  429  445  429  455  421  440 Among all entries (Table S1), 200 miRNAs were highly expressed (normalised reads number of the reported miRNA is higher than the average of the dataset), including 190 conserved miRNAs and 10 novel miRNAs. A total of 1435 miRNAs showed a medium expression level (normalised reads number of the reported miRNA is higher than 10 but less than the average of the dataset), including 1002 conserved miRNAs and 433 novel miRNAs. A total of 567 miRNA entries were lowly expressed (normalised reads number of the reported miRNA is less than 10), including 525 conserved miRNAs 42 novel miRNAs. In total, these 2202 MIR-miRNA entries represent 2177 unique mature miRNA products among four libraries, as shown in the Venn diagram ( Figure 2). Collectively, the conserved miRNAs identified in this study represent 70 MIR families (Table S2). The conservation profiles of all miRNAs across the reference plant species (Table S3) in miRBase are also shown ( Figure 3).

Differentially Expressed miRNAs (DEMs)
Differentially expressed miRNAs (DEMs) were identified based on their normalised reads number. A total of 1368 miRNAs (Table S4) showed differential expression across four libraries (TCG vs TWH vs SCG vs SWH). A heatmap was compiled to show DEMs across four libraries using their log10 (normalised count) value ( Figure S1). Some miRNAs exhibited a genotype-specific pattern (Table S5). A total of 144 miRNA were found to be exclusively expressed in DBA Aurora (TCG and/or TWH), while 139 miRNAs were only expressed in L6 (SCG and/or SWH). A number of miRNAs also showed Plants 2020, 9, 189 7 of 20 specificity to parent treatment (Table S5). A total of 129 miRNAs were only detected in the CG libraries (TCG and/or SCG) while 72 miRNAs were only found in the WH libraries (TWH and/or SWH).

Differentially Expressed miRNAs (DEMs)
Differentially expressed miRNAs (DEMs) were identified based on their normalised reads number. A total of 1368 miRNAs (Table S4) showed differential expression across four libraries (TCG vs TWH vs SCG vs SWH). A heatmap was compiled to show DEMs across four libraries using their log10 (normalised count) value ( Figure S1). Some miRNAs exhibited a genotype-specific pattern (Table S5). A total of 144 miRNA were found to be exclusively expressed in DBA Aurora (TCG and/or TWH), while 139 miRNAs were only expressed in L6 (SCG and/or SWH). A number of miRNAs also showed specificity to parent treatment (Table S5). A total of 129 miRNAs were only detected in the CG libraries (TCG and/or SCG) while 72 miRNAs were only found in the WH libraries (TWH and/or SWH).

Differentially Expressed miRNAs (DEMs)
Differentially expressed miRNAs (DEMs) were identified based on their normalised reads number. A total of 1368 miRNAs (Table S4) showed differential expression across four libraries (TCG vs TWH vs SCG vs SWH). A heatmap was compiled to show DEMs across four libraries using their log10 (normalised count) value ( Figure S1). Some miRNAs exhibited a genotype-specific pattern (Table S5). A total of 144 miRNA were found to be exclusively expressed in DBA Aurora (TCG and/or TWH), while 139 miRNAs were only expressed in L6 (SCG and/or SWH). A number of miRNAs also showed specificity to parent treatment (Table S5). A total of 129 miRNAs were only detected in the CG libraries (TCG and/or SCG) while 72 miRNAs were only found in the WH libraries (TWH and/or SWH).  When comparisons were made between two libraries within the genotype or within the treatment source (TWH vs TCG, SWH vs SCG, SCG vs TCG, SWH vs TWH), DEMs were identified using Chi-squared 2 × 2 test and Fisher's test (see Materials and Methods). The number of up-and down-regulated DEMs identified in each comparison using different p values as thresholds are shown in Figure 4. Interestingly, a genotypic pattern can be observed here for the number of DEMs. For the stress-tolerant genotype DBA Aurora (TWH vs TCG), the number of miRNAs that were up-regulated was higher than the ones that were down-regulated, irrespective of the p value used. However, for the stress-sensitive genotype L6 (SWH vs SCG), an opposite trend was observed, where there were more miRNAs that were down-regulated than up-regulated. When the same parent treatment source Plants 2020, 9, 189 8 of 20 was compared between two genotypes (SCG vs TCG, SWH vs TWH), there are more down-regulated miRNAs than up-regulated miRNAs (except for SCG vs TCG at p < 0.1). The expression levels of all DEMs are shown in Table S6. Highly expressed DEMs with a |log2 (fold change)| > 1 in each comparison (p < 0.05) are listed in Table 5.
using Chi-squared 2 × 2 test and Fisher's test (see Materials and Methods). The number of up-and down-regulated DEMs identified in each comparison using different P values as thresholds are shown in Figure 4. Interestingly, a genotypic pattern can be observed here for the number of DEMs. For the stress-tolerant genotype DBA Aurora (TWH vs TCG), the number of miRNAs that were upregulated was higher than the ones that were down-regulated, irrespective of the P value used. However, for the stress-sensitive genotype L6 (SWH vs SCG), an opposite trend was observed, where there were more miRNAs that were down-regulated than up-regulated. When the same parent treatment source was compared between two genotypes (SCG vs TCG, SWH vs TWH), there are more down-regulated miRNAs than up-regulated miRNAs (except for SCG vs TCG at P < 0.1). The expression levels of all DEMs are shown in Table S6. Highly expressed DEMs with a |log2 (fold change)| > 1 in each comparison (P < 0.05) are listed in Table 5.

Target Prediction and Enrichment Analysis
The target genes were predicted for highly expressed DEMs using in silico analysis with the recently constructed durum wheat genome (see Materials and Methods). A total of 311,964 transcripts (target genes) were identified as potential targets for 189 DEMs (Table S7). Gene ontology (GO) analysis was conducted to classify the biochemical functions of the predicted target genes. GO annotation grouped all targets in the following three GO categories (i) Biological Process (BP), (ii) Molecular Function (MF), and (iii) Cellular Component (CC). GO enrichment analysis of the predicted targets was performed (Table S8). The most represented GO terms are shown in Figure 5. Notably, a few biological processes with importance are highly represented, including the regulation of transcription, DNA-templating (GO:0006355), protein phosphorylation (GO:0006468), the oxidation-reduction process (GO:0055114), transmembrane transport (GO:0055085) and signal transduction (GO:0007165). Highly represented cellular components include the plasma membrane (GO:0005886), cytoplasm (GO:0005737), integral components of the membrane (GO:0016021), cytosol (GO:0005829) and chloroplast (GO:0009507). In the molecular function category, binding activities are highly represented, including protein binding (GO:0005515), ATP binding (GO:0005524), DNA binding (GO:0003677) and metal ion binding (GO:0046872).   KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was performed as previously described [35] to identify the biochemical pathways with which durum target genes were associated (Table S9). The top-ranked enriched pathways are shown in Figure 6. In KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was performed as previously described [35] to identify the biochemical pathways with which durum target genes were associated (Table S9). The top-ranked enriched pathways are shown in Figure 6. In particular, a total of 2406 targets were matched to protein processing in endoplasmic reticulum (ko04141), 2248 targets were matched to starch and sucrose metabolism (ko00500), 2171 targets were matched to the MAPK signalling pathway (ko04010), 1191 target were matched to the spliceosome (ko03040), and 1720 targets were matched to ABC transporters (Table S9). Figure 5. The most represented Gene Ontology terms among the predicted targets in three categories. GO terms were shown for each category based on the number of genes matched for each GO term (from higher to lower). Percent of genes (%) indicates the percentage of the gene number matched for each GO among the total number of genes within each category.
KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was performed as previously described [35] to identify the biochemical pathways with which durum target genes were associated (Table S9). The top-ranked enriched pathways are shown in Figure 6. In particular, a total of 2406 targets were matched to protein processing in endoplasmic reticulum (ko04141), 2248 targets were matched to starch and sucrose metabolism (ko00500), 2171 targets were matched to the MAPK signalling pathway (ko04010), 1191 target were matched to the spliceosome (ko03040), and 1720 targets were matched to ABC transporters (Table S9).

Discussion
Cereal crops are often challenged by multiple biotic and abiotic factors during their life cycle, affecting both production volume and quality [8,[36][37][38][39]. The detrimental effects on the crop plants are often exacerbated when several stress factors are combined, simultaneously or sequentially [8,40]. Increasing evidence has shown that the effects of environmental stress on the parental plants can be carried over to the subsequent generation(s). In the attempt to achieve sustainable crop production, it is essential that we understand the impact of abiotic stress within and across generations, as well as the fundamental networks that underpin crop performance. The influence of the maternal environment on the offspring early growth has rarely been explored in durum wheat before the current study.
For the two genotypes used in this study (DBA Aurora and L6), the parental exposure of water-deficit and heat stress during the reproductive stage significantly reduced the germination potential and germination index. These traits reflect the time-course and the cumulative status of seed germination. Similarly, in soybean, water shortage during the reproductive stage had a pronounced transgenerational effect on seed germination performance [11]. The cumulative percentage of germination negatively correlated with the severity of the drought stress that was imposed on the parents. Germination rate was also significantly reduced in the F1 generation, with reduced seed number and individual seed weight in both soybean cultivars. However, in our study, the germination rate was not affected by the parental stress treatment factor. Moreover, the genotype had a significant impact on all germination traits measured (except for the coleoptile length), with superior performance observed for DBA Aurora. In our previous study, the effect of stress on the parental plants, including physiological responses, yield component and particularly grain traits, was well documented [8]. DBA Aurora produced bigger grains that were full-looking under water-deficit stress and heat stress. The 1000-grain weight of DBA Aurora water-deficit plus heat stress group was significantly higher than that of the control group. On the other hand, L6 had significantly smaller seeds with reduced weight under the effect of such stress (that is, significantly lower 1000-grain weight compared with its control group). For cereal species, high temperature and water shortage during reproduction often leads to an accelerated grain growth rate and a shortened grain development period [40]. As a result, smaller grains with poor grain filling are often produced. However, elite stress-tolerant varieties like DBA Aurora could mitigate the adverse impact of stress and still manage to produce bigger grains, which could (at least in part) contribute to its superior germination performance over L6. It is also worth noting that there is no significant difference between the 1000-grain weight of L6WH seeds and AuWH seeds, suggesting that other factors could contribute to the transgenerational effect of stress rather than grain size.
Similar to germination traits, early durum seedling growth of the offspring was also significantly affected by the parental stress exposure. Genotype seems to play a significant role here. Although seedling vigour was significantly reduced under the transgenerational effect of stress in both genotypes, a greater impact was always observed in L6. In particular, a significant interaction was found between genotype and parent treatment for all root traits measured in the young seedlings. Only L6 has shown reduced root weight due to parental stress exposure. Such genotypic differences in the transgenerational effects of stress have also been noted before [41]. The findings suggest that breeders and growers should take care in selecting varieties under different parental environments for crossing or seed production (this is often overlooked, with a 'mother stock' source used for such purposes). Furthermore, the stress exposure of the parents increased the young seedling root length for both genotypes (with a higher increase in DBA Aurora), suggesting that root elongation could be promoted through transgenerational stress exposure irrespective of genotypes. Establishing a longer, deeper root system at earlier developmental stages might be an adaptive mechanism to prepare the offspring to improve water and nutrient uptake in any possible future stress occurrence [11]. Previous research has suggested that the impact of parental stress on early seedling growth is often associated with seed storage reserves such as carbohydrates and proteins. However, in our case, both the DBA Aurora and the L6 parents produced grains with increased grain protein content and reduced starch content under the same stress without any significant genotypic difference [8]. In the miRNA-target analysis, starch and sucrose metabolism came second for the most-enriched KEGG pathways (Table S9). The 1639 targets enriched for this pathway included many enzymes, such as sucrose synthase, starch synthase, beta-glucosidase, sucrose phosphatase, and isoamylase (Table S8). Altered enzyme activities regulated by differentially expressed miRNAs can affect how storage reserves were utilised during germination and early growth. The results suggest that, apart from seed quality traits being affected by stressful paternal environments, epigenetic mechanisms such as miRNA-target modules could be in play during the transmission across generations [11].
As a central class of epigenetic regulators, miRNAs take part in a wide range of cellular biological processes [17,18,42]. Hence, the identification and annotation of the miRNAome enables the epigenetic networks underlying the transgenerational effect of stress to be illustrated. The combination of NGS and bioinformatics analysis using genome reference has been the most effective approach to characterise the miRNA population on a genome-wide scale [17]. Such technologies have led to the complete discovery of small RNAs in many plant species on an unprecedented scale. However, despite a few reports in durum wheat (including our previous studies) [20,22,30,34], the discovery of novel miRNAs and a fully comprehensive miRNAome description in this important cereal has not been previously achieved due to incomplete genome information. Only recently, the ground-breaking genome assembly of the modern durum cultivar Svevo was made available [33]. The complete Triticum turgidum genome presented us with a timely opportunity to harness additional information from sRNA sequencing datasets. For the first time, to our knowledge, a systematic and complete archive of durum miRNAs has now been provided. In total, we successfully identified 1717 conserved and 485 novel miRNA entries matched with their MIR loci origins. This number exceeded all previous reports in durum wheat. Interestingly, among all the conserved miRNA entries, the top plant species with the highest miRBase conservation match was Glycine max, followed by Oryza sativa. Hexaploid wheat (Triticum aestivum), with tetraploid wheat as one of its genome progenitors, only came third in conservation ranking. This is likely due to soybean miRNAs being very well-documented (684 precursors, 756 mature) in the miRBase when compared to bread wheat (122 precursors, 125 mature), which is not surprising considering that the Glycine max genome was well-annotated ahead of all wheat species sequenced. Notably, there is currently only one Triticum turgidum miRNA entry (ttu-miR160, MI0016454) deposited in the miRBase [43]. It could be foreseen that this number will increase rapidly with the first durum genome assembly now being available to the scientific community.
Plant miRNAs can act rapidly in response to both environmental and developmental signals [17]. In particular, various studies have described the functional roles of miRNA-target modules in combatting water deficiency and extreme temperatures in cereal crops [8,23,24]. In our previous research, we have demonstrated that the expression of durum miRNAs can change differently with different timings of stress adaption within the life cycle [23,24]. Such patterns were also genotype-dependent among stress-tolerant and -sensitive varieties, suggesting that miRNA-mediated regulatory modules could contribute to the genotypic differences in the physiological performance, ultimately affecting yield production [23]. The present study further illustrated that durum miRNAs behave in a genotypic manner under the effect of transgenerational stress. Overall, there were more up-regulated miRNAs (275, p < 0.05) than down-regulated miRNAs (181, p < 0.05) in the stress-tolerant genotype, but more down-regulated miRNAs (242, p < 0.05) than up-regulated (153, p < 0.05) miRNAs in the stress-sensitive genotype. Durum miRNAs that were induced by the transgenerational effects could repress the expression of protein-coding genes with negative regulatory roles, while reduced miRNA expression could promote the expression of target genes to increase protein products with positive functions. A higher number of induced miRNAs in the stress-tolerant genotype and higher number of repressed miRNAs in the stress-sensitive genotype suggest that contrasting mechanisms were preferably deployed in response to transgenerational stress, resulting in the difference in young seedling growth.
Investigations into the expression pattern of specific durum miRNAs, together with the molecular functions of their predicted targets, provide further information on the possible pathways activated in different genotypes. As an example, tae-miR1847-5p was significantly down-regulated in DBA Aurora (−0.90 fold) under the transgenerational effects of stress (Table S6). In silico analysis revealed that tae-miR1847 targets several CREB-binding protein-like family members located on the 6A and 6B chromosomes (Table S7). CREB-binding proteins (CBPs) are a type of transcriptional coactivator that can switch on the cAMP-dependent signalling pathways, which are key components employed by many hormone-regulated activities to modulate gene expression [44,45]. CBPs initiate the function of many classes of transcription factors through their intrinsic transferase activity to acetylate histones, which is a key process associated with chromatin remodelling [44]. Stress-related chromatin modification plays a critical role in transgenerational memory in plants [15]. The repressed expression of tae-miR1847 in the stress-tolerant genotype would increase the expression of its target CBPs, which could ultimately promote histone acetyltransferase-regulated transduction pathways. However, this pathway might not be in play in the stress-sensitive genotype, as tae-miR1847-tp was up-regulated (0.23-fold) under the transgenerational effects of stress. Another example of a miRNA-regulated module that could contribute to the differential seedling performance involved mitogen-activated protein kinase (MAPK). MAPK signalling is one of the top KEGG pathways enriched among the predicted targets. MAPK cascades are a universal signal transduction mechanism that play diverse roles in connecting intraand extra-cellular crosstalk in eukaryotes [46]. MAPKs are involved in virtually all the important aspects of plant development and response to environmental signals. In particular, MAPKs facilitate plant adaptation and survival under many abiotic stresses, including drought, heat and salt stress, via coordinating the transduction of secondary messengers and hormone signals [46,47]. Here, a few durum miRNAs that target MAPKs exhibited genotype-dependent expression patterns. For example, conserved miRNAs such as ttu-miR160, hvu-miR444b, and tae-miR1847-5p were all down-regulated (−0.42-fold, −0.45-fold and −0.90-fold, respectively) in DBA Aurora under the effect of transgenerational stress (Table S6). This promotes the expression of MAPKs in the stress-tolerant variety to facilitate signal transduction. However, the expression of these miRNAs was either unchanged (ttu-miR160 and hvu-miR444b) or increased (tae-miR1847-5p) in L6, suggesting that miRNA-MAPK modules were not activated in the stress-sensitive genotype. Research in bread wheat has also demonstrated the involvement of MAPK in transgenerational tolerance to heat stress [48]. Increased expression of MAPK in the progenies of heat-stressed parents contributed to the enhanced activation of Ca 2+ signalling pathways to facilitate stress adaptation [48]. Moreover, stress-induced MAPK could promote root tip growth through controlling the organisation of the actin cytoskeleton [49]. Durum miRNA-promoted MAPK expression could possibly contribute to the adaptive root traits observed in the stress-tolerant genotype. Future research should investigate the co-expression of durum miRNAs and MAPK family members in different tissues types to further our understanding of miRNA-MAPK modules in transgenerational stress adaption.
The spliceosome (ko03040) is another major KEGG pathway enriched for the targets of differentially expressed durum miRNAs. The pre-mRNA splicing machinery is governed by the spliceosome, a large ribonucleoprotein complex that recognises the pre-mRNA splice sites to remove introns and join exons [50]. Research has shown that alternative splicing (AS), which enables the production of multiple mRNA isoforms from one single gene, is closely associated with abiotic stress adaption [51,52]. The AS process provides a fine-tuning mechanism to regulate the levels of active and inactive isoforms of stress-associated proteins. The specificity of splice-site selection in AS largely relies on serine/arginine-rich (SR) proteins (splicing enhancer) and heterogeneous nuclear ribonucleoproteins (hnRNPs, splicing suppressor) [53]. In our study, a large number of SR and hnRNP family members were targeted by differentially expressed durum miRNAs. Some miRNAs exhibited a genotype-dependent expression, while others seem to have universal regulatory roles. For example, ttu-miR160 targets several hnRNP-encoding genes located on the 4A and 4B chromosomes (Table S7). It also targets thirteen SR splicing factor-coding genes located on chromosome 6A and 6B. ttu-miR160 was down-regulated in DBA Aurora (−0.42-fold) in response to transgenerational stress, indicating that the expression of SR proteins and hnRNPs were possibly increased to coordinate the AS process. ttu-miR160-regulated AS was not in play in the sensitive genotype, based on the unaltered miRNA expression. However, a few other miRNA-regulated modules were active in both genotypes. For example, tae-miR9772, ata-miR396c-5p, and a miR396c variant (sbi-miR396c_R+1) all target hnRNP 1-like genes. The expression of these miRNAs was down-regulated in both genotypes, suggesting that hnRNP-modulated splicing events were activated by miRNAs without genotypic specificity. Furthermore, stress-induced AS can regulate the processing of non-coding RNA transcripts, so as to control the production of different pre-miRNA hairpins [54]. In barley, heat stress induced the splicing of introns hosting the miRNA primary transcripts (pri-miRNAs) of miR160a and miR5175a, which was correlated with the accumulation of mature miRNA products [55]. The results suggest that the durum miRNA-modulated AS mechanism could possibly provide a feedback loop to fine-tune the post-transcriptional regulation of miRNA precursor processing. Looking forward, the complexities in durum miRNA expression, spliceosomal protein activities, and the production of pri-miRNA isoforms require further exploration. In addition, it would be desirable for future work to study the upstream regulatory mechanisms of durum miRNAs so as to understand the underlying reasons for miRNA expression changes between stress-tolerant and stress-sensitive genotypes. Genomic analysis of the surrounding regions (such as transcription start sites and cis-regulatory motifs) of major stress-responsive miRNAs can compare the promoter features of the stress-tolerant and stress-sensitive durum genotypes to provide valuable information on transcriptional regulation of MIR genes. Our current research has paved the way for future investigation, having provided in-depth details of the MIR gene (genome ID, length, pre-miRNA sequence, secondary structure etc., Table S1).
In conclusion, the present study provides a valuable resource for an improved understanding of the transgenerational effects of water-deficit and heat stress on germination and seedling growth in durum wheat, as well as new insights of conserved and novel miRNAs that are involved in the epigenetic regulatory networks. Such knowledge provides a new foundation for durum researchers and breeders to use effective strategies for developing and selecting elite varieties that can mitigate or adapt to the negative impacts of stress that are transmitted across generations.

Plant Materials
Seeds of two Australian genotypes were used in this study. The seeds of both genotypes were collected from the parental plants in the control group and the water-deficit and heat stress group from the previous year's experiment [8]. Briefly, water-deficit stress was applied to the parents from the booting stage by maintaining the soil water content at half of the field capacity. Heat stress was applied post-anthesis by exposing the plants to 37/27 • C for 24 h at 5, 15, 25, 35, and 45 days after flowering (a total of five stress episodes on the same plants). In total, four different seed groups were used in this study-AuCG (seeds from DBA Aurora control group), AuWH (seeds from DBA Aurora stress group), L6CG (seeds from L6 control group), and L6WH (seeds from L6 stress group).

Germination Test
Germination tests of the four seed groups were performed according to the International Rules for Seed Testing 2019 (https://www.seedtest.org). The seeds used for germination were randomly selected. The seeds were sterilised by rinsing with 70% ethanol for 5 min, followed by a wash with sterile distilled water. The seeds were germinated on filter paper in petri dishes containing sterile distilled water for 96 h at 22 • C. Three replicates per group were used with 30 seeds per replicate. germination potential (Gp, %) was determined as number of seeds germinated within 24 h/total number of seeds × 100%. The germination rate (Gr, %) was determined as number of seeds germinated within 96 h/total number of seeds × 100%. The germination index (Gi) was determined as (n t /d t ), where n t represents the number of germinated seeds on the t day, and d t represented germination days [56]. The length of coleoptile and young roots was recorded at 96 h.

Measurement of Seedling Establishment
For each group, 16 randomly selected germinated seeds were transferred into pots (one seedling per pot). Each pot contained 1.2 kg of N40 sand with 0.5% CaCO 3 , and nutrient solution was provided as previously described [8]. The seedlings were grown under controlled glasshouse conditions with day/night temperature of 24/18 • C and a photoperiod of 16/8 h. Seedling height, shoot fresh weight and root fresh weight were measured on 12-day old seedlings. Shoot dry weight and root dry weight were measured after overnight oven drying at 65 • C. Fresh shoots were collected and snap-frozen in liquid nitrogen for RNA extraction. Eight biological replicates were used for seedling measurements and eight biological replicates were used for RNA extraction. Seedling vigour index I (SVI-I) was determined as germination rate × seedling fresh weight in mg [57]. Seedling vigour index II (SVI-II) was determined as germination rate × seedling dry weight in mg [57]. Two-way-ANOVA was used to analyse the statistical significance (p < 0.05) for germination and seedling data using Genstat (19th Edition, VSN International).

Target Prediction and Functional GO Analysis
To predict protein-coding genes targeted by durum miRNAs, in silico analysis using the GSTAr was employed. Gene Ontology terms of the miRNAs targets were annotated according to their biological process, cellular component and molecular function. KEGG Pathways (http://www.genome.jp/kegg) of miRNA targets were also analysed. The statistical enrichment of GO terms and KEGG pathways was performed as previously described [35,60].
Supplementary Materials: The following are available online at http://www.mdpi.com/2223-7747/9/2/189/s1. Figure S1: Heatmap of differentially expressed miRNAs (DEMs) identified in comparison across four libraries (TCG vs TWH vs SCG vs SWH). The colour key (ranging from blue to red) represents the miRNA abundance using the log10 (normalised count) value, Table S1: Summary of 2202 MIR-miRNA entries identified in this study, Table  S2: Conserved MIR families identified in this study, Table S3: List of plant species that were used as references in the miRBase, Table S4: List of miRNAs that were differentially expressed across TCG vs TWH vs SCG vs SWH, Table S5: List of miRNAs that were specifically expressed in DBA Aurora, L6, CG and WH seed sources, Table S6: List of differentially expressed miRNAs in comparisons made between two libraries (TWH vs TCG, SWH vs SCG, SCG vs TCG, SWH vs TWH), Table S7: Target genes predicted for differentially expressed miRNAs, Table S8: GO enrichment of predicted target genes in three categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), Table S9: KEGG pathway enrichment analysis of predicted target genes.