Integrated sRNA-seq and RNA-seq Analyses Reveal a microRNA Regulation Network Involved in Cold Response in Pisum sativum L.

(1) Background: Cold stress affects growth and development in plants and is a major environmental factor that decreases productivity. Over the past two decades, the advent of next generation sequencing (NGS) technologies has opened new opportunities to understand the molecular bases of stress resistance by enabling the detection of weakly expressed transcripts and the identification of regulatory RNAs of gene expression, including microRNAs (miRNAs). (2) Methods: In this study, we performed time series sRNA and mRNA sequencing experiments on two pea (Pisum sativum L., Ps) lines, Champagne frost-tolerant and Térèse frost-sensitive, during a low temperature treatment versus a control condition. (3) Results: An integrative analysis led to the identification of 136 miRNAs and a regulation network composed of 39 miRNA/mRNA target pairs with discordant expression patterns. (4) Conclusions: Our findings indicate that the cold response in pea involves 11 miRNA families as well as their target genes related to antioxidative and multi-stress defense mechanisms and cell wall biosynthesis.


Introduction
Grain legumes produce seeds that are a valuable source of protein and starch for both animal and human diets [1,2]. Additionally, these crops can establish a symbiosis with nitrogen-fixing soil bacteria, allowing a reduction in the use of nitrogen fertilizers in cropping systems and thus N 2 O emissions [3]. Due to these agronomical and ecological benefits, legumes are considered to be a key component of sustainable agriculture [4,5]. Among legumes, dry peas (Pisum sativum L., Ps) are widely grown in cool temperate areas. Although they are usually sown in spring in Northern Europe, autumn sowings (winter peas) are desirable as they allow for an increased and stabilized seed yield, as a result of a longer life cycle and the possibility of avoiding frequent drought during the reproductive period [6]. Frost is the main abiotic stress that must be overcome during the winter period. Thus, breeding for winter varieties has motivated several studies aiming at localizing the genetic determinants of frost tolerance and understanding the molecular bases of this trait [7][8][9][10][11]. In these studies, special attention has been given to the ability of frost tolerant genotypes to acclimatize to cold, that is to improve their frost tolerance level in response to the low above-zero temperatures generally encountered in autumn in field conditions. Understanding the mechanisms that regulate the cold acclimation process is of growing Total RNAs from each sample were sent to Eurofins MWG GmbH where small cDNA libraries with insert size of 19-30 bp were prepared. Subsequently, the 24 bar-coded libraries were sequenced with an Illumina HiSeq 2000/2500 machine, leading to a total of approximately 380 million single-end 50 bp reads. Raw sRNA-seq data were transferred to the NCBI-SRA database, in the BioProject PRJNA543764 [23], which also includes RNA-seq data previously deposited.

Identification of miRNAs and Their Targets
Prior to miRNA detection using miRkwood (v1.0; [28]), sRNA sequences were mapped without mismatch on the Cameor pea genome [22] and on the Ch/Te pea transcriptome [11] using Bowtie (-v 0 --all --best --strata). miRkwood was used with default parameters and only predictions with a score ≥5 (maximum of 6) were kept. The quality score provided by miRkwood (from 0 to 6) arises from the 6 criteria used to detect miRNAs. The first criterion is related to the stability of the stem-loop structure of the precursor, whereas the five others are related to the distribution of mapped reads along the putative pre-miRNA. They allow us to determine whether this distribution presents a typical profile with 2 peaks, corresponding to the guide and the passenger miRNAs. Candidates with a score of 6 fulfill all criteria, whereas those with a quality score of 5 fail a single criterion. miRNA predictions obtained using separate analyses for Ch and Te samples by aligning reads on either the pea genome or the pea reference transcriptome were then merged. Conserved miRNAs were identified by comparing miRNA predicted sequences to the 10,414 mature plant (Viridiplantae) miRNAs deposited in miRBase (v22; [29]) using Exonerate (v.2.2.0; [30]) selecting only alignments with at most three mismatches. Finally, only miRNA predictions with a miRkwood score of 6 or those with both a score of 5 and showing alignments with miRNA sequences from miRBase were considered for further analyses. miRNA sequences were additionally compared to the 68 miRNAs identified by Kreplak et al. [22] in developing seeds; only exact matches between miRNA sequences were considered.
For every conserved and non-conserved miRNA, targets were predicted using Tar-getFinder (v1.7; [31]) with default parameters and the pea predicted coding sequences as the target sequence database.

