The ability of organisms to produce distinct phenotypes under different environmental conditions is referred to as phenotypic plasticity (Schlichting, 1986; Sultan, 2000; Pigliucci, 2005), which is crucial for organisms to survive in a world that is constantly changing, especially when the environment changes quickly (Ghalambor et al., 2007). It is usually considered that plants have greater plasticity than animals for their sessile lifestyle (Bradshaw, 1965; Sultan, 2000; Borges, 2008). Plants are forced to cope with local environmental fluctuations by plastic modifications in morphology, physiology, behavior or life history traits, to overcome the challenges from changing environments and increase their niche breadth (Richards et al., 2006; Ghalambor et al., 2007). Over the past two decades, the molecular and cellular mechanisms of environmental perception and gene regulation underlying various plastic responses in plants have been an explicit focus of the ecological developmental (Eco-Devo) research, with the primary goal of characterizing the developmental genetic pathways supporting the ability to respond plastically to environment changes (Sultan, 2007, 2010; Gilbert, 2016). Understanding the regulatory events informed by both external and internal environmental signals and identifying genetic elements implicated in environmentally contingent development will provide key insights to the immediate tolerance and potential evolutionary resilience of organisms to environmental fluctuations (Sultan, 2007, 2010; Gilbert, 2016).

Investigating the development of plastic traits in ecological context has revealed an association between environmentally contingent phenotypic expression and the environmental sensitivity of gene expression (Morris and Rogers, 2014). This environmental sensitivity of gene expression is, in turn, driven by the flexibility of epigenetic gene regulation in fluctuating environments (Bossdorf et al., 2008; Richards et al., 2017). Epigenetic variations triggered by environmental perturbations can lead to inducible phenotypic changes by altering transcriptional profiles and eventually altering ontogenetic trajectory in response to environmental variation, even in the complete absence of genetic variability (Kalisz and Kramer, 2008; Aubin-Horth and Renn, 2009), and thus serve as a mechanistic link between genes and environment (Schlichting and Smith, 2002; Morris and Rogers, 2014). A suite of epigenetic mechanisms, including DNA methylation, histone modifications and small RNAs, have been shown to be able to control the temporal, spatial and abundance patterns of gene expression and RNA-translation under varied environments (Schlichting and Smith, 2002; Chinnusamy and Zhu, 2009). These epigenetic mechanisms may act separately or concomitantly (Bilichak et al., 2012; Richards et al., 2017).

microRNAs, a type of endogenous small-non-coding RNAs, have emerged as key post-transcriptional regulators of gene expression through base pairing with their complementary mRNA targets (Jones-Rhoades et al., 2006). Functional genomic studies have shown the involvement of plant microRNAs in a broad range of developmental processes (Eldem et al., 2013; Jin et al., 2013; Li and Zhang, 2016). miR165/166 have been proved to participate in SAM development, acting as mobile signal molecules for stem cell maintenance (Zhang and Zhang, 2012). miR319 mediates the change of plant leaf shape via targeting TCP transcription factors (Palatnik et al., 2003), while miR160 and miR393 play key roles in root architecture regulation via post-transcriptional modification of the auxin signal pathway genes (Wang et al., 2005; Bian et al., 2012). Both miR156 and miR172 are involved in regulation of flowering time by age-dependent pathway (Wang, 2014; Li and Zhang, 2016). Multiple microRNAs, including miR169, miR164, miR319, miR159, and miR167 function to control floral organ identity fate during flower development by regulating the spatial boundaries of expression of target genes, and to specify particular cell types during later stages of flower development (Nag and Jack, 2010). It has also been revealed that a single miRNA may have diverse regulatory roles. For example, miR159 is involved in flowering time control, leaf development, seed size and shape determination, respectively (Blazquez et al., 1998; Palatnik et al., 2003; Achard et al., 2004).

