White spruce (P. glauca), a keystone species of boreal forests in the North American taiga, is anemophilous (outcrossing) and its seed and pollen cones develop on the same tree and are diclinous (i.e., unisexual). According to seed orchards station records, four populations of P. glauca (Pop 1~4), characterized by different pollination timing and seed developmental duration, were chosen and 20 developing cones for each population were collected at early, middle, and late developmental stages for a total of four timepoints at the Kalamalka Research Station seed orchards (50°–51°37′N, 119°16′–120°29′W), British Columbia, Canada (Figure 1). Note that population 1 only had three sampling timepoints due to its short seed developmental span. Accordingly, the climatic data were retrieved from the Station. Likewise, as per contrasting seed-set durations, we selected two wild strains of a model annual selfing plant A. thaliana, Cvi-0 and Col-0, originating from the Cape Verde Islands and Columbia (Missouri, USA), respectively. The two Arabidopsis ecotypes exhibit contrasting dormancy intensities (Koornneef et al., 2000; Ali-Rachedi et al., 2004). They were cultivated in 2-in planting trays containing soil mixed with slow-release fertilizer 14-14-14 in a growth chamber with 16/8 h day/night photoperiod, PPFD (photosynthetic photon flux density) of 250 μmol·m²·s⁻¹, and constant temperature of 22°C. Individual flowers were tagged on the day of flowering and developing seeds were sampled manually every day for 9 consecutive days after pollination (DAP) (Figure 1). Five more sampling time points (10–14 DAP) were performed for Cvi-0 due to its slow seed development (Figure 1). Developing seeds at each timepoint were sampled three times (i.e., three biological replicates). After extraction from white spruce cones, direct collection of siliques at Arabidopsis embryogenesis (0–6 DAP) and quick hand-dissection in water from ~30 inflorescences at maturation stages (7 DAP onward), entire developing seeds or siliques were immediately frozen in liquid nitrogen and stored at −80°C until further use. Note that the seed samples from the same developing timepoint were pooled in the same vial for the subsequent analyses.

Total RNAs were extracted and divided into two aliquots (~15 μg each) from developing seed samples of P. glauca and Arabidopsis using PureLink Plant RNA Reagent (Ambion) according to the manufacturer's instructions. The intergrity and quantity of the RNA samples were assessed on a BioAnalyser 2100 (Agilent Technologies) and a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific). The sRNA-seq libraries were constructed using a strand-specific and plate-based protocol. To enrich sRNAs, total RNA samples underwent polyA selection using Miltenyi MultiMACS mRNA isolation kit (cat. 130-092-519) following the manufacturer's protocol and the flowthrough (i.e., containing sRNA species without mRNAs) was used for plate-based sRNA construction. A 3′ adapter that is an adenylated single-stranded DNA was selectively ligated to the sRNA template using a truncated T4 RNA ligase 2 (NEB Canada, cat. M0242L). A 5′ adapter was then added using a T4 RNA ligase (Ambion USA, cat. AM2141) and ATPs. After ligation, first strand cDNA was synthesized using a Superscript II Reverse Transcriptase (Invitrogen, cat. 18064 014) and one RT primer. This was the template for the final library PCR, into which 6-nt mers index was introduced to identify libraries (i.e., demultiplexed) from a sequenced pool. Constructed libaries were pooled by phylum; that is, 15 P. glauca and 25 Arabidopsis samples were pooled seperately, and both lanes were 31 base SET lanes. Sequencing (Illumina HiSeq™ 2500) was implemented using one short SET indexed lane per pool (BC Cancer Agency, Genome Sciences Centre, Vancouver, Canada).