RNA-seq Data Processing
For each RNA-seq library from Bahrman et al. [11], the overall quality of reads was checked using FastQC (v0.11.4; [26]). Adapter sequences were removed using Cutadapt (v1.0; [24]). For every sequenced fragment (i.e., pair of reads), overlapping ends were merged into a single sequence and nucleotides with the lowest quality score were replaced by an N until the overall error rate of the read dropped below 2%. N at read ends were also discarded.

Differential miRNA and Gene Expression Analysis
Differential expression analyses were performed with the R (v4.0.3)/Bioconductor (v3.12) package edgeR (v3.32.1; [36]). Non-uniquely mapped reads and genes with less than one count per million in at least half of the twenty-four libraries were filtered out. mRNA and miRNA library sizes were normalized using the TMM approach [37]. Differentially expressed genes (DEGs) and differentially expressed miRNAs were identified by a negative binomial generalized linear model (GLM) using a group-means parameterization (where each combination of genotype, treatment, and time is defined as a group) and the Cox-Reid profile-adjusted likelihood method for estimating dispersions. Comparisons of interest were then estimated using contrasts, using a false discovery rate (FDR) adjusted p-value < 0.05 as a significance threshold. T1 and T2 DEGs and differentially expressed miRNAs were identified within each time point by comparing the LT to the N condition for each genotype (Ch or Te) (Additional File S1).
For both periods, we identified four sets of DEGs and differentially expressed miRNAs: (1) between LT and N in Ch at T1; (2) between LT and N in Ch at T2; (3) between LT and N in Te at T1; and (4) between LT and N in Te at T2. These sets were further used to pinpoint Ch-and Te-specific cold responsive genes and miRNAs as well as those involved in cold response common to both genotypes. T0 was treated as a control time point, since environmental conditions for both the LT and N experiments were identical until this sampling date. As such, mRNAs or miRNAs with unexpected differential expression between the LT and N conditions at T0 were removed from the lists of differential elements of interest.

Gene Ontology Enrichment Analysis
The GO enrichment analyses for DEGs and miRNA targets were performed with the R (v4.0.3)/Bioconductor (v3.12) package topGO (v2.42.0; [38]). The elim algorithm with a Fisher test and a nodeSize at 5 was used to perform the analysis. The GO annotation file described in [22] was used to build the mapping between GO terms and Cameor genes. In the visual representation of enrichment analysis results, GO terms were clustered according to their semantic similarity calculated with GOGO [39] using hierarchical clustering with complete linkage.