In addition to control various developmental processes, plant microRNAs also act as node molecules to coordinate the balance between plant development and environmental clues, playing critical roles in response to biotic and abiotic stresses (Sunkar et al., 2012; Zhang, 2015; Li and Zhang, 2016; Noman et al., 2017), microRNAs involved in responses to water-logging/flooding have been demonstrated in different plants. For instance, Zhang et al. (2008) and Zhai et al. (2013) identified microRNAs responding to waterlogging in maize roots. Khan et al. (2014) predicted the associated microRNAs under waterlogged conditions in sugarcane. Ren et al. (2012) found significant changes in the expression of seven conserved microRNA families and five novel microRNAs in response to flooding stress in Populus tomentosa. Jin et al. (2017) revealed the microRNA-mediated regulatory networks in lotus by identifying submergence-responsive microRNAs and their targets. Flooding exposes plants to severe hypoxia. Thus, microRNAs mediating hypoxic stress responses in plants were also thought to facilitate the survival and adaptation to the flooding condition by modifying gene expression, metabolic networks, and developmental processes (Moldovan et al., 2009; Licausi et al., 2011).

Alternanthera philoxeroides is an exotic amphibious weed that is native to South America but has now invaded into the temperate and tropical regions across the world (Julien et al., 1995; Pan et al., 2007). A. philoxeroides rarely produces viable seeds in its introduced range, and propagates mainly via vegetative regeneration (Julien and Stanley, 1999; Geng et al., 2007; Pan et al., 2007). Plants produced by vegetative propagation are genetically identical. However, these plants can grow in diverse habitats, from purely aquatic environments (ponds, rivers) to dry lands (Pan et al., 2007). Individuals growing in aquatic and terrestrial habitats exhibited notable morphological differences, with the plants growing in aquatic environments showing significantly longer internode length and greater stem diameter (Geng et al., 2007; Gao et al., 2010). Plasticity in internode elongation, aerenchyma formation, and other morphological and physiological traits allows the amphibious plants of A. philoxeroides to withstand alternative wetting and drying conditions and to proliferate in both aquatic and terrestrial environments (Geng et al., 2006, 2007; Pan et al., 2007).

As the most conspicuous responses of plants subject to submergence, enhanced internode/petiole elongation and secondary aerenchyma formation have been found in many species (Jackson and Armstrong, 1999; Van Veen et al., 2013; Yamauchi et al., 2013; Voesenek and Bailey-Serres, 2015). The molecular and cellular mechanisms underlying these responses, such as the external and internal environmental cues, signaling pathways, and genetic elements implicated in developmental outcomes, have been explored in deepwater rice (Nishiuchi et al., 2012), Rumex (Peeters et al., 2002), and other species (Bailey-Serres and Voesenek, 2010; Voesenek et al., 2014; Loreti et al., 2016). However, the epigenetic mechanisms underpinning submergence-induced plastic development, including the potential roles of microRNA-mediated regulation of gene expression, remain elusive. Since the transition from aquatic to terrestrial environments, or vice versa, requires timely and flexible adjustments in morphological traits, physiological processes, and developmental machinery to extremely different conditions (Boeger and Poulson, 2003; Nielsen and Nielsen, 2006), the epigenetic control system is most likely to provide an effective short-term strategy for amphibious plants to plastically respond to environmental fluctuations by reversibly controlling effector gene expression states depending on environmental conditions.

The genetically identical but phenotypically distinct plants of A. philoxeroides provide a suitable model for us to rule out genetic contributions to observed differences, to explore the correlations between environmental stimuli, epigenetic modification and phenotypic variation. Previously, we have investigated the gene expression reaction norms of A. philoxeroides in contrasting hydrological conditions, and identified 11 transcriptionally coordinated gene groups associated with environment-induced phenotypic variation in A. philoxeroides (Gao et al., 2015). We have also examined genome-wide DNA methylation patterns of A. philoxeroides in natural and manipulated habitats, and found that plants of different source populations not only underwent significant morphological changes in common garden environments, but also underwent genome-wide DNA methylation alterations in response to different water treatments, demonstrating a correlation between epigenetic reprogramming and the reversible phenotypic response of A. philoxeroides to particular environmental factors (Gao et al., 2010; Deng et al., 2014). In this study, we aim to detect microRNAs potentially involved in the regulation of the plastic responses of A. philoxeroides to contrasting hydrological conditions, especially microRNAs associated with enhanced internode elongation growth under submergence. By preliminarily profiling the time-dependent changes of A. philoxeroides miRNome in alternate wetting and drying conditions using deep sequencing technology, and validating the most differentially expressed microRNAs by independent quantitative reverse transcription-PCR (qRT-PCR), we identified the submergence-sensitive microRNAs of A. philoxeroides, and determined microRNA:mRNA functional pairs and their temporal expression dynamics possibly associated with submergence-induced internode elongation in A. philoxeroides. This expanded investigation enhanced our understanding of the role of microRNAs in mediating adaptive developmental response to environmental changes, promoting the study of epigenetic regulation and phenotypic variation in an ecological context.