The sequence data were partitioned into individual libraries based on the index read sequences, and the reads underwent an initial QC assessment. After being preprocessed to clean reads by trimming adapters and barcode sequences using an internal matching algorithm (BC Cancer Agency), the raw sequencing data (in bam format) were converted into sam, fastq, fasta, and txt formats under Linux in a command-line environmemt for subsequent use. The sRNA toolbox was used to profile sRNAs and size distribution, perform conserved miRNA analysis using miRNAs for Arabidopsis or a high confidence set of miRNAs on miRBase for P. glauca, and their consensus miRNA families and unique sequences (Rueda et al., 2015). sRNAs in sequencing libraries were computationally predicted against the P. glauca genome assemblies (highly framented, containing 3.1 million scaffolds longer than 500 bases and a N50 of 43.5 kbp; PG29 v3, 20Gb divided into 30 Mb per file) (Warren et al., 2015) and the Arabidopsis genome (TAIR10), respectively, using miRPlant (An et al., 2014) at default settings with a sequence length cutoff of 18~24 nt. As the size of our P. glauca sequence libraries exceeded the maximal single load on miRPlant, we divided each library into several sub-ones with a maxium size of 120 Mb under Linux and after combining the output files in the same library (e.g., summing up raw abundance for the same unique reads from different sub-files), duplicated reads were removed and we only retained one copy of unique reads with the highest prediction score for each library using R 3.2.2 (The R Project for Statistical Computing), followed by a visual inspection to ensure the R coding has attained our objectives.