Small RNA Analysis
High-throughput Illumina sequencing yielded between 9.9 and 27.6 million reads for each library, for a total of approximately 380 million reads. After removing adapters and filtering out low quality sequences, 123 million clean 18-25 nt reads were counted, representing 20 million unique sequences. The majority of the sRNAs were 21-24 nt in size (94%), with 21 and 24 nt being the most frequently represented for all reads and unique reads (without reads redundancy), respectively ( Figure 1a). These results were consistent with the representative size range of Dicer-like (DCL) cleavage products [40] and with the typical small RNA distribution of higher plants (e.g., in At, Medicago truncatula and Paeonia suffruticosa [41][42][43][44], supporting the reliability of the Illumina sequencing data for these sRNA libraries.
Among the 136 miRNA predictions, 102 (75%) presented alignments with known miRNAs from miRBase belonging to 28 distinct families (+11 that were not assigned to a miRNA family) ( Table S1). The three most represented families were miR156, miR159 and miR166, with 12, 10 and 13 miRNAs, respectively. Seven of the 28 families represented in our dataset (miR156, miR159, miR160, miR166, miR171_1, miR390, and miR396) are highly conserved among plants. Indeed, they are represented in each plant lineage, from mosses to higher plants, for which miRNAs have been deposited into miRBase. The miRNAs belonging to these conserved families were represented by 5 million reads on average (Additional File S3). On the contrary, several miRNAs (including families miR169_6, miR862_2, miR1507, miR5037 and miR5770) that were retrieved only in Fabaceae species were represented by a much smaller number of reads, on average 13,336 (mean comparison using Mann-Whitney test: p = 0.001892). These results are consistent with previous studies indicating that conserved miRNAs are highly expressed and that most have several family members, whereas non-conserved ones are often weakly expressed and encoded by single loci [45].
The predicted miRNAs were also compared with the set of 68 miRNAs identified in pea developing seeds by Kreplak et al. [22]. Among the 68 miRNAs expressed in seeds, 47 were found in our analysis, indicating they are expressed in both seeds and aerial parts.

miRNAs Involved in Cold Response
Based on a comparison between the samples subjected to a cold acclimation treatment (4.5 °C; LT) or not (N) before applying frost, the differential analysis allowed the identification of 59 miRNAs involved in cold response (43.3% of the miRNAs) (Tables S2-S7, Additional File S4). Among them, 40 miRNAs presented a differential expression (LT vs. N) in only one of the two genotypes (i.e., genotype-specific responses) at T1 and/or T2, 16 presented a differential expression in the two genotypes (common responses) at T1 and/or T2, and 3 were found to be related to both specific and common responses at different time points (Figure 2).
Among the 136 miRNA predictions, 102 (75%) presented alignments with known miRNAs from miRBase belonging to 28 distinct families (+11 that were not assigned to a miRNA family) ( Table S1). The three most represented families were miR156, miR159 and miR166, with 12, 10 and 13 miRNAs, respectively. Seven of the 28 families represented in our dataset (miR156, miR159, miR160, miR166, miR171_1, miR390, and miR396) are highly conserved among plants. Indeed, they are represented in each plant lineage, from mosses to higher plants, for which miRNAs have been deposited into miRBase. The miRNAs belonging to these conserved families were represented by 5 million reads on average (Additional File S3). On the contrary, several miRNAs (including families miR169_6, miR862_2, miR1507, miR5037 and miR5770) that were retrieved only in Fabaceae species were represented by a much smaller number of reads, on average 13,336 (mean comparison using Mann-Whitney test: p = 0.001892). These results are consistent with previous studies indicating that conserved miRNAs are highly expressed and that most have several family members, whereas non-conserved ones are often weakly expressed and encoded by single loci [45].
The predicted miRNAs were also compared with the set of 68 miRNAs identified in pea developing seeds by Kreplak et al. [22]. Among the 68 miRNAs expressed in seeds, 47 were found in our analysis, indicating they are expressed in both seeds and aerial parts.

miRNAs Involved in Cold Response
Based on a comparison between the samples subjected to a cold acclimation treatment (4.5 • C; LT) or not (N) before applying frost, the differential analysis allowed the identification of 59 miRNAs involved in cold response (43.3% of the miRNAs) (Tables S2-S7, Additional File S4). Among them, 40 miRNAs presented a differential expression (LT vs. N) in only one of the two genotypes (i.e., genotype-specific responses) at T1 and/or T2, 16 presented a differential expression in the two genotypes (common responses) at T1 and/or T2, and 3 were found to be related to both specific and common responses at different time points (Figure 2). ily gathered 27% of the differentially abundant miRNAs. Four miRNA families were equally represented within Te-specific miRNAs: members of the miR164 and miR167_1 families were up-regulated at T1 and T2, a member of the miR159 family was down-regulated at T1 and up-regulated at T2, and a member of miR156 was up-regulated at T1 and down-regulated at T2. Among miRNAs related to the common response of the two genotypes, a member of the miR156 family was up-regulated at T1 and down-regulated at T2, whereas one of the miR166 family was up-regulated at T2.

Identification of miRNA Targets
We then predicted putative target genes for the 136 predicted miRNAs by looking for sequence complementarity with the coding sequences extracted from the pea genome annotation [22]. Among the 136 previously selected miRNAs, 106 (77.9%) were found to have putative targets (1789 transcripts representing 1253 distinct genes, Table S8). A Gene Ontology (GO) analysis performed on the set of target genes pointed out an enrichment of biological processes (BP) related to development and defense ( Figure 3). Among miRNAs presenting a modulation of expression only in Ch, the miR159 family gathered 27% of the differentially abundant miRNAs. Four miRNA families were equally represented within Te-specific miRNAs: members of the miR164 and miR167_1 families were up-regulated at T1 and T2, a member of the miR159 family was down-regulated at T1 and up-regulated at T2, and a member of miR156 was up-regulated at T1 and downregulated at T2. Among miRNAs related to the common response of the two genotypes, a member of the miR156 family was up-regulated at T1 and down-regulated at T2, whereas one of the miR166 family was up-regulated at T2.

Identification of miRNA Targets
We then predicted putative target genes for the 136 predicted miRNAs by looking for sequence complementarity with the coding sequences extracted from the pea genome annotation [22]. Among the 136 previously selected miRNAs, 106 (77.9%) were found to have putative targets (1789 transcripts representing 1253 distinct genes, Table S8). A Gene Ontology (GO) analysis performed on the set of target genes pointed out an enrichment of biological processes (BP) related to development and defense ( Figure 3).

Differentially Expressed Genes in Response to Cold
The recent release of a pea reference genome (cv Cameor, [21]) enabled us to perform a new assessment of gene expression using RNA-seq data previously produced by Bahrman et al. [11] thanks to a pipeline similar to that used for the miRNA-seq data. This updated gene expression quantification facilitated comparisons and integration of RNA-seq and miRNA-seq results. Illumina sequencing of the twenty-four mRNA libraries previously produced by Bahrman et al. [11] yielded a total of 420 million pairs of reads. After filtering out low quality reads, more than 99.99% of these were aligned to the Cameor genome sequence. 393 million (94.3%) pairs could be aligned to the genome assembly, including 358 million with a single genomic location. 94% and 1% of these alignments corresponded to the coding and non-coding strands of annotated genes, respectively, and the remaining 5% were located outside any annotated gene.

Differentially Expressed Genes in Response to Cold
The recent release of a pea reference genome (cv Cameor, [21]) enabled us to perform a new assessment of gene expression using RNA-seq data previously produced by Bahrman et al. [11] thanks to a pipeline similar to that used for the miRNA-seq data. This updated gene expression quantification facilitated comparisons and integration of RNAseq and miRNA-seq results. Illumina sequencing of the twenty-four mRNA libraries previously produced by Bahrman et al. [11] yielded a total of 420 million pairs of reads. After filtering out low quality reads, more than 99.99% of these were aligned to the Cameor genome sequence. 393 million (94.3%) pairs could be aligned to the genome assembly, including 358 million with a single genomic location. 94% and 1% of these alignments corresponded to the coding and non-coding strands of annotated genes, respectively, and the remaining 5% were located outside any annotated gene.
The transcriptomic responses of the two contrasted pea lines (Ch, frost tolerant; Te, frost susceptible) subjected to a cold acclimation treatment (4.5 °C; LT) or not (N) before applying frost was studied after different periods of cold exposure T1 and T2 [11]. For The transcriptomic responses of the two contrasted pea lines (Ch, frost tolerant; Te, frost susceptible) subjected to a cold acclimation treatment (4.5 • C; LT) or not (N) before applying frost was studied after different periods of cold exposure T1 and T2 [11]. For both periods, we identified four sets of DEGs: (1) between LT and N in Ch at T1 (Table S9); (2) between LT and N in Ch at T2 (Table S10); (3) between LT and N in Te at T1 (Table S11); and (4) between LT and N in Te at T2 (Table S12). These sets were then used to pinpoint Ch and Te-specific cold responsive genes (Tables S9-S12) and genes involved in cold response common to both genotypes (Tables S13 and S14). Among the various comparisons considered, we identified between two and four thousand differentially expressed genes under cold stress (Table 1). A total of 3663 and 3666 DEGs at T1 and T2, respectively, were involved in cold response in Ch and Te. A small number of genes was found to be involved in cold response at T1 and T2 in both genotypes but with discordant differential expressions. Some DEGs were found to have a consistent response regardless of the cold exposure period: 459 DEGs for the Ch-specific response, 397 DEGs for the Te-specific response, and 1610 DEGs for the common response to cold. We also observed cases of offsets in cold response between the two genotypes: for example, 220 DEGs in Te at T1 were also differentially expressed in Ch at T2, and 331 DEGs in Ch at T1 were differently expressed in Te at T2 (Figure 4). Table 1. Summary of DEGs identified in Champagne (Ch) and Térèse (Te) in response to cold in time periods T1 and T2 (adjusted p-values < 5%).

List of DEGs DEGs at T1 DEGs at T2
Te specific cold responsive genes S9); (2) between LT and N in Ch at T2 (Table S10); (3) between LT and N in Te at T1 (Table  S11); and (4) between LT and N in Te at T2 (Table S12). These sets were then used to pinpoint Ch and Te-specific cold responsive genes (Tables S9-S12) and genes involved in cold response common to both genotypes (Tables S13 and S14). Among the various comparisons considered, we identified between two and four thousand differentially expressed genes under cold stress (Table 1). A total of 3663 and 3666 DEGs at T1 and T2, respectively, were involved in cold response in Ch and Te. A small number of genes was found to be involved in cold response at T1 and T2 in both genotypes but with discordant differential expressions. Some DEGs were found to have a consistent response regardless of the cold exposure period: 459 DEGs for the Ch-specific response, 397 DEGs for the Tespecific response, and 1610 DEGs for the common response to cold. We also observed cases of offsets in cold response between the two genotypes: for example, 220 DEGs in Te at T1 were also differentially expressed in Ch at T2, and 331 DEGs in Ch at T1 were differently expressed in Te at T2 (Figure 4).