Due to their sessile lifestyle, plants have developed a remarkable plasticity for shaping their body plan to adapt to changes in their immediate environment. The high sensitivity of microRNAs to environmental perturbations and their significant role in readjusting the expression levels of genes and phenotypic profiles under changing environments have made them the best candidates as dynamic regulators of phenotypic plasticity (Rubio-Somoza and Weigel, 2011).

Genome-wide perspective of miRNAome in response to pond and upland treatments revealed 146 submergence-responsive microRNAs belonging to 49 microRNA families in A. philoxeroides. The members from 19 families (such as MIR159, MIR167_1, MIR169, MIR396, MIR482, MIR530, and MIR3630) were not only involved in responses to water-logging/flooding in A. philoxeroides, but also in other plant species. For instance, gene expression profiling and qRT-PCR analyses revealed the varied expression patterns of miR169u and miR396a-3p in A. philoxeroides during submergence. The similar patterns were also found in Nelumbo nucifera (Jin et al., 2017) and Zea mays (Zhang et al., 2008; Liu et al., 2012). Three members of the miR167_1 family (miR167b-3p, miR167d-5p, and miR167d-3p) showed significantly altered expression between the pond- and upland-treated samples in A. philoxeroides. Similarly, the members from this family were differentially expressed under flooding condition in P. tomentosa (Ren et al., 2012). These results suggest the common roles of these microRNAs in different plant species. The submergence-responsive microRNAs identified in A. philoxeroides can also be induced by components of submergence, such as hypoxic stress (e.g., miR159, miR169, miR319, miR394, and miR395) (Licausi et al., 2011), elevated CO₂ concentration (e.g., miR156, miR157, miR172, and miR160) (May et al., 2013), darkness (e.g., miR160, miR167, miR171, miR396, and miR398) (Huo et al., 2015) and ethylene entrapment (e.g., miR399, miR408, miR482, miR4414, and miR5139) (Pei et al., 2013).

Target prediction and functional annotation of the submergence-responsive microRNAs recognized microRNAs potentially associated with plastic internode elongation in A. philoxeroides upon submergence. The target of novel_mir_26 was homeobox protein 22 (HB22), which has been shown to be associated with the synthesis of gibberellic acid (GA) in Arabidopsis (Bueso et al., 2014). The primary action of GA in stem growth was on cell elongation (Arney and Mancinelli, 1966). The expression of A. philoxeroides HB22 (Contig539) was significantly up-regulated at 12 h of the pond treatment accompanied by down-regulation of novel_mir_26, that might promote internode elongation in submerged plants by enhancing GA biosynthesis (Voesenek and Bailey-Serres, 2015). Both cytokinin and brassinosteroid (BR) were also implicated in stem elongation growth (Azpiroz et al., 1998; Hansen et al., 2009). The target of novel_mir_19 in A. philoxeroides was contig16589, which was homologous to Arabidopsis AHP4 and potentially involved in cytokinin signaling pathway (Hutchison et al., 2006). Novel_mir_19 and its target exhibited a significant inverse expression pattern under the pond treatment (Figure 4E), and the target cleavage site was verified by 5′ RLM-RACE. A. philoxeroides contig2589, a homolog of Arabidopsis CDL1 that positively regulates BR signaling (Kim et al., 2011) and stem elongation (Peeters et al., 2002), was predicted as the target of miR398a-5p. The expression levels of contig2589 and miR398a-5p were also changed accordingly in an inversely correlated pattern according to the expression profile.