The libraries substracted r/t/sn/snoRNAs (i.e., ribosomal RNA, transference RNA, small nuclear ribonucleic RNA, and small nucleolar RNA, originally sourced from Sanger RNA family database 12.0, Nawrocki et al., 2015) and non-sRNAs from the total library size, and the resulting number was used for abundance normalization. Abundances were expressed throughout this study in reads per millon (RPM) unless otherwise indicated. Due to the unavailability of complete r/t/sn/snoRNA sequences in P. glauca, we utilized Arabidopsis r/t/sn/snoRNA annotations to classify non-sRNAs in P. glauca, which include reads of non-predictable secondary RNA structure and non-mapped to the genome. After genome-wide identification of sRNAs, the unique sRNA sequences with raw abundance above 10 at least in one library, along with raw and nomalized counts and precursors, were archived and compared with previously identified spruce sRNAs (Källman et al., 2013; Xia et al., 2015). To isolate conserved miRNAs by using a homolog search, sRNAs in this study were aligned versus a Viridiplantae-specific miRBase reference file containing 100,014 miRNA sequences from 1,397 miRNA families (Chávez Montes et al., 2014), curated in Table S1. Note that miRBase v21 (http://www.mirbase.org/) has archived c. 8,000 miRNAs from 73 plant species.

Heat maps graphically depicted the most conserved and differentially expressed miRNAs at seed set across populations/ecotypes, whereby cluster analyses were performed for seed-set phases and key miRNAs within phylum using an euclidean method. Note that raw processed data (i.e., log2 of count per million) underwent a log10-transformation before being fitted on the map.

Top-one targets of the archived unique sRNAs were predicted using transcripts without miRNA genes on psRNATarget (Dai and Zhao, 2011). We focused on the sRNAs of high strength of prediction (score ≥ 0) in Arabidopsis, while in P. glauca, on sRNAs that were detected in at least 14 of 15 libraries across four populations, termed “most conserved” sRNAs hereafter (N.B. total reads in one library are much less than the others). To annotate target mRNA functions, the top predicted target gene for each sRNA of interest was aligned against the Gene Ontology (GO) protein database for GO term classification and KEGG pathway enrichment (Ashburner et al., 2000; Kanehisa and Goto, 2000).

In a parallel analysis to search for potential sRNAs targeting genes involved in adaptation to climate, we manually refined sRNAs through identifying their target genes that function in stress responses for Arabidopsis. Reportedly, there were 73 key adaptive genes in P. glauca (17 overlapped between two articles, Hornoy et al., 2015; Yeaman et al., 2016). We adopted the pipeline of user-submitted small RNAs and transcripts on psRNATarget (Dai and Zhao, 2011) to detect the sRNAs that may target the 73 genes. As most conifer genes were unannotated, we employed a reciprocal BLAST to identify homologs. Specifically, miRNA-targeted genes in P. glauca were retrieved via a BLASTN search against the Arabidopsis genome on EnsemblPlants (http://plants.ensembl.org) and then putative proteins in Arabidopsis were searched via a tBLASTN (six frames) against the P. glauca PlantGDB Putative Unique Transcripts (PUTs) database on ConGenIE (http://congenie.org/). Each PUT underwent another BLASTN search against the Arabidopsis genome. Homologs were identified only if the same pair of sequences could be found among the top three candidates in the two BLASTNs. Adaptive genes functions were classified based on the records in TAIR10 (https://www.arabidopsis.org/).

To extract expression structures of genotypes (represented by miRNA and mRNA relative expression) that can be explained by environments (temperature at seed set and phenology represented by developmental phase and pattern), redundacy analysis (RDA) was conducted (Oksanen et al., 2015) and visualized by a constrained ordination triplot with both response and explanatory variables in the same coordinate using a scaling 2 method. RDA combines multivariate regression with PCA of multiple dependent variables and we used it to test and quatify the overall contribution of a primary climate variable to the expression pattern of different types of miRNA and mRNAs.This analysis was carried out under R 3.2.2.

With high confidence, genes targeted by conserved miRNAs of interest were experimentally validated using a quantitative RT-PCR (qRT-PCR) assay as follows. Two microgram of the other aliquot of total RNAs was reverse-transcribed into cDNAs using the EasyScript PlusTM kit (abmGood) with oligo-dT primers following the manufacturer's instructions and first-strand cDNA synthesis products were diluted fivefold as qRT-PCR templates. qRT-PCR was run in 15 μl reaction volumes on an ABI StepOnePlus™ machine (Life Technologies) using the PerfeCTa® SYBR® Green SuperMix with ROX (Quanta Biosciences). The reaction components and procedure were carried out as previously described (Liu et al., 2015). We adopted three technical replicates for each pooled sample of three biological replicates. Reference genes were used as previously descibed (Czechowski et al., 2005; Liu et al., 2015). Gene homolog identifications and primer pairs for the qPCR amplification were listed in Tables S2, S3.

Small RNA libraries were prepared from four seed developmental stages of four populations in P. glauca (N.B. only three timepoints for population 1 due to its short seed developmental duration) and 10 consecutive days after pollination (0~9 DAP) of Cvi-0 and Col-0 in Arabidopsis plus another 5 days (10~14 DAP) for Cvi-0 due to its slow seed development (Figure 1). After a quality filtering, the Illumina deep sequencing yielded 15.1 M reads, on average, in 15 P. glauca libraries (Table S4). By contrast, there were on average 10.3 M reads in 25 A. thaliana libraries (Table S4). Mapping the quality-filtered reads to their corresponding genome with predictable hairpin RNA secondary structures generated an average of 3.33 M (24 ± 6%) and 7.33 M (84.3 ± 5%) reads in P. glauca and Arabidopsis libraries, respectively (Figure 2A and Table S4). Low proportion of sRNAs for P. glauca may be due to the state of scaffold-level genome assemblies in P. glauca compared with the complete assembled genome by chromosomes in the model organism, Arabidopsis. Using the archived miRNA reads under Arabidopsis species on miRBase, the known miRNAs among quality-filtered sRNA reads occupied on average 0.3% (0.05 M) and 10% (0.65 M) in P. glauca populations and Arabidopsis ecotypes over time, respectively (Table S4 and Figure 2B). This indicates that an appreciable amount of other types of sRNAs (e.g., siRNAs) is largely produced at seed set of P. glauca. The libraries of Arabidopsis and P. glauca showed, in consistency, that 24-nt sRNAs were dominantly produced at seed set (Figure 2C). In addition, sRNA reads in P. glauca classified into r/t/sn/snoRNAs took up 34.4, 5.5, 10.9, and 0.08% of the total sequences, respectively (Figure 2A); analogously, these small RNA categories were 27.3, 1.3, 4.4, and 0.1%, respectively, in Arabidopsis (Figure 2B).

After a suite of filters for P. glauca libraries, we identified 5,969 sRNA sequences, including 149 miRNAs from 62 miRNA families (Figure 3), curated in Tables S5, S6, respectively. Compared with previous reports (Källman et al., 2013; Xia et al., 2015), 90 were known miRNAs, whilst 59 were novel (Table S6). Size distribution showed that miRNAs are mainly 21-nt long throughout the plant kingdom (Figure 3-192) and identified miRNAs in spruce consistently exhibited that 21- and 22-nt were enriched in previous and this study (Figures 3-193,194). As the 21-/22-nt size class mainly consists of miRNA and phasiRNAs (or tasiRNAs), this result suggests that tasiRNAs or phasiRNAs (triggered by miRNA pairing) may play an important role in conifers, as already examined in P. abies (Xia et al., 2015). Additionally, the size distribution of all sRNAs showed that both 21- and 24-nt were abundant (Figure 3-195), indicating that hc-siRNAs for the reinforcement of silencing chromatin marks were highly generated typically at seed set in spruce. This observation from developing seeds (~30%) is significantly different from that in buds (~1%) (Källman et al., 2013), which confirms that the 24-nt long sRNA class is specific to reproductive tissues in conifers (Nystedt et al., 2013).

The top 40 conserved and differentially expressed miRNA sequences between different timepoints at seed set were respectively selected from P. glauca and Arabidopsis libraries and their expression was showcased in bivariate plots (Figure S1). In small panels of Figure S1, histograms along the diagonal were distributed in a right-skewed manner and scatter plots for P. glauca and Arabidopsis generally showed a monotonic and linear relationship, respectively. This indicates that the expression of conserved miRNAs is intrinsically correlated among different phases of seed set and populations/ecotypes.

The expression of the top 27 conserved miRNA families in P. glauca and Arabidopsis seed set was respectively employed to create a heat map and perform cluster analyses for seed-set phases and normalized miRNA expressions over time (Figures 4A,B). In P. glauca, the seed-set phases in the same population were not completely classified in one of the four major “clades” (Figure 4A). However, in Arabidopsis, the phases at 3DAP onward were neatly clustered into two primary groups by ecotype (Figure 4B). The expression pattern for Cvi-0 at flowering was different from other phases and it constituted a single “clade” (Figure 4B), while early seed-set for Cvi-0 (i.e., Cvi_1~2) and Col-0 (i.e., Col_0~2) were within the same group (Figure 4B). This indicates that the pattern of relative expression of conserved miRNAs was highly correlated with seed-set phases in Arabidopsis but not in P. glauca, and miRNAs significantly regulate seed development at least as early as 3DAP in globular embryos of Arabidopsis ecotypes. Regarding the most highly expressed and conserved miRNAs, miR166g and −319b were most abundant, followed by miR396b-5p, −156f-5p, and −167a-5p, in P. glauca (Figure 4A), while in Arabidopsis, most enriched ones were miR159b-3p, −161.1, −166a(b,e)−3p/c/d/f/g, −163, and −159c, 167a-5p/b (Figure 4B). Taken together, miR166, −319, and −396 families were highly expressed in both P. glauca and Arabidopsis and reportedly, they play regulative roles in seed/cell development (see a summary for selected miRNAs and their functions in Table S8). In general, the most abundant miRNA families (under the “clade” marked by a red dot in Figure 4) were overlapped in P. glauca and Arabidopsis but abundant miRNA families were less in number in P. glauca (Figure 4). Moreover, 18 conserved miRNA families across vascular plants were identified in both P. glauca and Arabidopsis (Figure 4C), in which some were engaged in auxin and GA signaling pathways, including miR159, −160, and −167 (Figure 4 and Table 1). In addition, there were conserved miRNAs uniquely detected at seed set in the Col-0 ecotype or in the population of late maturation (Pop 4) in P. glauca (Table 2).

Because the majority of miRNA families has a very limited taxonomic distribution (Cuperus et al., 2011), it is meaningful to carry out a lineage- and population-specific sRNA analysis. Here, lineage-specific sRNAs were represented by sRNAs identified throughout the studied populations in respective species. In P. glauca, there was no significant difference in expression abundance and pattern for lineage-specific sRNAs (p > 0.05; a total of 187 sequences used) (Figure 5IA). However, the expression pattern was significantly different between populations, when population-specific sRNAs (435, 1181, 433, and 1065 sequences for Pop 1~4, respectively) were used for comparison (all p-values < 0.05) (Figure 5IB). The sequence length peaked only in the 21-nt size class (no longer at 24-nt size; compare Figures 3-195, 5IC), supportive to miRNAs and/or tasiRNAs that were abundantly generated and may be selected for.

Analogously, there were 91 conserved miRNAs identified from sequence libraries (Table S9), in which 85 were expressed in both early (Cvi_0~7 and Col_0~4) and late (Cvi_8~14 and Col_5~9) seed-set phases (Figure S3). Conserved ath-miRNAs had a dominant length of 21 nt (Figure 5IIC) and there was no significant difference in the expression pattern for known ath-miRNAs identified across the seed-set period in both ecotypes (p = 0.786; totalizing 23 sequences used) (Figure 5IIA), while significantly different when known ath-miRNAs in either ecotype were compared (p = 0.029; 23 and 47 sequences for Cvi-0 and Col-0, respectively) (Figure 5IIB). Interestingly, there was a significant difference for sRNAs in both ecotypes (p = 0.027; totalizing 69 sequences; known ath-miRNAs excluded) (Figure 5IID) and no significant difference was found for the ones in either ecotype (p = 0.21; 85 and 99 sRNAs for Cvi-0 and Col-0, respectively) (Figure 5IIE). The size of these sRNAs peaked at the 21- and 24-nt size classes (Figure 5IIF). This indicates that hc-siRNAs, as well as tasiRNAs and novel miRNAs may have undergone selection in Arabidopsis.

In Picea glauca, considerable sRNAs were generated across seed-set phases in populations (Figure S2). The number of sRNAs detected in all four populations (1,318 reads) was as many as that of unique sRNAs in different populations (1,200 reads on average) (Figure S2A). The unique sRNAs that were expressed across developing phases in all populations were less in number than those that were only detected in one single population (Figure S2B). Burgeoning sets of sRNAs are analogous to the frequent emergence and decay previously reported in novel miRNAs in Arabidopsis (Fahlgren et al., 2007). In Arabidopsis, the number of sRNA reads without conserved miRNAs was almost 10 times as many as conserved miRNA sequences (Figure S3A). Different from conserved miRNAs that were expressed in both early and late seed-set phases, these unknown sRNAs were differently produced in phases as well as between ecotypes (Figures S3A,B). A total of 705 unknown sRNA reads out of 5,749 reads were detected in both early and late seed-set phases in the two ecotypes (Cvi_0-7 vs. Cvi_8-14 and Col_0-4 vs. Col_5-9; Figure S3B), in which 44 sequences were expressed across all seed-set phases in two ecotypes.

In light of the tight correlation between sRNA conservation and expression abundance (Chávez Montes et al., 2014), GO classification for target genes of the 107 “most conserved” sRNAs in P. glauca (Table S7) showed that sRNAs were mostly involved in the regulation of pathways related with metabolic and cellular processes (Figure 6A), which was in line with the observation in conserved miRNAs (Table S9) in Arabidopsis (Figure 6B). The experimentally validated miRNA-gene pairs were listed in Figure 6C and their target genes are mainly involved in seed development (Figure 6C).

We adopted an Illumina sequencing approach to get a deep coverage of mature sRNAs and predicted sRNA candidates that may target the reported genes involved in adaptation to climate in P. glauca. We found 24 sRNA candidates (Table 3) targeting genes putatively functioning mainly in abiotic and biotic stresses (Figure 7A). In model organisms, a group of key and conserved miRNA sequences involved in plant seed development and phase transitions was summarized in Table S8, some of which are involved in adaptation to stresses (see description in the “function” column). They were identified by fold changes in miRNAs between stress-treated and control samples. With reference to previous reports, we identified 9 miRNAs involving in adaptation during Arabidopsis seed development and they were miR159, −167, −169, −393, −395, −397, −398, −399, and −408 (Figure 7B). These miRNA families also have functionalities in seed development via regulating transcriptional factors in hormone-based GRNs, such as miR159, −167, −393 targeting GAMYB, AUXIN RESPONSE FACTORS (ARFs), auxin-receptors and cyclin-like F-box, respectively (see a summary in Table S8).

We chose the seed dormancy phenotype to investigate how environmental signals impinge on the expression of conserved miRNAs targeting genes related to seed dormancy and also that of key conservative genes conducive to the phenotype, thus collectively manipulating phenotypic variation. Gene phylogeny for ARF10/16 showed that gymnosperm and model angiosperm species were separated into different clades except for Picea, which has a higher sequence similarity with angiosperms than species in the same taxonomic category (Figure S4). This indicates that ARF10/16 is ancient and evolutionarily conserved within the plant kingdom. Conserved domain analysis showed that putative ARF10 in P. glauca harbors an Aux_IAA super family domain (Figure S5). The pivotal activation function of ARF proteins is conferred by their four-domain architecture, including DNA binding region (a B3 and an ARF domain) and protein dimerization motifs (Ulmasov et al., 1999; Tiwari et al., 2003). Loss of the canonical four-domain structure has promoted functional shifts within the ARF family by disrupting either dimerization or DNA-binding capacities (Finet et al., 2013). As such, the putative ARF10 in P. glauca may not correspond to its counterpart in Arabidopsis or the essential Aux_IAA domain is sufficient to have the function of ARF10 in P. glauca. No homolog hits for ABI3 and DOG1 were obtained in some species (Table S2) and the phylogeny showed that they might have undergone substantial selection (Figure S4). The miR160 targets ARF10 and its hairpin structure was computationally predicted (Figure S6).

The relative expression of aforementioned genes was displayed in Figure S7. As predicted, AtDOG1, AtABI3, AtARF10, and AtARF16 were highly expressed in Cvi-0 than Col-0 (Figure S7, left panels). However, such pattern was not observed for gene counterparts in populations of P. glauca (Figure S7, first three panels on the right side), indicating that these genes in P. glauca may have different functions or that the studied phenotype in P. glauca is regulated by other unknown mechanisms. In the RDA triplot, the percentage of accumulated constrained eigenvalues showed that the first axis explained 43.8% variance (Figure S7, last panel), indicating that the major trends have been modeled by RDA. In addition to species and dev_phase, developmental temperature played an important role in the dispersion of developmental phases along the first axis and had high correlation with miR160 and ARF10/16 (Figure S7, last panel). As transcripts of ARF10/16 are targeted by miR160 (Liu et al., 2013), the expression patterns of ARF10/16 and miR160 were highly positively correlated with each other but negatively correlated with that of DOG1 (Figure S7). In addition, the projection of the same developmental phases on the first axis was overlapped across populations in P. glauca (Figure S7). This suggests that the same phase between populations has more similarities in gene expression pattern than different phases within populations, and in turn, prompts the conservation of gene regulation at a temporal scale across populations.

Of 37 miRNA families that are deeply conserved in plant development throughout the plant kingdom (Willmann and Poethig, 2007), we identified 12 miRNA families at seed set between phyla (Table 1). Some are conserved across spermatophytes (i.e., miR163, −394, and −396), tracheophytes (i.e., miR159), and embryophytes (i.e., miR160, −166, −167, −171, −319, −390, and −408). The ubiquitously conserved miRNAs had significantly differential expression abundance between P. glauca and Arabidopsis as well as between populations within species (e.g., Figure 5IA vs. Figure 5IIA), and this observation is correlated to genome evolution (Hodgins et al., 2016). These mature miRNA sequences have a well-conserved consensus with minor variation at the 5′ or 3′ end between P. glauca and Arabidopsis (Table 1), and their targets may not always be conserved partially due to no homologs found in P. glauca (Table 1). This indicates that the cognate miRNA-target pairs may be acquired prior to the split of gymnosperms and angiosperms and have undergone natural selection. Conserved miRNAs in Arabidopsis (Table S8) were enriched and expressed across seed set (Figure 5IIA), while abundant sRNAs expressed throughout seed set in P. glauca did not belong to miRNA families documented on miRBase (Table S7). Nonetheless, their target genes were classified into the same GO categories (Figure 6). This indicates that spruce and Arabidopsis may use different sRNAs to attain similar regulations at molecular levels. Note that the known miRNA-gene pairs in Arabidopsis are associated with seed development and nodes in hormone-based GRNs (Figure 6C).

In general, ancient miRNAs are more highly and broadly expressed than younger ones (Fahlgren et al., 2007); while newly emerged miRNAs may be used as substrates for natural selection and form specific miRNA circuitry interplaying with GRNs, whereby adaptation and speciation occur (Chen and Rajewsky, 2007). During evolution, deleterious miRNA-target pairs are selected against and miRNAs by gradually decreasing the number of target genes are retained into GRNs (Chen and Rajewsky, 2007). Through constraining variance or mean gene expression, miRNAs render phenotypic traits evolvable as well as heritable (Wagner and Altenberg, 1996) and after selection, these miRNAs may improve fitness of phenotypes (Wu et al., 2009). However, most of young miRNAs appear to be selectively neutral [S = 0] (Fahlgren et al., 2010), because the majority of them does not have target genes or does not regulate targets in such a way that plant fitness can be enhanced, and they were quickly lost during evolution. In regulation of target genes with similar functions, different species may rely on variants of the same MIRNA families or produce distinct miRNAs through changes in miRNA precursor sequences and/or binding sites. From this study, the conserved miRNAs throughout the plant kingdom and lineage-specific sRNAs (Figures 4, 5) represent different levels of sRNA conservation serving for basic needs in plants and particular requirements in a given lineage, respectively. On the other hand, although members of plant MIRNA families are often highly similar, the same gene family varies in size and genomic structure between species, indicating dosage effects and spatiotemporal differences in MIRNA regulation (Li and Mao, 2007). There are miRNA orthologs that function quite differently and several miRNA isoforms have specific tissue expression patterns (Lelandais-Brière et al., 2009), indicative of functional divergence. We detected some miRNAs solely expressed in one population/ecotype, which were conserved in miRNA families on miRBase (Table 2). As per target gene function, target genes do not seem to be key components in GRNs and they primarily target repeats, intergenic sequences, and genes in cell-cell communication and signaling (Table 2). This suggests that altering mean miRNA expression in different plant populations is subject to fine-tuning to regulate the expression of target genes, especially for non-key ones in GRNs. The dosage balance hypothesis, a mechanism in maintaining MIRNA duplicates after WGDs (possibly also in small-scale tandem or segmental duplication events due to little impact on maintaining optimal stoichiometry among gene products; e.g., Birchler and Veitia, 2012), supports our inference, as key genes in GRNs are dosage sensitive, such that deletion of one copy may disrupt the whole network.

As requirements for a functional MIRNA are less demanding than a protein-encoding gene (PEG), MIRNA can easily evolve from various sources of unstructured transcripts, such as gene duplication, intergenic regions, transposable elements (Nozawa et al., 2012). As well, like PEGs, MIRNAs are subject to the same evolutionary processes, such as, substitution, insertion, recombination, and natural selection. There are three proposed models that explain the origination of MIRNAs, that is, transposable elements, inverted duplication of target genes, and random hairpin sequences with subsequent mutations, reviewed by Cui et al. (2016). A study in spruce supports the second model, evidenced by the precursor sequence of a copy of miR482 highly similar to nucleotide-binding site and leucine-rich repeat domains (NBS-LRR) genes (Xia et al., 2015). The third model emphasizes that a serendipitous DNA sequence transcribed and Dicer-processed yields a sRNA that interacts with a homolog target. Coevolution of miRNA sequence and target genes refines such interaction (Felippes et al., 2008). The miRNA-target complementarity in plants is extensive and stringent (Rhoades et al., 2002), suggestive of coevolution of miRNAs and their targets. Taken together, miRNA-target systems have undergone evolution in the context of genome evolution; despite divergences in miRNA sequences, the function of their target genes is largely of commonality (e.g., indirect evidence from Figure 6) and this supports similar molecular mechanisms and signaling cascades at seed set across phyla.

To cope with the vagaries of environmental perturbations, plants have evolved mechanisms of stress avoidance (acclimation) and stress tolerance (adaptation) via well-honing gene expression. Reprogramming of gene expression via sRNAs is a major defense mechanism in plants as response to stresses (see a list of references in the Figure 7 legend). The sRNA pathways through (P)TGS intersect with mechanisms regulating different steps in the life of an mRNA, starting from transcription and ending at mRNA decay, in both nuclear and cytoplasmic compartments. Stresses usually trigger destabilization of chromatin states, manifested as changes in DNA methylation and reductions in nucleosome occupancy, and some of these processes involve miRNAs via the RdDM pathway and TE-associated 21-nt siRNAs (Dowen et al., 2012). Recently, a type of stress-induced, UTR-derived siRNAs in 24-nt long has been identified in Brachypodium and proposed to mask cis-elements (i.e., pre-mRNA) and direct the splicing machinery to use the correct splice site (Wang et al., 2015). Under stress, alternative splicing results in production of different splice variants and thus proteins with different functions required for survival, reviewed by Mastrangelo et al. (2012).

Studies in genomic basis of adaptation in long-lived organisms have largely focused on identifying intraspecific DNA variation and comparative gene expression that are associated with adaptation, reviewed by Prunier et al. (2016). In boreal and temperate regions, conifers with complex life cycles have been found to rely on the regulation of certain gene expression to defense against climatic stresses (Hornoy et al., 2015; Yeaman et al., 2016). These adaptive genes mainly participate in the process against abiotic stresses (e.g., cold injury, resistance to aridity, high salinity, UV-B, etc.) and fungal pathogens, growth cessation and bud set to synchronize with seasonal climatic changes. As such, divergence in gene expression possibly through stabilizing or directional selection (because studies in animals indicate that most gene expression variability may be deleterious and selected against [S = −∞]) is important to species adaptability and evolution. However, the long generation time (i.e., old age of sexual maturity) imposes limitations on natural selection and adaptive evolution for genes to cope with rapidly changing climate. Consequently, these organisms confront a real need to use epigenetics, owing to its merits of high flexibility, to deal with a wide range of environmental conditions throughout their life time. Our reported sRNAs (Figure 7A and Table 3) have opened a venue into better understanding evolutionary trajectories among small RNAs targeting adaptive genes (albeit not experimentally confirmed), as natural selection of sRNAs in different lineages and among populations within lineage may exhibit different evolutionary patterns and delving into these patterns for different organismal complexity helps uncover sRNA evolution in line with the complexity of organisms.

The phenology or temporal control of the life cycle provides adaptive strategies to avoid adverse consequences in harsh environments at seed set and seedling establishment (Krämer, 2015). In the life cycle, temperature is of utmost importance in seed set, germination, and seedling stage due to epigenetic imprinting and their fragile state (Liu et al., 2016). At seed set, temperature signals are a critical selective pressure and have a strong influence on life history traits, such as timing of seed set and seed dormancy depth (Springthorpe and Penfield, 2015; Vidigal et al., 2016). Specifically, the temperature-mediated control of flowering has evolved to constrain the maternal environment for setting seeds to a specific temperature window, thus yielding seeds with dormancy variation. Maternal environments (e.g., temperature cues) have a persistent and transgenerational effect on the expression pattern of regulatory molecules in GRNs (e.g., DOG1, Chiang et al., 2013) and on the biogenesis of versatile miRNAs in life histories. The latter was confirmed by miRNAs (i.e., miR160) targeting key genes (i.e., ABI3 and ARF10/16) that modulate timing of seed set and the dormancy state (Huo et al., 2016). In this study, the chosen populations of P. glauca and ecotypes of Arabidopsis are characterized by different fertilization timepoints, contrasting seed-set durations and ensuing phenotypic variation (i.e., dormancy intensities; Figure 1). We tested how the role of miRNA-gene interactions is at play for one study phenotype, seed dormancy, in which we chose a miRNA-gene pair (i.e., miR160-ARF10/16) and key genes (i.e., ABI3 and DOG1) (Figure S7). Although the chosen gene may have different functions in P. glauca, our results reinforce the notion that environment factors, such as temperature have great impact on the expression of genotypes (miRNAs and genes).

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 RN and handling Editor declared their shared affiliation.