List of DEGs DEGs at T1 DEGs at T2
Te specific cold responsive genes  To investigate the function of cold-responsive DEGs, a GO enrichment analysis was performed on three separate gene lists ( Figure 5, Tables S15-S17). To constitute the three lists, we first combined the results for the two different periods of cold exposure (T1 and T2); we then focused on Te-specific, Ch-specific, and common concordant cold-responsive genes. Signal transduction, protein dephosphorylation and response to light stimulus were the most enriched biological processes for DEGs common to both genotypes. For Ch-specific DEGs, photosynthesis, protein folding and metal ion transport were most highly enriched, while cellular modified amino acid and cellular glucan metabolic processes were most highly enriched for Te-specific DEGs. Based on a p-value significance threshold of 1%, no biological processes were found to be enriched across all three lists.
To investigate the function of cold-responsive DEGs, a GO enrichment analysis was performed on three separate gene lists ( Figure 5, Tables S15-S17). To constitute the three lists, we first combined the results for the two different periods of cold exposure (T1 and T2); we then focused on Te-specific, Ch-specific, and common concordant cold-responsive genes. Signal transduction, protein dephosphorylation and response to light stimulus were the most enriched biological processes for DEGs common to both genotypes. For Chspecific DEGs, photosynthesis, protein folding and metal ion transport were most highly enriched, while cellular modified amino acid and cellular glucan metabolic processes were most highly enriched for Te-specific DEGs. Based on a p-value significance threshold of 1%, no biological processes were found to be enriched across all three lists.  . Significantly (p-value < 0.01) enriched Gene Ontology (GO) biological process terms among differentially expressed genes (DEGs) specific to Ch, Te, and common to the two genotypes. Counts indicate the number of DEGs annotated with each GO term, and dots are colored by p-value. GO terms are clustered using hierarchical clustering (complete linkage) using GOGO semantic similarity [38] between terms.