Abscisic acid (ABA) antagonizes the effects of GA and other phytohormones in numerous processes during plant growth and development (Bailey-Serres and Voesenek, 2010; Finkelstein, 2013). Growth-regulating factors (GRFs), a conserved class of plant-specific transcription factors, participate in ABA signaling pathway as a transcriptional repressor of ABA- and osmotic-stress responsive genes (Kim et al., 2012). The expression of several GRFs is controlled by miR396 (Omidbakhshfard et al., 2015), functioning as a regulator of hormone homeostasis and signaling under flooding stress (Zhang et al., 2008). In our study, miR396-3p presented fluctuate expression patterns opposite to contig34390 (GRF6) during the pond treatment (Figure 4C), indicating the interplay between ABA and other phytohormones in A. philoxeroides upon submergence. Previous studies have shown that interactions between different plant hormones play a central role in generating plastic phenotypic responses to environmental variation (Sultan, 2000; Lema, 2014), including in promoting shoot elongation upon submergence (Van Der Straeten et al., 2001; Voesenek et al., 2003; Cox et al., 2004, 2006; Jackson, 2008; Van Veen et al., 2013). The detection of dynamic expressions of microRNAs and their corresponding targets associated with different hormone signaling pathways in this study possibly suggests the potential role of A. philoxeroides microRNAs in coordinating multiple phytohormone responses during the submergence treatment. The crosstalk between microRNAs and hormone signaling cascades deserves better attention while studying submergence-induced elongation growth in flooding-tolerant plants.

Previous gene expression profiling showed that a large number of genes associated with different processes of cell wall strengthening, such as the xylan biosynthetic process, the lignin biosynthetic process, secondary cell wall biogenesis, etc., exhibited a significant “up-down” pattern in the early stage of the pond treatment, suggesting a “Just-in-time” cell wall thickening in response to submergence though its function remaining unclear (Gao et al., 2015). In coincidence with the transient cell wall strengthening, we found in this study that NF-YA3 (A. philoxeroides contig18972), a transcript factor which has been shown to impair cell elongation (Leyva-González et al., 2012), was transiently up-regulated in the early stage of the pond treatment (Figure 4B), inversely correlated with the expression of the corresponding miR169u. In contrast, the expression of SKS5 (A. philoxeroides contig45839), a gene that has been identified as a cell wall-related gene (Mutwil et al., 2009) and was expressed most strongly in expanding tissues (Sedbrook et al., 2002), was highly repressed in the early stage of the pond treatment accompanied by induced and transient up-regulation of the corresponding miR167d-5p (Figure 4A). The transient down-regulation of the cell wall expansion-promoting gene SKS5 was consistent with the transient up-regulation of cell wall strengthening genes. microRNAs, while functioning as a trigger of hormone response to external environment, may also coordinate the expression of genes associated with different processes of cell wall modification under submergence.

However, the relationships between microRNAs and their target genes was not always straightforward. One microRNA seemed to be able to regulate more than one target gene and, in turn, one gene might be regulated by more than one microRNA. For example, both contig49354 (DUF579) and contig18934 (UDP-Glycosyltransferase superfamily protein) were identified as the targets of novel_mir_25 and showed similar inverse relationships with novel_mir_25 in expression (Figure 4F), while the expression of contig2773 (SPL3) was regulated by miR157a and miR4387b, respectively, at different time points of the treatment. We also found a positive relation between the expression of novel_mir_29 and its target Contig23974 (vacuolar invertase, VI2), a gene that has long been considered as a major player in cell expansion (Wang et al., 2010) (Figure 4H). The most striking feature revealed by time course expression profiling was the fluctuating pattern in microRNA expression. Most of the differentially expressed microRNAs were transiently up-/down-regulated at certain time points during the treatments. More than half of the submergence-responsive microRNAs (e.g., miR1023a-3p, miR3630-3p, and novel_mir_25) responded to the pond treatment at the early stage of the treatment (1–12 h). About one-third of the microRNAs were differentially expressed at the 1st hour from the start of treatments, while 19 microRNAs (e.g., miR1439, miR7542, and novel_mir_48) responded to the treatment 120 h later. The fine-scale temporal changes in microRNA expression highlight the importance of time-series sampling in identifying stress-responsive microRNAs and analyzing their role in stress response/tolerance. Incomplete sampling might give a biased conclusion about the microRNA’s role in readjusting the expression levels of genes and phenotypic profiles under changing environments.

GL performed the wet lab work and data analysis, and drafted the manuscript. YD, YG, and CZ helped to set up experiments. YW, WZ, and ZS participated in data analysis. LG and JY conceived the idea, participated in the design of the study and finalized the manuscript. All of the authors read and approved the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer OV-L and handling Editor declared their shared affiliation.