miRNAs and Target mRNAs with Discordant Expression Profiles
Among the 1253 miRNA target genes identified, 252 were included in the lists of DEGs. In total, 84 miRNA/target DEG pairs were identified, with 39 of them presenting discordant expression patterns ( Figure 6). For Ch-specific miRNA/mRNA pairs at T1, four miRNAs were up-regulated and six of their targets were down-regulated, including three multicopper oxidase. No miRNA/mRNA pairs were identified for Ch-specific genes and miRNAs involved in cold response at T2. For Te-specific cold response at T1, two miRNAs were up regulated, one was down regulated. At T2, six pairs including five target genes were found. In this case, 15 pairs of miRNAs/mRNAs were identified with a cold response common to the two pea genotypes. Interestingly, four miRNAs were up-regulated at T1 but down-regulated at T2 with a unique and different target DEG at T2, suggesting they have different targets at different time points. miRNAs were up-regulated and six of their targets were down-regulated, including three multicopper oxidase. No miRNA/mRNA pairs were identified for Ch-specific genes and miRNAs involved in cold response at T2. For Te-specific cold response at T1, two miRNAs were up regulated, one was down regulated. At T2, six pairs including five target genes were found. In this case, 15 pairs of miRNAs/mRNAs were identified with a cold response common to the two pea genotypes. Interestingly, four miRNAs were up-regulated at T1 but down-regulated at T2 with a unique and different target DEG at T2, suggesting they have different targets at different time points. Figure 6. Discordant patterns of significant differential expression for miRNAs and their putative mRNA targets. Ch, Te and Com, respectively represent Ch-specific, Te-specific and common differential cold-response. Down-and up-regulated mRNAs and miRNAs are shown in blue and red, respectively. Lines joining miRNAs and mRNAs represent putative miRNA-target relationships.

Metabolic Responses to Cold in Pea
The GO enrichment analysis we present here provides useful insight into the pea response to cold stress ( Figure 5). For common differential cold-response in gene expression, signal transduction and protein dephosphorylation are the most highly enriched biological processes. Signal transduction and protein dephosphorylation have a key role in miRNA family miRNA  Figure 6. Discordant patterns of significant differential expression for miRNAs and their putative mRNA targets. Ch, Te and Com, respectively represent Ch-specific, Te-specific and common differential cold-response. Down-and up-regulated mRNAs and miRNAs are shown in blue and red, respectively. Lines joining miRNAs and mRNAs represent putative miRNA-target relationships.

Metabolic Responses to Cold in Pea
The GO enrichment analysis we present here provides useful insight into the pea response to cold stress ( Figure 5). For common differential cold-response in gene expression, signal transduction and protein dephosphorylation are the most highly enriched biological processes. Signal transduction and protein dephosphorylation have a key role in plant responses to environmental stresses such as drought, high salinity, cold, and pathogen attack. Kinase proteins are known as a key feature of signal transduction [46]. In this study, 177 DEGs (4.83%) and 199 DEGs (5.43%) encode protein kinases at T1 and T2, respectively, and 75 of these DEGs are common between T1 and T2. In the common differential coldresponse, trehalose biosynthetic process is also a significantly enriched term. Trehalose is known as a membrane stabilizer and an inhibitor of membrane fusion [47]. For example, seven genes involved in trehalose synthesis are reported to be more expressed in At under cold exposure [48,49]. Polyamines can also stabilize the membranes and inhibit protein and mRNA denaturation. In the common differential cold response, several biological processes are linked to polyamines: amine catabolic process, amine transport, methionine metabolic process and sulfur amino acid biosynthetic process. Methionine and cysteine are substrates for the synthesis of various polyamines with important roles in stress tolerance, the most prominent being putrescine, spermidine and spermine [50]. Comparing two chickpea (Cicer arietinum) genotypes, higher levels of putrescine, spermidine and spermine were observed after six days of cold stress at 4 • C [51].
Metal ion transport, which appears in Ch-specific cold response ( Figure 5), could be linked to the photosynthesis changes during cold stress. Fe, Cu, and Mn readily change their oxidative state, which enables these transition metals to participate in vital cellular processes such as photosynthesis or the electron transport chain [52]. In Ch, a higher inherent photosynthetic potential at the beginning of the cold exposure, combined with an early ability to start metabolic processes aimed at maintaining the photosynthetic capacity, has already been described [9]. This ability to start metabolic processes in Ch is also observed via other biological processes: secondary metabolic process, glycerol ether metabolic process or protein folding. Protein folding in vivo is mediated by an array of proteins that act as molecular chaperones. Molecular chaperones share the property that they bind substrate proteins that are in unstable, non-native structural states [53]. Among the chaperone proteins, heat shock proteins (HSPs) are known to enhance stress tolerance such as temperature shift. In At, four genes encoding HSP70 proteins are up-regulated after 48 h at 4 • C [54]. In Nt, 10 HSP70 genes out of 61 are up-regulated during cold treatment [55]. In Ch, three HSP70 genes are up-regulated: Psat7g103280 and Psat3g071280 at T1, Psat7g103280 and Psat1g222760 at T2.
Several biological processes involved in Ch-specific cold response and common differential cold-response are linked to Te-specific cold response biological processes. Similar to trehalose, glucans stabilize the membranes and inhibit membrane fusion [56]. As in Ch, cold stress induces photosynthetic changes in Te with a modification of the electron transport in photosystem I. In Te, the cytoskeleton undergoes modification during exposure to low temperature with a significant enrichment of the GO term actin filament depolymerization. It has been shown that the rearrangement of the cytoskeleton is also linked to the induction of genes related to cold acclimation in At cell suspensions and Brassica napus leaves [57,58]. In Nt cell suspensions at 0 • C, the microtubules and filaments form bundles [59] and, in Vitis rupestris cell suspensions at 0 • C, microtubules totally disassemble [60].

Regulation of the Cold Response by miRNAs in Pea
Eleven miRNA families were represented in the 39 miRNA/target DEG pairs identified in this study, including miR156, miR159, miR166, miR167_1, miR396, miR397, miR482, miR2111, miR5770, miR5232 and a miRNA sequence which was not annotated. Their corresponding putative target genes provide a first insight into the main miRNA regulatory pathways involved in the pea response to cold.

Antioxidative Processes Are Predominantly Highlighted in Both Genotypes
Eight of the 22 miRNA targets differentially expressed in this study are annotated as genes potentially implied in antioxidative processes under abiotic stresses in plants. In our study, they appeared to be either commonly regulated in both genotypes or regulated in the sensitive genotype Te.
A first group of 3 genes, namely Psat4g012280, Psat6g124520 and Psat5g056680, respectively, annotated as UDP-glucoronosyl and UDP-glucosyl transferase, POT family and WD40/YVTN repeat-like-containing domain, are known to participate in the synthesis, regulation or transport of flavonoids. Flavonoids, including anthocyanins, play an important role in scavenging reactive oxygen species (ROS) implicated in many biotic and abiotic stresses. In At, for example, the expression of a citrus UDP glucosyl transferase gene enhanced the accumulation of anthocyanins and flavonoids, which was related to a better antioxidant potential under high light stress [61]. In tea leaves (Camellia sinensis L.), the CsUGT75C1 gene, encoding UDP-glucose:anthocyanin 5-O-glucosyl transferase, was upregulated in response to cold, along with an increase of the anthocyanin content [62]. In rice, a POT protein was mentioned as a component of the regulation of flavonoid metabolism during pollen intine formation [63]. In pea, many genetic loci are known to affect the synthesis of anthocyanins and were described for their role in the pigmentation of flowers and other organs [64,65]. Only three of these loci have been identified until now [66,67] and their sequence allowed us to note that none of them corresponded to the target genes mentioned here.
The gene Psat5g196400 is annotated as squalene/phytoene synthase, an enzyme which regulates the first step of the carotenoid biosynthetic pathway in plants. Similar to flavonoids, the carotenoid pigments are known to protect photosynthesis from ROS under various environmental stresses [68]. Phytoene synthase expression was shown to be regulated under abiotic stresses such as heat and drought stress in wheat [69] and heat stress in maize [70]. In these species, the allelic variation of PSY genes has been explored to breed for enhanced carotenoids levels under stresses.
In addition to the regulation of pigments, four other genes are potentially implied in ROS scavenging. They are Psat7g006120, Psat7g216840, Psat6g046440 and Psat1g035040, respectively, annotated as glutathione peroxydase, ankyrin repeats and ring finger domain for the last two. Glutathione peroxidase participates in ROS enzymatic scavenging systems in plants by reducing H 2 O 2 and organic and lipid hydroperoxides [71]. In winter rapeseed, for instance, the coordinated level of superoxyde dismutase, ascorbate peroxydase and glutathione peroxydase activities in cold acclimated young leaves of the cultivar Hansen was shown to play an essential role in order to avoid oxidative damage during cold exposure [72]. Plant ankyrin repeat (ANK) proteins are involved in various important biological processes, including response to biotic and abiotic stresses [73]. For example, the GmANK114 gene was overexpressed in soybean hairy roots in response to drought or salt stress, and lower ROS contents were also observed comparatively to control plants. The ring finger domain proteins are a large family, among which the ring-type E3-ubiquitin ligases have received a particular attention due to their potential impact in plant development, hormonal control and defense against biotic and abiotic stresses [74]. Regarding the response to cold stress, Fang et al. [75] functionally characterized the rice OsSRFP1 gene, coding the stress-related ring finger protein 1. They evidenced the increased cold tolerance of OsSRFP1 knock-down rice plants, associated with higher amounts of free proline and activities of antioxidant enzymes. miR156, miR166 and miR159, which are miRNAs conserved among plants, are the most strongly represented miRNAs in the present study. Five among the eight genes discussed above were regulated by members of the miR156 or miR166 families, as were almost all the genes commonly regulated by Ch and Te. The differential expression of miR156 and miR166 under low temperature conditions has been shown in various plant species [18] in organs of plantlets which are comparable to those studied here for pea. In Brassica rapa, for instance, miR156 and miR166 appeared among the 11 miRNA families whose members showed differential expression patterns under freezing stress (−4 • C), compared to control [76]. In Triticum aestivum, miR156, miR166 and miR159 were differentially expressed between cold (0 • C) and control [77]. In Medicago sativa, which such as Ps belongs to the Fabaceae family, the same three miRNAs were represented among the 27 miRNA families differentially expressed under cold (4 • C) and freezing (−8 • C) conditions [78].
Among cold responsive miRNAs, miR156 is also known to play an important part in plant morphology and in transition from the vegetative to the reproductive phase. Moreover, a coordinating role of miR156 in the relationship between plant development and abiotic stress tolerance has been demonstrated in At under salt and drought signals [79]. This previous study evidenced that the expression level of miR156 regulated, in a coordinated way, downstream genes implied in floral transition and in anthocyanin biosynthesis. It is interesting to notice that, in the present study, miR156 simultaneously regulated Psat0s740g0280, annotated as Fip1 motif, and Psat4g012280 discussed above, annotated as UDP-glucoronosyl and UDP-glucosyl transferase. In At, FIP1 is thought to mediate alternative polyadenylation and thus to contribute to regulation of the transcriptome [80,81]. The fip1-2 mutant of At exhibited multiple phenotypes affecting plant height, leaf size, flowering time and response to Cadmium and salt stress [81]. Additionally, as stated above, UDP glucosyl transferase is implied in the production of anthocyanins under cold conditions [62], which could be related to ROS scavenging. Thus, even if the miRNA-target modules vary among species, it would be of great interest to further study in pea the potential co-regulation of plant morphology, phenology and response to cold by miR156.
Lastly, in addition to their response to cold, there is also evidence of the induction of miR156, miR159 and miR166 by multiple stresses, including biotic stresses. For instance, in wheat, these three miRNAs were induced in response to cold [77], but also in response to powdery mildew [82]. In pea, the co-regulation by miRNAs of the cold and biotic stress responses merits additional exploration, as the question of a common genetic determinism of frost tolerance and Ascochyta resistance also arises from the colocalization of main QTLs (quantitative trait loci) for both traits [83].

Laccases Are Specifically Involved in Cold Adaptation of the Tolerant Genotype Ch
Six pairs of miRNAs/mRNAs with discordant expression patterns were identified at T1 for the frost tolerant pea accession Ch. Half of the miRNA targets were annotated as multicopper oxidases, which are an enzyme superfamily widely distributed in avian, bacteria, fungi, insects, mammalians and plants [84]. Using InterProScan [85], Psat2g059600, Psat3g137280 and Psat5g161560 were identified more precisely as laccases. In the present study, they were regulated by members of miR167, miR397 and miR5770 families. The miR397-laccase module has been particularly reported as an actor of plant response to cold in various plant species. As laccases are involved in lignin polymerization, one hypothesis is that miR397 directly regulates cold tolerance via reduction of lignin content in cell wall, thereby increasing cell wall permeability and elasticity, which can be a mechanical advantage under cold and freezing stress [18]. A second hypothesis, based on an indirect role of miR397 and miR397 targets in the CBF (C-repeat Binding Factor) regulon, an important regulatory pathway of the cold response, was also evoked by Dong and Pei [19]. These authors showed that transgenic plants overexpressing miR397a exhibited an improved tolerance to chilling (4 • C) and freezing (−1 • C) stresses and that overexpression of miR397a was associated with enhanced transcript levels of cold-regulated CBF genes and downstream COR (Cold Responsive) genes. They thus hypothesized that the transcription of CBF and COR genes was regulated by the miR397 target genes (3 laccases and a casein kinase β subunit 3 in this study). This complicated regulation would require further study in pea where CBF genes have been identified among candidate genes underlying a major QTL for frost tolerance located on linkage group VI [8]. Interestingly, the favorable allele at this QTL is contributed by the frost tolerant genotype Ch, which also showed a specific regulation of laccases in response to cold in this study. Lastly, it is interesting to consider that miR397 has also been shown to regulate the response to biotic stress, as in Malus hupehensis, a Malus rootstock species, where transient overexpression of miR397b in leaves triggered an increased sensitivity to Botryosphaeria dothidea, owing to the reduced Laccase7 and lignin contents [86]. As such, the miR397-laccase module should be studied in pea as a potential determinant of co-regulation of cold and biotic stress responses.

Conclusions
This study contributes new insight on the mechanisms implemented in pea during cold stress response by elucidating miRNA-mediated regulations. Our functional investigation of miRNA-target modules particularly highlighted antioxidative mechanisms in both a tolerant and a sensitive pea genotype. The tolerant line moreover showed specific expression of the miR397-laccase module, which could have either a direct role in cell wall remodeling under cold exposure or an indirect role in cold acclimation via the CBF regulon. Finally, the expression of miRNA families which are conserved among plant species, such as miR156, miR159, miR166 and miR397, suggests that these miRNAs could participate in the co-regulation of the responses to cold and biotic stresses in pea, as well as in the regulation of plant morphology and phenology under low temperature.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13071119/s1, Additional File S1 contains a schematic representation of the study design. Additional File S2 comprises Table S1 (miRNAs identified in this study), Table S2 (Ch specific differentially expressed miRNAs at T1), Table S3 (Ch specific differentially expressed miRNAs at T2), Table S4 (Te specific differentially expressed miRNAs at T1), Table S5 (Te specific differentially expressed miRNAs at T2), Table S6 (Common miRNAs differentially expressed at T1), Table S7 (Common miRNAs differentially expressed at T2), Table S8 (Identification of miRNA targets), Table S9 (Ch specific differentially expressed genes at T1), Table S10 (Ch specific differentially expressed genes at T2), Table S11 (Te specific differentially expressed genes at T1), Table S12 (Te specific differentially expressed genes at T2), Table S13 (Common genes differentially expressed at T1), Table S14 (Common genes differentially expressed at T1), Table S15 (GO terms enrichment for the Ch specific differentially expressed genes), Table S16 (GO terms enrichment for the Te specific differentially expressed genes) and Table S17 (GO terms enrichment for the common genes differentially expressed). Additional File S3 represents the expression of conserved and Fabaceae-specific miRNAs.