Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 870446 10.3389/fgene.2022.870446 Genetics Original Research Small RNA Sequencing Revealed that miR4415, a Legume-Specific miRNA, was Involved in the Cold Acclimation of Ammopiptanthus nanus by Targeting an L-Ascorbate Oxidase Gene and Regulating the Redox State of Apoplast Zhu et al. miRNAs in Ammopiptanthus nanus Zhu Ming 1 2 3 Wang Xue 4 Zhou Yanqiu 1 2 3 Tan Jinhua 1 2 3 Zhou Yijun 1 2 3 * Gao Fei 1 2 3 * 1 Key Laboratory of Ecology and Environment in Minority Areas, National Ethnic Affairs Commission, Minzu University of China, Beijing, China 2 Key Laboratory of Mass Spectrometry Imaging and Metabolomics, National Ethnic Affairs Commission, Minzu University of China, Beijing, China 3 College of Life and Environmental Sciences, Minzu University of China, Beijing, China 4 Beijing Center for Disease Prevention and Control, Beijing, China

Edited by: Yuepeng Song, Beijing Forestry University, China

Reviewed by: Jiaxing Tian, Beijing Academy of Agricultural and Forestry Sciences, China

Cunfu Lu, Beijing Forestry University, China

 *Correspondence: Yijun Zhou, zhouyijun@muc.edu.cn; Fei Gao, gaofei@muc.edu.cn

This article was submitted to RNA, a section of the journal Frontiers in Genetics

 04 04 2022 2022 13 870446 06 02 2022 07 03 2022 Copyright © 2022 Zhu, Wang, Zhou, Tan, Zhou and Gao. 2022 Zhu, Wang, Zhou, Tan, Zhou and Gao

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

MicroRNAs (miRNAs) are small endogenous single-stranded RNAs that regulate plant growth, development, and environmental stress response posttranscriptionally. Ammopiptanthus nanus, a rare evergreen broad-leaved shrub in the temperate area of Central Asia, can tolerate freezing stress as low as –30 degrees centigrade in winter, and miRNA might be involved in the cold acclimation which enables A. nanus to obtain tolerance to freezing stress. Systematic identification and functional analysis of the miRNAs involved in the cold acclimation in A. nanus may promote understanding of the miRNA-mediated gene regulation network underlying cold acclimation. Here, based on small RNA and degradome sequencing, 256 miRNAs and 1,808 miRNA-target pairs were identified in A. nanus. A total of 39 cold-responsive miRNAs were identified, of which 29 were upregulated and ten were downregulated. These cold-responsive miRNAs may participate in the cold acclimation by regulating redox homeostasis (miR398, miR4415, and miR408), calcium signaling (miR5225 and miR5211), growth and development (miR159 and miR390), and small RNA–mediated gene silencing (miR168 and miR1515). We found that miR4415, a legume-specific miRNA, is involved in the cold acclimation of A. nanus by targeting an L-ascorbate oxidase gene and then regulating the redox state of the apoplast. Our study provides important data for understanding the regulatory role of miRNA in the cold acclimation of A. nanus.

 miRNA miR4415 cold acclimation Ammopiptanthus nanus L-ascorbate oxidase 31770363 31670335 National Natural Science Foundation of China10.13039/501100001809 Introduction

Low temperature is a common unfavorable factor that negatively affects plant growth, development, and productivity and limits the geographical distribution of plant species (Khan et al., 2017; Ding et al., 2020). According to the temperature range, low-temperature stress can be divided into cold stress (0–10°C, chilling temperature) and freezing stress (<0°C, freezing temperature). Cold stress can directly affect the activity of various biological macromolecules, including various enzymes and bio-membranes, resulting in metabolic disorders, photosynthesis inhibition, and increased levels of reactive oxygen species (ROS) (Liu et al., 2018). Under freezing stress, ice crystals are first formed in the intercellular space or apoplast of plants. The ice crystals gradually increase, leading to severe cell dehydration and plasma membrane damage and rupture, further causing intracellular protein denaturation, precipitation of various biological molecules, disorder of cell metabolism, large accumulation of ROS, and even cell death (Bredow and Walker, 2017). Plants in temperate regions can gradually develop tolerance to freezing stress by exposure to low nonfreezing temperatures in autumn, a process known as cold acclimation. During the long evolutionary process, plants have evolved fine-tuned molecular regulatory mechanisms at transcriptional, posttranscriptional, translational, and posttranslational levels to cope with low-temperature conditions (Krutzfeldt and Stoffel, 2006; Liu et al., 2019). Recent studies have shown that in addition to signal transduction components such as mitogen-activated protein kinase (MAPK) pathway and transcription factors such as dehydration responsive element-binding (DREB) that play key roles in stress signal transduction and gene expression regulation and noncoding RNAs, such as microRNA (miRNA), also play essential regulatory roles in cold acclimation (Yu et al., 2019).

miRNAs are a class of endogenous single-stranded non-coding small RNAs (sRNAs) of about 18–24 nt. miRNAs play vital regulatory roles in plant growth and development and respond to environmental factors by negatively regulating protein-coding genes at the posttranscriptional level, thus regulating various biological processes (Song et al., 2019). It has been shown that miRNAs are involved throughout the development of reproductive organs, leaves, and stems and respond to high salinity, drought, and low-temperature stress (Ma et al., 2019). The application of high-throughput sequencing techniques to investigate the effects of cold stress on the expression of the conserved and nonconserved miRNAs in plants has been reported in Arabidopsis thaliana (Sunkar and Zhu, 2004), Oryza sativa (Lv et al., 2010), Zea mays (Yang et al., 2011), Glycine max (Xu et al., 2016), and Triticum aestivum (Wang et al., 2014). Several miRNAs were found to be involved in plant response to cold stress, including miR156, miR159, miR167, miR168, miR169, miR172, miR319, miR393, miR396, miR398, and miR408 (Megha et al., 2018), and further functional studies revealed the regulatory roles of cold-responsive miRNAs in cold stress response in plants. For example, it was demonstrated that miR156 enhanced cold stress tolerance in rice by targeting the expression of transcription factor genes such as SQUAMOSA promoter binding protein-like (SPL) (Zhou and Tang, 2019). Furthermore, miR528 decreased the expression of transcription factor MYB30 by regulating an F-box domain–containing protein gene (Os06g06050), leading to enhanced cold stress tolerance in rice (Tang and Thompson, 2019).

Ammopiptanthus nanus is a desert shrub belonging to Ammopiptanthus, Leguminosae. This plant species is mainly distributed in the junction of the western Tianshan Mountains and Kunlun Mountains in Xinjiang, Wuqia county of Xinjiang Uygur Autonomous Region, and Uzbekistan. A. nanus originates in the Paleotropical region and is a tertiary relic that survived through Mediterranean retreat and climate aridity. A. nanus is of considerable scientific value in investigating the changes of paleogeography, paleoclimate, and paleoflora in Central Asia. In addition, A. nanus plays a vital role in the maintenance of the fragile ecosystems in arid areas in Central Asia. It can also be used as an appropriate plant species for windbreak, sand fixation, and greening in arid areas. A. nanus has evolved high tolerance to low temperature, drought, and other adverse conditions in temperate deserts. It is noteworthy that A. nanus is a rare, evergreen, broad-leaved woody plant in the desert area of Central Asia. A. nanus can tolerate low-temperature environments as low as –30°C. Thus, A. nanus was proposed to be a crucial material for analyzing the mechanism of low-temperature tolerance and screening candidate genes used for genetic engineering for improving low-temperature tolerance of economic plants (Gao et al., 2016; Cao et al., 2019; Cao et al., 2020; Zhang et al., 2021). However, system identification of miRNAs in A. nanus has not been reported, and how miRNAs respond to cold acclimation in A. nanus is not clear till now.

In the current study, we conducted genome-wide identification of the conserved and species-specific miRNAs in A. nanus based on deep sequencing and performed the expression profiling of miRNAs and their targets under cold acclimation. Several conserved and species-specific miRNAs and their targets were identified, and some miRNAs were demonstrated to play roles in cold acclimation of A. nanus via negatively regulating their targets. We further analyzed the biological functions of a lineage-specific miRNA, miR4415, in cold acclimation. Our results provided important data for understanding the miRNA-mediated gene regulation network underlying the cold acclimation in A. nanus.

 Materials and Methods Plant Materials and Cold Acclimation Treatment

Seeds of A. nanus were collected from Wuqia county (75°1′E, 39°43′N), Xinjiang Autonomous Region, China. After surface sterilization with 75% (v/v) ethanol and soaking in water for 48 h at room temperature, the seeds were planted in commercial pots containing a 3:1 (v/v) mixture of perlite and vermiculite. Seedlings were grown in a plant incubator with a photoperiod of 16 h, at 22–28°C and relative humidity of 35%. The seedlings were irrigated every 4 days with a solution of half-strength Hoagland. Eight weeks after germination, seedlings with similar heights were randomly divided into two groups. The first group of seedlings were kept in the plant incubator at 22–28°C for 7 days, and half of the plants were used as the control group (CK), and the other half were exposed to freezing stress treatment (–30°C) for 1 day and used as the freezing-stressed control group (SCK). The other group of seedlings were exposed to cold treatment (5°C/4°C, light/dark) for 7 days for cold acclimation, and half of the plants were used as the cold-acclimated group (CA), and the other half were exposed to freezing stress treatment (–30°C) for 1 day and used as the freezing-stressed cold-acclimated group (SCA). Leaf samples were harvested from each group and used for biochemical parameter measurement. Leaf samples from CK and CA were used for small RNA sequencing, and leaf samples from CK were used for degradome sequencing.

 Biochemical Parameter Measurement

Malondialdehyde (MDA) and relative electrolyte leakage (REL) were measured using the methods described by Yasar and Yang, respectively (Yang et al., 1996; Yasar et al., 2008). The activities of malate dehydrogenase (MDH) and ascorbate oxidase (AO) were measured using kits manufactured by Solarbio Science & Technology Co., Ltd. (Beijing, China). The content of ascorbate (ASA) and dehydroascorbate (DHA) was measured using kits manufactured by Solarbio Science & Technology Co., Ltd. (Beijing, China). Five independent plant samples were used (n = 5) for each biochemical parameter measurement.

 Small RNA and Degradome Sequencing

RNA samples were prepared using the TRIzol reagent. The concentration and purity of total RNA samples were assayed using Bioanalyzer 2100 and the RNA 6000 Nano LabChip Kit (Agilent, CA, United States) with a RIN value of >7.0. For small-RNA sequencing, approximately 1 μg of total RNA was used to construct a small-RNA library according to the protocol of the TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, United States). Single-end sequencing (36 bp) was conducted on an Illumina Hiseq2500 at the LC Sciences (Hangzhou, China) by following the manufacturer’s instructions. Three independent biological replicates were set for each group, and a total of six high-throughput small RNA sequencing data were surveyed. The raw data of small RNA sequencing have been submitted to the NCBI SRA database with accession numbers SRR10906487, SRR10906486, SRR10906485, SRR10906484, SRR10906483, and SRR10906482.

For degradome sequencing, 150 ng of the poly (A) RNA sample was used as input RNA and annealed with biotinylated random primers. After RNA fragments were captured using biotinylated random primers, a 5ʹ adapter was ligated to the RNAs with monophosphate at the 5ʹ end. Next, reverse transcription and PCR were performed. Finally, the resulting library was single-end sequenced using an Illumina Hiseq2500 platform at the LC Sciences (Hangzhou, China) by following the vendor’s recommended protocol. The raw data of degradome sequencing have been submitted to the NCBI SRA database under the accession number SRR17711209.

 Identification of the Conserved and Species-Specific miRNAs

miRNAs were predicted using the ACGT101-miRprogram (version 4.2, LC Sciences), and the predicted miRNAs were then checked manually to remove false miRNAs. The genome sequence of A. nanus was downloaded from GigaBase (Gao et al., 2018). The conserved miRNAs and species-specific miRNAs were determined by aligning to the mature sequences of miRNAs from all green plants in the miRBase database (http://www.mirbase.org) (Kozomara et al., 2019). Species-specific miRNAs are miRNAs that are not aligned to the miRNAs deposited in miRBase (mismatches >4 nt) and meet the criteria of high-confidence miRNAs according to the latest miRNA identification protocol (Kozomara and Griffiths-Jones, 2014).

 miRNA Target Prediction and Degradome Analysis

The targets of miRNAs were predicted by using psRNATarget with the settings (maximum expectation: 4.0; length for complementary scoring: 20 bp; target accessibility–allowed maximum energy to unpair the target site: 25.0; flanking length around target site for target accessibility analysis: 17 bp in upstream and 13 bp in downstream; range of central mismatch leading to translational inhibition: 9–11 nt) (Dai et al., 2018). CleaveLand 3.0 was used for analyzing sequencing data with default parameters (Addo-Quaye et al., 2009). CPC2.0 was used to evaluate the coding potential of miRNA targets (Kang et al., 2017).

 Differential Expression Analysis of miRNAs

The read counts from CK and CA were normalized to transcripts per million reads (TPM) (Dewey, 2010), and the differentially expressed miRNAs were determined by comparing the conserved or species-specific miRNA expression between two groups. Genes with a fold change ≥1.5 or ≤0.67 and a p-value <0.05 were identified as differentially expressed miRNAs. The heatmap was generated by using TBtools (Chen et al., 2020).

 qRT-PCR Analysis of miRNAs and Their Targets

RNA samples were prepared from the leaves using TRIzol reagent. qRT-PCR analysis of miRNAs and their targets were performed using the methods described previously (Abla et al., 2019). All primers were designed based on the sequences of selected miRNAs and their targets (Supplementary Table S1). The internal control genes are U6 snRNA (for miRNA) and 18S rRNA (for target gene). Three technical replicates were performed for each analyzed gene, and the expression levels of miRNAs and their target genes were calculated using the 2−ΔΔCt method (Varkonyi-Gasic et al., 2007).

 Promoter Analysis

Plant promoter prediction software PlantCARE (http://intra.psb.ugent.be: 8080/PlantCARE) was used to predict the cis-acting elements involved in stress response or hormone response, and 1000 bp upstream sequences of the stem-loop precursors of each miRNA were used for promoter prediction.

 Identification, Multiple Sequence Alignment, and Phylogenetic Analysis of miR4415

To identify miR4415 and their precursors, mature sequences of know miR4415 were downloaded from miRBase and used to search for homologous sequences using the BLAST program, allowing four mismatches. Then, 400-bp flanking sequence was retrieved from reference sequences and folded using the online tool mfold to evaluate their secondary structure. Only those genomic loci with a stem-loop structure with an MFE value less than 30 kcal/mol and having the typical secondary structure for miRNA precursors were retained for subsequent analysis. All plant genomes in Phytozome v12 (https://genome.jgi.doe.gov/portal/PhytozomeV12) were downloaded and used for identifying miR4415. Multiple sequence alignment was performed using Clustal and visualized with Weblogo. The phylogenetic tree of the precursor sequences of miR4415 was built by MEGA X using the neighbor-joining (NJ) method.

 Dual-Luciferase Reporter Assay

Dual-luciferase reporter assay was performed according to the method for studying miRNA function (Clement et al., 2015). The precursor sequence of A. nanus miR4415 was ligated to the pGreen II 62 SK vector, and the coding sequence of L-AO was ligated to the pGreen II 0800-LUC vector. All primers were designed based on the sequences of Pre-miR4415, L-AO, pGreen II 62 SK, and pGreen II 0800-LUC vector (Supplementary Table S1). Extraction of A. thaliana protoplasts was performed as described by Yoo et al. (2007). The luciferase signals were detected using a Bio-Lite TM Luciferase Assay system kit manufactured from Vazyme Biotech Co., Ltd. (Nanjing, China). Five independent repeats were used for this assay.

 Subcellular Localization Analysis

Subcellular localization analysis of L-AO protein was performed as described in the literature (Yin et al., 2021). The coding sequence of the L-AO gene was ligated to pCAMBIA1303, and L-AO was transiently expressed in tobacco leaves for 48–72 h and imaged by using a confocal laser-scanning microscope (Leica TCS Sp8).

 Extraction of the Leaf Apoplast Fluid of <italic>Ammopiptanthus nanus</italic>

The leaf apoplast fluid of A. nanus was obtained using the infiltration–centrifugation methods (O'Leary et al., 2014). Water was used as the buffer, a 60-ml syringe was used for leaf infiltration, and the resulting leaf apoplast fluid was recovered by centrifugation for 10 min at 3,000 g.

 Results Cold Acclimation Enhances the Tolerance to Freezing Stress in <italic>Ammopiptanthus nanus</italic>

To evaluate the effect of cold acclimation on the tolerance to freezing stress in A. nanus, cold-acclimated and untreated A. nanus seedlings were subjected to freezing stress treatment, and MDA and REL, two biochemical indicators reflecting the degree of cell damage, were measured. The seedlings in SCK exhibited significant wilting, while the seedlings in SCA showed only slight wilting (Figure 1A). The MDA and REL values in SCA were significantly lower than those of SCK (Figures 1B,C), indicating that cold acclimation alleviated the damage to seedlings of A. nanus caused by the freezing stress treatment. In other words, high freezing tolerance can be obtained by cold acclimation in A. nanus seedlings.

Cold acclimation enhanced the freezing stress tolerance of A. nanus. (A) Phenotype of A. nanus seedling in CK, SCK, CA, and SCA. The change of MDA (B) and REL (C) in the leaves of A. nanus before (CK and CA) and after freezing treatment (SCK and SCA). Values are expressed as means ± SD. Duncan’s method was used for multiple comparison analysis of variance, n = 3.

 Cold Acclimation Resulted in a Significant Alteration in the Small RNA Population in <italic>Ammopiptanthus nanus</italic>

To evaluate the influence of cold acclimation on the sRNA population and identify miRNAs involved in cold acclimation in A. nanus, six sRNA libraries were constructed by using leaf samples isolated from CK and CA seedlings. CK and CA samples from A. nanus were sequenced by high-throughput sequencing technology. Approximately, 18 M raw sequence reads were obtained for each library. After removing adapters, low-quality sequences and sequences smaller than 18 nt, 12.28–16.10 M clean reads with sizes ranging from 18 to 25 nt were obtained for these sRNA libraries (Table 1).

Statistics of the clean reads generated by high-throughput sequencing.

Library Number of raw reads Number of high-quality reads Reads without 3′ adapter Reads without insert fragment 5′ adapter contaminants Reads with polyA Reads smaller than 18nt Number of clean reads
miCK_1 176,958,37 (100%) 173,638,56 (98.12%) 94,823 (0.54%) 158,912 (0.91%) 71,408 (0.41%) 728 (0.01%) 4,759,128 (27.40%) 12,278,857 (70.71%)
miCK_2 18,724,118 (100%) 18,316,541 (97.82%) 70,928 (0.38%) 147,893 (0.80%) 52,832 (0.28%) 1340 (0.01%) 3,316,221 (18.10%) 14,727,327 (80.40%)
miCK_3 19,377,347 (100%) 19,001,586 (98.06%) 61,387 (0.32%) 88,535 (0.46%) 59,142 (0.31%) 1097 (0.01%) 2,690,773 (14.16%) 16,100,652 (84.73%)
miCA_1 18,647,262 (100%) 18,277,927 (98.01%) 108,694 (0.59%) 128,250 (0.70%) 53,975 (0.29%) 918 (0.01%) 3,387,959 (18.53%) 14,598,131 (79.86%)
miCA_2 19,485,083 (100%) 19,126,024 (98.15%) 66,229 (0.34%) 109,935 (0.57%) 82,378 (0.43%) 576 (0.01%) 5,235,308 (27.37%) 13,631,598 (71.27%)
miCA_3 18,046,481 (100%) 17,705,703 (98.11%) 66,387 (0.37%) 130,810 (0.73%) 75,991 (0.42%) 504 (0.01%) 4,288,459 (24.22%) 13,143,552 (74.23%)

The size distribution pattern of the three sRNA biological replicates in each group (CK or CA) was similar (Figure 2). In group CK, the top four abundant small RNA size categories were 24, 21, 22, and 20 nt groups, and such a size distribution pattern was similar to that of the typical plant mature sRNAs reported in several plant species, for example, A. thaliana, Medicago truncatula (Szittya et al., 2008), O. sativa (Morin et al., 2008), and Citrus trifoliate (Song et al., 2010). While there was a significant difference in the sRNA size distribution pattern between the CK and CA groups, the top two abundant sRNA size categories were 21 nt and 24 nt. These data suggested that cold acclimation resulted in a significant alteration in the sRNA population in A. nanus leaves.

Length distribution of sRNA sequences in the six sRNA libraries prepared from leaves of the A. nanus seedlings cultured in untreated conditions (CK) (A) and under cold acclimation (CA) (B).

 Identification of the Conserved and the Species-Specific miRNAs in <italic>Ammopiptanthus nanus</italic>

To identify miRNAs, all sRNAs of the six libraries were mapped to A. nanus genome sequences, and the flanking sequences were folded to determine whether a qualified hairpin could be formed. Consequently, 188 conserved mature miRNA sequences and 107 corresponding stem-loop precursors were identified (Supplementary Table S2). All the stem-loop precursors showed similarity to the miRNA precursors of green plants in the miRBase (Supplementary Table S2), indicating that these conserved miRNAs are generated from conserved stem-loop precursors. Most of the conserved miRNAs are supported by the simultaneous presence of the hairpin 5p and 3p arms, but no miR* sequences have been found for 80 stem-loop precursors. Among the 80 precursors, 35 with only 5p arm miRNAs and 45 with 3p arm miRNAs were found.

By aligning the stem-loop sequences obtained above with the miRBase database, 41 conserved miRNA families were identified (Figure 3A). MIR159 and MIR156, with up to 11 precursor loci, are the largest miRNA families, followed by MIR166 (10), MIR395 (9), and MIR171-1 (9). Twenty miRNA families, such as MIR162, MIR394, MIR477, MIR862-2, MIR1507, MIR1515, and MIR5559, had only one predicted precursor.

Statistics on the number of precursors of the miRNA families (A) and the secondary structure of the miRNA precursors of ana-miRN4 (B), ana-miRN9 (C), and ana-miRN33 (D). Mature sequences of miRNAs are shown in uppercase, with red representing 5p and blue representing 3p. The structure diagrams were generated by using the UNAFold Web Server.

In the current study, all predicted miRNAs that have no homologs (with less than four mismatches) in miRBase were classified as species-specific miRNAs. According to the criteria for distinguishing the high- and low-confidence miRNA annotations (Kozomara and Griffiths-Jones, 2014), a total of 68 nonconserved miRNAs and 39 corresponding stem-loop precursors were identified with high confidence (Figures 3B–D; Table 2; Supplementary Tables S3).

Species-specific miRNAs predicted from A. nanus.

ID miRNA name Sequence Length (nt) Total reads ΔG (kcal/mol)
N1 ana-miRN1-5p AAG​CTC​TTG​CTA​GGT​TGA​TTG 21 12 −70.20
ana-miRN1-3p ATT​GCC​TAC​AGT​TAG​ATC​CTT​A 22 281
N2 ana-miRN2-5p AGG​ATA​TTG​CTG​GGT​TGA​TTG 21 57 −76.80
ana-miRN2-3p ATT​ACC​TAT​AGT​TTT​ATC​CT 20 3
N3 ana-miRN3-5p TAC​TAA​TTT​GGT​CCT​CAA​GAA 21 11 −66.80
ana-miRN3-3p TCT​TGA​GGA​TCA​AAT​TTG​TAT​T 22 34
N4 ana-miRN4-5p TGC​CAA​CTT​GAA​GAA​CAA​ACA 21 150 −56.60
ana-miRN4-3p TTT​GTT​CTT​CAA​GTT​GGC​ATC 21 40,865
N5 ana-miRN5-5p TGT​TTG​GGA​GAA​GTG​GGA​AAG 21 243 −55.80
ana-miRN5-3p TTT​TCC​GCT​TCT​TTC​AAA​CAG 21 1791
N6 ana-miRN6-5p TCA​TTT​AAT​CCC​TAT​ACG​TAG 21 25 −38.80
ana-miRN6-3p TAA​CGG​GAG​GGA​CCA​AAT​GGG 21 101
N7 ana-miRN7-5p TAT​GGT​GAT​ACA​TAT​CTG​ATG 21 709 −36.10
ana-miRN7-3p TTA​GAT​ATG​TAT​CTC​CTT​ATT 21 772
N8 ana-miRN8-5p TCC​GTG​GCC​AAG​ATG​ACG​AAG​AA 23 1437 −40.70
ana-miRN8-3p CTC​TTC​ATC​TTC​CTC​GGA​TTC 21 43
N9 ana-miRN9-5p TTT​AGG​ATA​GGC​TTT​GAT​ACC 21 630 −54.30
ana-miRN9-3p TAT​CAA​AGC​CTA​TCC​GAG​ATC 21 16
N10 ana-miRN10-5p ATT​CTT​CAT​TTA​GTC​TCT​ATA 21 334 −40.70
ana-miRN10-3p TAG​GGA​CTA​AAT​GAA​GAA​TCT 21 10
N11 ana-miRN11-5p CAT​TAA​CGA​AAT​TCC​AAA​AAC 21 69 −93.40
ana-miRN11-3p TTT​TGG​AAT​TTC​GTT​AAT​ATC 21 393
N12 ana-miRN12-5p TCC​GCC​GAG​ACT​CGT​CTG​CA 20 172 −81.90
ana-miRN12-3p CAG​ACG​GGT​CAT​GGC​AGA​CGA 21 21
N13 ana-miRN13-5p TTT​CCT​GTC​ACA​TTC​TCT​GTC 21 49 −93.70
ana-miRN13-3p CAG​AGA​ATG​TGA​CAG​GAA​AGT​G 22 12
N14 ana-miRN14-5p TTG​TAC​GAT​TTT​GGT​CCC​TCA 21 12 −85.00
ana-miRN14-3p TGA​GGG​ACC​AAA​ATT​ACA​CAA​T 22 43
N15 ana-miRN15-5p TGG​GAT​TAC​AGG​GGT​TCT​CTT 21 10 −82.30
ana-miRN15-3p AGA​GAA​CCC​CTG​TAA​TCC​CAG 21 36
N16 ana-miRN16-5p TTA​CGG​GAC​ATT​CTT​ATG​TGG​C 22 14 −43.40
ana-miRN16-3p CAC​ATA​GGA​ATG​ACA​CGT​AAG​CT 23 20
N17 ana-miRN17-5p TTT​TGT​TTT​TGG​TCC​CTG​TCA 21 18 −54.00
ana-miRN17-3p CAG​GGA​CCA​AAA​ACA​AAA​TTT​T 22 10
N18 ana-miRN18-5p CCC​ACC​CAC​ACA​GTT​GAC​CTA​A 22 1057 −37.10
ana-miRN18-3p AGG​TCA​AAT​TTG​TGG​GGG​GTG 21 67
N19 ana-miRN19-5p GAG​GGA​CTA​AAT​GAA​GAA​TTT 21 105 −123.00
ana-miRN19-3p ATT​CTT​CAT​TTA​GTC​TCT​ATA 21 347
N20 ana-miRN20-5p TAA​GGA​CTA​AAT​GAA​GAA​TTT 21 32 −63.90
ana-miRN20-3p ATT​CTT​CAT​TTA​GTC​TCT​ATA 21 347
N21 ana-miRN21-5p CAT​TAA​CGA​AAT​TCC​AAA​AAC 21 67 −78.70
ana-miRN21-3p TTT​TGG​AAT​TTC​GTT​AAT​ATC 21 365
N22 ana-miRN22-5p CGT​CCT​GGT​GAA​ACG​CGC​CAC​T 22 307 −70.10
ana-miRN22-3p TGG​CGC​GTT​CCA​CCA​GGA​TGC 21 11
N23 ana-miRN23-5p TCC​ACT​GAA​CGT​AAT​TAA​CCA 21 227 −57.70
ana-miRN23-3p ATT​AAT​TAC​GTT​CAG​TGG​ATG 21 13
N24 ana-miRN24-5p TTG​ATC​AAC​TCC​TCC​ACC​GTG​A 22 25 −34.90
ana-miRN24-3p CAC​GGT​GAG​AGT​TGG​ACT​TTG​C 22 200
N25 ana-miRN25-5p TCC​ACT​GAA​CGT​AAT​TAA​CCA 21 210 −90.90
ana-miRN25-3p ATT​AAT​TAC​GTT​CAG​TGG​ATG 21 15
N26 ana-miRN26-5p TAC​CAA​TTT​GAT​CCT​CAA​GAA​T 22 10 −68.90
ana-miRN26-3p TTT​GAG​GAT​CAA​ATT​GGT​ATT 21 194
N27 ana-miRN27-5p TGG​CGC​GAC​GCA​GTG​GAA​GGC 21 72 −114.50
ana-miRN27-3p TCC​ACT​GCG​TCG​CGC​CAC​GTG 21 67
N28 ana-miRN28-5p TGG​CGC​GAC​GCA​GTG​GAA​GGC 21 43 −83.10
ana-miRN28-3p TCC​ACT​GCG​TCG​CGC​CAC​GTG 21 67
N29 ana-miRN29-5p TCA​CTC​ATC​GTT​GGA​TCA​ATC 21 75 −32.50
ana-miRN29-3p TTG​ATC​TAA​CGG​TGC​GTG​AAT 21 18
N30 ana-miRN30-5p TAG​CAC​ATC​ATG​TCC​CAA​TCA 21 12 −47.80
ana-miRN30-3p ATT​GGG​ATA​TAA​TGT​GAT​ACA​T 22 55
N31 ana-miRN31-5p CAT​TAA​CCA​CGA​TTT​TGA​ACG 21 37 −91.40
ana-miRN31-3p TTC​AAA​ATC​GTG​GTT​AAT​GAA 21 27
N32 ana-miRN32-5p TCC​GGA​TCC​TCT​AAC​TTT​AGG 21 49 −55.70
ana-miRN32-3p TAA​AGT​TAG​AGG​ATC​CAA​ATC 21 20
N33 ana-miRN33-5p CCA​AAT​CTC​AAT​CGT​TGG​ATT 21 39 −58.90
ana-miRN33-3p TCC​AAC​GGT​TGA​GAT​TTG​AGT 21 34
N34 ana-miRN34-5p CAT​TAA​CCA​CGA​TTT​TGA​ACG 21 37 −72.40
ana-miRN34-3p TTC​AAA​ATC​GTG​GTT​AAT​GAA 21 21
N35 ana-miRN35-5p TAT​CAA​TTT​GGT​CCT​CAA​GAA 21 32 −56.07
ana-miRN35-3p TCT​TGA​GGA​TCA​AAT​TTG​TAT​T 22 34
N36 ana-miRN36-5p TTT​GGA​TCC​TCT​AAC​TTT​AGG 21 15 −56.58
ana-miRN36-3p TAA​AGT​TAG​AGG​ATC​CGG​ATC 21 26
N37 ana-miRN37-5p TTT​GGA​TCC​TCT​AAC​TTT​AGG 21 15 −50.50
ana-miRN37-3p TAA​AGT​TAG​AGG​ATC​CGG​ATC 21 26
N38 ana-miRN38-5p TAT​CAA​TTT​GGT​CCT​CAA​GAA 21 32 −50.10
ana-miRN38-3p TTG​AGG​ACC​AAA​TTA​ATA​TTT​T 22 13
N39 ana-miRN39-5p CGG​AAA​ATT​GTT​GCA​GTT​AAG​C 22 41 −83.10
ana-miRN39-3p TTA​ACT​GCA​ACA​ATT​TGT​CCA​T 22 2655
Target Identification of miRNAs from <italic>Ammopiptanthus nanus</italic> Using Bioinformatics Analysis and Degradome Sequencing

To dissect the function of the predicted miRNAs in A. nanus, psRNAtarget, a miRNA target predicting tool, was used to identify miRNA target genes, and a total of 15,005 miRNA target pairs were predicted. Among them, 159 conserved miRNAs formed 10,159 miRNA-target pairs (Supplementary Table S4), and 64 nonconserved miRNAs formed 4,846 miRNA-target pairs (Supplementary Table S5).

Through Gene Ontology (GO) functional classification analysis, the top 10 of the biological process (BP), cellular component (CC), and molecular function (MF) categories were obtained (Figure 4A). Among these groups, the term regulation of programmed cell death (GO: 0012501), cell death (GO: 0008219), CCAAT-binding factor complex (GO: 0016602), protein serine/threonine phosphatase complex (GO: 0008287), nucleoside triphosphatase activity (GO: 0017111), and pyrophosphatase activity (GO: 0016462) was the top two abundant terms in BP, CC, and MF categories. These results indicated miRNAs might be involved in pathways such as programmed cell death, chromatin assembly or disassembly, signal transduction, enzyme activity, transcription, and secondary metabolism in A. nanus by negatively regulating their target genes.

GO analysis and T-plots of miRNA targets in A. nanus. GO analysis of the miRNA targets (A). Red columns represent BP, green columns represent CC, and blue columns represent MF. T-plots of miR160-ARF (B), miR159-MYB (C), miR169-NF-YA (D), miR1507-NB-ARC (E), miR1509-comp35748_c0_seq1 (F), and miR1511-comp15186_c1_seq10 (G) pairs identified by degradome sequencing. Cleavage sites on target genes are indicated by red dots.

To experimentally identify the target genes of the predicted miRNAs in A. nanus in batch, degradome sequencing was conducted. High-throughput sequencing produced 11,048,980 reads representing the 5’ ends of uncapped, poly-adenylated RNAs. The total number of signatures matching the transcriptome sequences was 10,151,291 (91.88%), and the number of distinct sequences matching A. nanus transcriptome sequences was 2,440,623. CleaveLand pipeline was used to identify the cleaved miRNA targets in the A. nanus transcriptome. A total of 1,808 nonredundant putative splicing sites with p-values ≤1 were identified confidently (categories ≤ 4) (Supplementary Table S6). Among these predicted targets, 81 were classified in category 0, and they were the most credible targets. Of the 81 targets, 47 were known targets that were reported in other plant species previously, and 34 were noncanonical targets that were identified for the first time. Most of the 34 targets (21/34) were protein-coding genes, including comp13662_c0_seq4 and comp15069_c7_seq8, which encode a global transcription factor group E8 and a basic-leucine zipper (bZIP) transcription factor family protein, respectively. Of the other 13 noncanonical targets, ten were predicted to be noncoding transcripts by CPC2.0, a coding potential calculator (Kang et al., 2017). (Figures 4B–G; Table 3).

Noncanonical targets of the predicted miRNAs identified by degradome sequencing in A. nanus.

miRNA Non-canonical target Annotation/coding potential
ana-miR1507a-3p comp13515_c0_seq1 Pre-rRNA-processing protein TSR2, conserved region
ana-miR1507a-3p comp15069_c7_seq8 Basic-leucine zipper (bZIP) transcription factor family protein
ana-miR1507a-3p comp19088_c0_seq1 TIM23-2 | translocase inner membrane subunit 23–2
ana-miR1509a-5p comp15920_c0_seq10 hydrolases; protein serine/threonine phosphatases (PP)
ana-miR1509a-5p comp15920_c0_seq1 hydrolases; protein serine/threonine phosphatases
ana-miR1511a-3p comp15499_c1_seq10 Bromo-adjacent homology (BAH) domain-containing protein
ana-miR1511a-3p comp15499_c1_seq2 Bromo-adjacent homology (BAH) domain-containing protein
ana-miR1511a-3p comp15499_c1_seq9 Bromo-adjacent homology (BAH) domain-containing protein
ana-miR1511a-3p comp15745_c1_seq1 GAUT6 | galacturonosyltransferase 6
ana-miR1511a-3p comp19214_c0_seq1 ARM repeat superfamily protein
ana-miR1511a-3p comp6070_c1_seq1 SWIB complex BAF60b domain-containing protein
ana-miR156d-3p comp13662_c0_seq4 GTE8 | global transcription factor group E8
ana-miR156i-3p comp15773_c2_seq28 DNA glycosylase superfamily protein (DGS)
ana-miR156j-5p comp11270_c0_seq3 BPM2 | BTB-POZ and MATH domain 2
ana-miR156j-5p comp16323_c0_seq1 Glycine cleavage T-protein family
ana-miR159h-3p comp4558_c0_seq1 Mitochondrial transcription termination factor family protein
ana-miR159i-3p comp14691_c0_seq1 Ypt/Rab-GAP domain of gyp1p superfamily protein (Gyp1p)
ana-miR171-1b-5p comp15661_c0_seq18 Phox-associated domain, Phox-like, Sorting nexin, C- terminal (PX)
ana-miR171-1f-3p comp14326_c1_seq1 Nucleotide/sugar transporter family protein
ana-miR396c-5p comp11385_c0_seq1 LACS1 | AMP-dependent synthetase and ligase family protein
ana-miR396c-5p comp19381_c0_seq1 BSD domain-containing protein
ana-miR1509a-5p comp15155_c0_seq14 Noncoding
ana-miR1509a-5p comp35748_c0_seq1 Noncoding
ana-miR1509b-5p comp14487_c0_seq2 Noncoding
ana-miR1509b-5p comp15155_c0_seq6 Noncoding
ana-miR1511a-3p comp15186_c1_seq10 coding
ana-miR1511a-3p comp15876_c0_seq40 Noncoding
ana-miR1511a-3p comp390_c0_seq1 Noncoding
ana-miR156d-5p comp12590_c1_seq4 coding
ana-miR156h-5p comp12590_c1_seq2 coding
ana-miR164a-3p comp14847_c0_seq4 Noncoding
ana-miR164c-3p comp21396_c0_seq1 Noncoding
ana-miR396b-5p comp10567_c1_seq1 Noncoding
ana-miR396b-5p comp10567_c1_seq2 Noncoding
Identification of the Cold Stress-Responsive miRNAs in <italic>Ammopiptanthus nanus</italic>

We analyzed miRNA expression by comparing the abundance of each miRNA in the group CK and CA base on the high-throughput miRNA sequencing (miRNA-seq) (Figure 5A), and 39 cold-responsive miRNAs were identified, including 29 conserved miRNAs and ten nonconserved. Twenty-one conserved miRNAs were upregulated under cold acclimation, including ana-miR156a-5p, ana-miR159a-3p, ana-miR162a-3p, ana-miR166a-3p, ana-miR168a-5p, ana-miR390a-5p, ana-miR395a-3p, ana-miR399a-3p, ana-miR1507a-3p, ana-miR1515a-5p, ana-miR2118a-3p, ana-miR5211a-5p, and ana-miR5225a-3p, and eight conserved miRNAs were downregulated, that is, ana-miR1507a-5p, ana-miR156g-3p, ana-miR160a-3p, ana-miR171-1g-3p, ana-miR396b-3p, ana-miR398b-3p, ana-miR408a-5p, and ana-miR4415a-3p. Among the differentially expressed nonconserved miRNAs, eight were upregulated under cold acclimation, that is, ana-miRN8-5p, ana-miRN9-5p, ana-miRN18-5p, ana-miRN29-5p, ana-miRN6-3p, ana-miRN8-3p, ana-miRN15-3p, and ana-miRN39-3p and two were downregulated, that is, ana-miRN7-5p and ana-miRN19-3p.

Expression of cold-responsive miRNAs in A. nanus. (A) Expression of cold-responsive conserved miRNAs and non-conserved miRNAs. (B) qRT-PCR validation of the cold-responsive conserved miRNAs. (C) qRT-PCR validation of the cold-responsive nonconserved miRNAs. The Y-axis shows the relative expression levels of conserved miRNAs and non-conserved miRNAs, with U6 as an internal reference, averaged over three technical replicates. Values are expressed as means ± SD (n = 3).

qRT-PCR experiments were conducted to validate the results of high-throughput sequencing. Twelve conserved miRNAs, including ana-miR159a-3p, ana-miR168a-5p, ana-miR390a-5p, ana-miR399a-3p, ana-miR1507a-3p, ana-miR1515a-5p, ana-miR2118a-3p, ana-miR5211a-5p, ana-miR5225a-3p, ana-miR398b-3p, ana-miR408a-3p, and ana-miR4415a-3p, and five nonconserved miRNAs, that is, ana-miRN9-5p, ana-miRN18-5p, and ana-miRN39-3p, ana-miRN7-5p, and ana-miRN19-3p were selected for qRT-PCR analysis (Figures 5B,C). The change patterns of A. nanus miRNAs revealed by qRT-PCR analysis were consistent with those calculated based on the miRNA-seq results.

 <italic>Cis</italic>-Elements Involved in Abiotic Stress Responses in the Promoters of the Cold-Responsive miRNA Genes

Cis-acting elements have been reported to be associated with plant responses to abiotic stresses, including cold stress (Zhang et al., 2005). To identify the cis-acting elements associated with the expression patterns of the cold-responsive miRNAs, we analyzed the 1000-bp upstream promoter sequence of the 12 cold-responsive miRNAs, which represented 12 cold-responsive miRNA families with higher abundance, by using a plant promoter prediction software, PlantCARE (http://intra.psb.ugent.be: 8080/PlantCARE For) (Lescot et al., 2002), and consequently, several stress response- and hormone response–related cis-acting elements, including the low-temperature responsive elements (LTRs), were found.

Among the various cis-acting elements predicted from the promoters of the 12 cold-responsive miRNAs (Figure 6), the AU-rich element (ARE) was the most frequently detected cis-acting element, which was found in promoters of ten miRNAs, followed by ethylene response factor (ERE), appeared in promoters of seven miRNA genes. ARE has been shown to respond to hypoxic, low-temperature, and dehydration stress (Jang, 2016). Ethylene is an important gaseous phytohormone, playing a crucial role in plant growth and development and the response to biotic or abiotic stresses, and genes with ERE located in the promoter may be involved in ethylene response in plants (Rajesh et al., 2017). We also found several other cis-acting elements, including MBS, LTR, TC-rich repeat, TGA element, AuxRR core, GARE motif, P Box, ABRE, CGTCA motif, TGACG motif, TCA element, and W-box. ABREs were found in the promoter regions of most ABA-responsive genes and were reported to respond to abiotic stress in A. thaliana (Nagatoshi et al., 2018; Chang et al., 2019), and MBS is the MYB-binding site involved in drought inducibility (Lv et al., 2020).

Abiotic stress–responsive cis-elements predicted from the promoters of cold-responsive miRNA genes.

 Expression Analysis of the Targets of the Cold-Responsive miRNAs in <italic>Ammopiptanthus nanus</italic>

To investigate whether the expression of the predicted target negatively correlated with the expression of the corresponding miRNAs, we selected 12 predicted targets and analyzed their expression in cold-acclimated seedlings by qRT-PCR (Figure 7). After cold acclimation, seven miRNAs, including ana-miR159a-3p, ana-miR390a-5p, ana-miR2118a-3p, ana-miR5211a-5p, ana-miR5225a-3p, ana-miRN9-5p, and ana-miRN39-3p were upregulated, and the corresponding targets E3 UFM1-protein ligase-like protein (E3 UFM1), TAS3, ABC2 homolog 13 (ABC2), ATMRP3, Ankyrin repeat family protein (AFP), BTB-POZ and MATH domain 4 (BTB), and LFK L-fructokinase (LFK) were downregulated. The expression levels of ana-miR398b-3p, ana-miR408a-3p, ana-miR4415a-3p, ana-miRN7-5p, and ana-miRN19-3p were downregulated, and the corresponding targets LPR1 Cupredoxin superfamily protein (LPR1), UCC1 uclacyanin 1 (UCC1), L-AO, UBC23 ubiquitin-conjugating enzyme 23 (UBC23), PC-MYB1 homeodomain-like protein (PC-MYB1), auxin response factor 3 (ARF3), TIR1-NBS-LRR (TIR1), TIR2-NBS-LRR (TIR2), cold shock domain protein 1 (CSDP1), and DICER-LIKE2 (DCL2) were upregulated. The expression of these 12 selected targets displayed the opposite pattern to that of the corresponding miRNAs after cold acclimation.

qRT-PCR analysis of the targets of the 12 selected cold-responsive miRNAs. Red columns represent miRNAs and blue columns represent target genes. The Y-axis shows the relative expression with U6 or actin as the internal reference, and data are averaged over three technical replicates. Values are expressed as means ± SD.

 miR4415 is a Legume-Specific miRNA

To understand the roles of miRNAs in cold acclimation, miR4415 was selected for further functional study. In the miRBase database (miRBase 21), there are two mature sequences of miR4415, namely, gma-miR4415a and gma-miR4415b. The mature sequences of miR4415 and the corresponding precursor sequences from A. mongolicus and Astragalus membranaceus have been identified by small RNA sequencing in previous studies (Gao et al., 2018; Abla et al., 2019). miR4415 is absent in both A. thaliana and rice, and thus we speculated that miR4415 might be a lineage-specific miRNA. We used the BLAST program to compare the mature sequences of the abovementioned miR4415 with the genomic sequences of all plant species in Phytozome v12 to detect the mature sequences of miR4415 in other plant species. Consequently, 11 additional mature sequences of miR4415 (guide strand) that were not included in miRBase 21 were identified from nine plant species, that is, Glycine soja, Trifolium pratense, Cicer arietinum, Cajanus cajan, Vigna unguiculata, Vigna angularis, Medicago truncatula, Trifolium subterraneum, and Vigna radiata. It is worth noting that all miR4415 (guide strand) are located at 3p of the stem-loop precursors.

Multiple sequence alignment analysis indicated that the mature forms of miR4415 from different species differed significantly at nucleotide 21, and minor differences were also found at nucleotide 7, 9, and 18 (Figure 8A). Multiple sequence alignment analysis of the precursor sequences revealed that in addition to the strong conservation of miR4415-3p and miR4415-5p, some nucleotides in other sites of the precursor also showed a high degree of conservation, such as the two nucleotides before the 5′ end of miR4415-3p (Figure 8C). These conserved nucleotides may play important roles in the formation of hairpin structures.

Bioinformatics analysis of miR4415 sequences. Multiple sequence alignment of the mature forms of miR4415 (guide strand) (A). Phylogenetic tree of the precursors of miR4415 from several plant species (B). Multiple sequence alignment of miR4415 precursors from several plant species (C). cca-miR4415, miR4415 from Cajanus cajan; tsu-miR4415, miR4415 from Trifolium subterraneum; vun-miR4415, miR4415 from Vigna unguiculata; gso-miR4415a and gso-miR4415b, miR4415 from Glycine soja; tpr-miR4415a and tpr-miR4415b, miR4415 from Trifolium pratense; van-miR4415, miR4415 from Vigna angularis; mtr-miR4415, miR4415 from Medicago truncatula; vra-miR4415, miR4415 from Vigna radiata; pvu-miR4415, gma-miR4415a and gma-miR4415b, miR4415 from Glycine max; amo-miR4415, miR4415 from Ammopiptanthus mongolicus; ana-miR4415, miR4415 from Ammopiptanthus nanus.

Phylogenetic analysis using the precursors of the 12 miR4415 indicated that all miR4415 were clustered into three clades. The first clade contained van-miR4415, vra-miR4415, vun-miR4415, pvu-miR4415, and cca-miR4415, the second contained gso-miR4415a, gma-miR4415a, gso-miR4415b, gma-miR4415, and amo-miR4415, while mtr-miR4415, tsu-miR4415, tpr-miR4415a, and tpr-miR4415b were clustered into the third branch. Gma-miR4415a and gso-miR4415a in soybean and wild soybean were clustered into one clade, while gma-miR4415b and gso-miR4415b were clustered into another clade (Figure 8B). Thus, miR4415a and miR4415b were already present before the species differentiation in soybean and wild soybean.

 miR4415 Targets an L-Ascorbate Oxidase

An L-AO gene (EVM0030483.1) was predicted to be the target of ana-miR4415a-3p by using the online software psRNAtarget (Figure 9A). Dual-luciferase reporter assay, a reliable method to detect the targeting relationship between miRNAs and their possible targets, was conducted to validate the prediction experimentally. Compared with 62 SK + L-AO, the Luciferase/Renilla (LUC/REN) ratio was significantly decreased in the Pre-miR4415 + L-AO cells (Figure 9B), indicating that miR4415 negatively regulated the expression of L-AO (Figure 9C). The opposite change pattern between ana-miR4415a-3p and L-AO (EVM0030483.1) in cold-acclimated A. nanus revealed by qRT-PCR analysis also supported the targeting relationship of miR4415 and L-AO (EVM0030483.1) (Figure 7).

L-AO gene was targeted by ana-miR4415a-3p. Binding sites of ana-miR4415a-3p on L-AO as predicted by psRNAtarget (A). Dual-luciferase reporter vector construction for dual-luciferase reporter assay (B). Ratio of LUC/REN in the 62 SK + L-AO and Pre-miR4415 + L-AO cells (C). Values are expressed as means ± SD (n = 5). Student’s t-test, *p < 0.05.

 <italic>Ammopiptanthus nanus</italic> L-Ascorbate Oxidase is Located in the Apoplast

The L-AO (EVM0030483.1) protein was predicted to be located in the apoplast bioinformatically. To determine the subcellular localization of A. nanus L-AO (EVM0030483.1) experimentally, the L-AO (EVM0030483.1) gene was fused with the GFP gene and transiently expressed in Nicotiana tabacum. The GFP signal was mainly observed in the intercellular space, validating the apoplastic localization of L-AO protein (EVM0030483.1) (Figure 10A).

Subcellular localization of L-AO-GFP and analysis of apoplastic redox homeostasis in leaves of A. nanus under cold acclimation. (A) Subcellular localization of L-AO-GFP. Bars = 50 µm (B) AO activity. (C) Content of ASA. (D) Content of DHA. (E) ASA/DHA ratio. Values are expressed as means ± SD. Student’s t-test, *p < 0.05.

 L-AO Regulated the Contents of ASA and DHA in the Leaf Apoplast of <italic>Ammopiptanthus nanus</italic> Under Cold Acclimation

L-AO catalyzes the oxidation of ascorbic acid, which may affect the contents of ASA and DHA in the leaf apoplast of A. nanus. To assess the effect of upregulated L-AO on the concentration of ASA and DHA in the apoplast under cold acclimation, we extracted the apoplast fluid and measured the AO activity and ASA and DHA levels in the apoplast fluid of A. nanus before and after cold acclimation. Leaf extracts of A. nanus have MDH activity (4,755 u/g FW), while the MDH activity of the leaf apoplast liquid of A. nanus was not detected, indicating that the apoplast liquid we prepared had no cytoplasmic contamination and can be used for measurement of AO activity and the levels of ASA and DHA. The AO activity assay indicated that the AO activity in the apoplast of A. nanus was elevated after cold acclimation (Figure 10B). ASA and DHA measurements indicated a significant decrease in ASA (Figure 10C), a marked increase in DHA (Figure 10D), and a remarkable decrease in ASA/DHA ratio (Figure 10E) in the leaf apoplast liquid of A. nanus after cold acclimation.

 Discussion

A. nanus is a rare temperate evergreen broad-leaved plant, and it is of great significance to study its tolerance mechanism to freezing stress. As a class of endogenous regulatory RNA molecules, miRNAs not only participate in the regulation of plant growth and development but also in the response to environmental stresses (Song et al., 2019). Our study demonstrated that cold acclimation can improve the tolerance of A. nanus to the freezing stress, and miRNAs might play regulatory roles in cold acclimation in A. nanus. Analyzing the expression changes of miRNAs and their targets under cold acclimation is helpful in understanding the gene expression regulation network associated with the cold acclimation process. Our study systematically identified the miRNA and the corresponding targets that might be involved in cold acclimation and provided important data for understanding the molecular regulation mechanism of cold acclimation of A. nanus.

The Conserved miRNAs Predicted in <italic>Ammopiptanthus nanus</italic>

The completion of genome sequencing of A. nanus makes it possible to systematically identify the miRNA in A. nanus. Most of the identified miRNAs in A. nanus were conserved miRNAs, and the precursor sequences of these miRNAs were distributed on different assembled genome contigs. Of them, several conserved miRNA loci were clustered in genome contigs (Supplementary Table S2). For example, three precursors of miR2111, that is, ana-MIR2111a, ana-MIR2111b, and ana-MIR2111c were clustered in Contig00573, and the other three miR2111 precursors, that is, ana-MIR2111d, ana-MIR2111e, and ana-MIR2111f were clustered in Contig00782. There are nine precursors in the MIR395 family, and seven of them were clustered in Contig00419, including ana-MIR395a, ana-MIR395b, ana-MIR395c, ana-MIR395d, ana-MIR395e, ana-MIR395f, and ana-MIR395g, and the other two precursors, that is, ana-MIR395h and ana-MIR395i were distributed closely in Contig00478. Some of these clustered miRNA loci are probably generated via gene duplication.

Based on their conservation among plant species, the 41 conserved miRNA families in A. nanus can be classified into four groups. Ten miRNA families, including MIR156, MIR159, MIR160, MIR166, MIR167_1, MIR 171, MIR390, MIR395, MIR396, and MIR408, showed a high level of conservation among multiple plant species and are presented in Bryophyta, Lycopodiophyta, Coniferophyta, and Magnoliophyta, but not in Chlorophyta. This group contains the most conserved miRNAs in the plant kingdom. Group 2 comprised 12 miRNA families, that is, MIR162, MIR164, MIR168, MIR169, MIR172, MIR393, MIR394, MIR397, MIR398, MIR399, MIR2111, and MIR2118. These miRNAs are mainly identified in Magnoliophyta and they represent classical higher plant miRNAs. miRNAs in Group 3 (MIR403, MIR858, MIR1509, MIR1511, and MIR5225) are present only in eudicotyledons; hence, they might originate before the emergence of eudicotyledons. The remaining miRNA families, including MIR1507, MIR1515, MIR4415, MIR5037, and MIR862_2, were identified only in legume plants such as G. max and M. truncatula and might be legume-specific miRNAs.

 The Species-Specific miRNAs Predicted in <italic>Ammopiptanthus nanus</italic>

Besides the 68 nonconserved miRNA, we also identified approximately 264 miRNA candidates in A. nanus. Although their putative precursors can be folded into qualified hairpins (Axtell and Meyers, 2018), the reads mapped to 5′ arm or 3′ arm of the hairpins were not found or the abundance of 5p or 3p in sRNA libraries is less than ten; thus, these sRNA cannot be annotated as miRNA with high confidence. Maybe increasing the sequencing depth will result in identifying more novel miRNAs.

The sequencing frequencies of the 68 identified species-specific miRNAs varied greatly (Supplementary Table S3). Although the expression levels of most of the species-specific miRNAs were very low, some of them were expressed in a considerable amount. More than one-third of the species-specific miRNAs (24/64, 35.29%) were sequenced more than 100 times, and five species-specific miRNAs (ana-miRN8-5p, ana-miRN18-5p, ana-miRN4-3p, ana-miRN39-3p, and ana-miRN5-3p) have read numbers of more than 1,000 (Supplementary Table S3). Among all the identified species-specific miRNAs, ana-miRN4-3p was the most abundant, with a read number of more than 40,000 (Supplementary Table S3). Some of the species-specific miRNAs with high expression levels might have developed under environmental pressure, and thus might be involved in the adaptive evolution of A. nanus.

 The Noncanonical Targets Predicted by Degradome Sequencing Analysis

Degradome sequencing has been proven to be an efficient method to identify plant miRNA target genes (Ma X. et al., 2015). In A. nanus, many miRNA-target pairs identified by degradome sequencing have been reported in other plant species, including SQUAMOSA promoter binding protein-like (SPL) genes as the targets of miR156, the MYB genes as the targets of miR159, ARFs as the targets of miR160, miR167, miR393, endoribonuclease Dicer-like 1 (DCL1) as the target of miR162, nuclear factor Y genes (NF-Y) as the target of miR169, Integrase-type DNA-binding superfamily protein (AP2) as the target of miR172, ATP sulfurylase 1 (APS1) as the target of miR395, Laccase/Diphenol oxidase family protein (LAC) as the target of miR397, phosphate 2 as the target of miR399, Argonaute family protein (AGO2) as the target of miR403, NB-ARC domain-containing disease resistance protein as the targets of miR1507 and miR2118, and autoinhibited Ca 2+ -ATPase (ACA) as the target of miR5225. These results demonstrated the effectiveness of our degradome sequencing analysis.

Moreover, our degradome sequencing analysis also identified some noncanonical target genes of A. nanus miRNAs, which were not reported in other plant species. Among the 34 highly reliable target genes, which were classified in category 0, 21 had homologs in A. thaliana protein, and 13 cannot be annotated by the A. thaliana protein dataset, which were likely to be noncoding RNAs. Most of the noncanonical targets predicted in our study have no homology with the targets of the conserved miRNAs; thus, further investigations need to be conducted to understand the biological roles of these novel miRNA-target regulatory motifs.

 The miRNAs Were Involved in the Cold Acclimation in <italic>Ammopiptanthus nanus</italic>

A total of 39 cold-responsive miRNAs were identified in A. nanus, including 29 conserved miRNAs and ten nonconserved. The 29 conserved miRNAs belonged to 20 miRNA families such as MIR159, MIR168, MIR399, MIR5211, MIR5225, MIR398, MIR408, and MIR4415. Twelve conserved miRNAs that represent most of these cold-responsive conserved miRNAs and five nonconserved miRNAs with higher abundance were further validated by qRT-PCR. These cold-responsive miRNAs might participate in cold acclimation by regulating multiple biological processes, including ROS homeostasis (miR398, miR408, and miR4415), development (miR159 and miR390), defense (miR2118 and miR1507), calcium signaling (miR5225 and miR5211), and small RNA (sRNA)–mediated gene silencing (miR168 and miR1515) (Figure 11).

miRNA-mediated gene regulation network in response to cold acclimation in A. nanus leaves. miRNAs are indicated by red squares and targets by blue squares.

Calcium is a ubiquitous signal and an essential macronutrient in plants, and calcium signaling plays a vital role in environmental stress signal transduction. Calcium ATPases catalyze the hydrolysis of ATP coupled with the translocation of calcium from the cytosol into the endoplasmic reticulum or out of the cell and are essential components of stress signaling. In A. nanus, miR5225 and miR5211 were predicted to target genes encoding calcium-transporting ATPases, and both were upregulated after cold acclimation, indicating that these miRNAs might be involved in cold acclimation by regulating the calcium signaling pathway.

Three miRNAs, that is, miR398, miR408, and miR4415 might be involved in cold acclimation by regulating their targets, copper/zinc superoxide dismutase gene (CSD), ceruloplasmin/laccase, and L-AO, respectively. It is intriguing that all three targets encode multicopper oxidases (MCO). Multicopper oxidases are a class of enzymes that catalyze the oxidation of various substrates via the reduction of O2 to H2O without releasing activated oxygen species (Komori and Higuchi, 2015). Of them, ascorbate oxidase (AO, EC 1.10.3.3) is localized in the cell wall and widespread in the plant kingdom (Sanmartin et al., 2007). L-AO oxidizes AsA to DHA by monodehydroascorbate (MDHA) with the participation of O2− and H2O2 while reducing oxygen to water, thereby regulating the redox state of the apoplast.

Laccases are classified as benzene diol oxygen reductases (EC 1.10.3.2) and are involved in lignin degradation and detoxification of lignin-derived product wound repair, disease resistance, stress resistance, and pigment synthesis in plants (Ranocha et al., 1999). Ceruloplasmin plays a vital role in copper metabolism and is also involved in scavenging anion radicals in plants (Bielli and Calabrese, 2002). miR408 was found to be downregulated in our study, which is consistent with previous reports in wheat and grape (Wang et al., 2014) but different from the observations in A. thaliana (Liu et al., 2008). Overexpression of miR408 promoted drought tolerance in chickpea (Hajyzadeh et al., 2015) and cold and oxidative stresses tolerance in A. thaliana (Ma C. et al., 2015).

CSDs are copper- and zinc-dependent enzymes that catalyze the conversion of superoxide radicals to H2O2, which can subsequently be reduced to H2O, and CSDs are proposed to play essential roles in scavenging the excessive ROS generated in plants under stressful conditions. Several reports have shown that miR398 plays an important role in plant antioxidant stress by regulating the expression of its target gene CSD (Jones Rhoades and Bartel, 2004; Sunkar and Zhu, 2004), although the expression pattern of miR398 was distinctly different. It has also been found that miR398 is upregulated in wheat but downregulated in rice when subjected to heat stress (Sailaja et al., 2014). Our data showed that miR398, miR408, and miR4415 were downregulated and their targets were upregulated in cold-acclimated A. nanus seedlings, and these miRNAs might contribute to cold acclimation by regulating ROS homeostasis in A. nanus.

miR159 has been reported to be induced by various abiotic or biotic stresses (Millar et al., 2019) and transient or mild upregulation of miR159 was observed in A. thaliana under cold stress (Liu et al., 2008). In A. nanus, we found that ana-miR159a-3p was upregulated under cold acclimation. The targets of ana-miR159a-3p in A. nanus encoded an MYB65 and an E3 UFM1-protein ligase–like protein (E3 UFM1). MYB TFs were involved in the regulation of various biological processes in plant growth and development in several plant species such as A. thaliana (Dubos et al., 2010; Chen et al., 2013) and rice (Yang et al., 2012). E3 ligase modulates plant abiotic stress responses through posttranslational modifications (Sun and Chen, 2004; Lyzenga and Stone, 2012). Ana-miR159a-3p may be involved in cold acclimation of A. nanus by negatively regulating ATMYB65 and E3 UFM1.

 MiR4415 is involved in Cold Acclimation in <italic>Ammopiptanthus nanus</italic> by Regulating the Apoplastic Redox Status by Targeting an L-AO

Most miRNAs are conserved, to some extent, among different species, but not all miRNAs are highly conserved. Some miRNAs exist only in a specific evolutionary branch, and these miRNAs are called lineage-specific miRNAs. The lineage-specific miRNAs were evolved/developed in a specific evolutionary branch, hence contributing to developmental novelties during evolution (Zhang et al., 2008), and might be closely associated with some unique regulatory mechanism of the plants in this evolutionary branch. Our results demonstrated that miR4415 is a legume-specific miRNA family, which has only been reported in legumes such as G. max (Lei et al., 2019), G. soja (Zeng et al., 2012), Cajanus cajan (Nithin et al., 2017), Phaseolus vulgaris (Formey et al., 2015), A. membranaceus (Abla et al., 2019), and A. mongolicus (Gao et al., 2018).

miR4415 has been reported to be induced by aluminum treatment and drought stress but downregulated by biotic stress in G. soja (Zeng et al., 2012). miR4415 was observed to be downregulated by cold stress treatment in A. nanus and A. membranaceus (Abla et al., 2019) and by osmotic stress in A. mongolicus (Gao et al., 2018). These results indicated its involvement in the abiotic and biotic stress response.

AO was a key enzyme regulating the ROS status of the apoplast and an important regulatory enzyme for plant growth, development, environmental stress perception, and stress signal transduction (Foyer et al., 2020). AO can catalyze the oxidation of ASA to unstable MDHA and then form dehydroascorbic acid through nonenzymatic disproportionation. The redox status of the apoplast is largely determined by AO (Karpinska et al., 2018), and the change of ASA/DHA ratio caused by environmental stress may regulate gene expression and protein abundance through the change of redox state of the apoplast.

In A. nanus, expression of the miR4415 precursor was inhibited after cold acclimation, leading to downregulation of the ana-miR4415a-3p and upregulation of the L-AO gene. An elevated level of L-AO promotes the oxidation of ASA to form DHA and lowers the AsA/DHA ratio, resulting in alteration of the redox state of the apoplast. Change in the apoplastic redox state regulates the tolerance to freezing stress in A. nanus seedlings (Figure 12). A similar regulatory pathway has been reported in rice, in which miR528 regulated the apoplastic ROS state by targeting the L-AO gene, leading to enhanced tolerance to salt stress (Wang et al., 2021) and viral infection (Wu et al., 2017).

Proposed model of the regulatory role of miR4415 in cold acclimation in A. nanus. Under cold acclimation, expression of the miR4415 precursor was inhibited, resulting in downregulation of the ana-miR4415a-3p and upregulation of L-AO. Upregulated L-AO catalyzes the oxidation of ASA to form DHA, lowering the AsA/DHA ratio, thereby changing the redox status of the apoplast. Change in apoplastic redox state regulates the tolerance to freezing stress in A. nanus seedlings.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, SRR10906482https://www.ncbi.nlm.nih.gov/, SRR10906483https://www.ncbi.nlm.nih.gov/, SRR10906484https://www.ncbi.nlm.nih.gov/, SRR10906485https://www.ncbi.nlm.nih.gov/, SRR10906486https://www.ncbi.nlm.nih.gov/, SRR10906487https://www.ncbi.nlm.nih.gov/, SRR17711209.

 Author Contributions

YiZ and FG contributed to the conception and design of the study. MZ performed the statistical analysis. MZ, JT, XW, YaZ, and FG performed the investigation. MZ and FG wrote the first draft of the manuscript. FG and YiZ wrote sections of the manuscript. All the authors contributed to manuscript revision and read and approved the submitted version.

 Funding

This study was funded by the National Natural Science Foundation of China (Grant Number 31770363 and 31670335), the Ministry of Education of China under “Double First-Class” projects (Number Yldxxk201819), Key Laboratory of Ecology and Environment in Minority Areas (Minzu University of China), the National Ethnic Affairs Commission (No. KLEEMA202106), and Beijing Advanced Discipline for Mass Spectrometry Imaging and Metabolomics (No. 104-01900403).

 Conflict of Interest

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.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.870446/full#supplementary-material

Primers used for qRT-PCR analysis.

Conserved miRNAs identified from A. nanus.

Species-specific miRNAs identified in A. nanus.

Targets of the conserved miRNAs predicted in A. nanus by psRNATarget.

Targets of the species-specific miRNAs predicted in A. nanus by psRNATarget.

Targets of the miRNAs identified in A. nanus by degradome analysis.

 References Abla M. Sun H. Li Z. Wei C. Gao F. Zhou Y. (2019). Identification of miRNAs and Their Response to Cold Stress in Astragalus Membranaceus . Biomolecules 9, 182. 10.3390/biom9050182 Addo-Quaye C. Miller W. Axtell M. J. (2009). CleaveLand: A Pipeline for Using Degradome Data to Find Cleaved Small RNA Targets. Bioinformatics 25, 130–131. 10.1093/bioinformatics/btn604 Axtell M. J. Meyers B. C. (2018). Revisiting Criteria for Plant microRNA Annotation in the Era of Big Data. Plant Cell 30, 272–284. 10.1105/tpc.17.00851 Bielli P. Calabrese L. (2002). Structure to Function Relationships in Ceruloplasmin: A 'Moonlighting' Protein. Cell Mol. Life Sci. (Cmls) 59, 1413–1427. 10.1007/s00018-002-8519-2 Bredow M. Walker V. K. (2017). Ice-Binding Proteins in Plants. Front. Plant Sci. 8, 2153. 10.3389/fpls.2017.02153 Cao S. Wang Y. Li X. Gao F. Feng J. Zhou Y. (2020). Characterization of the AP2/ERF Transcription Factor Family and Expression Profiling of DREB Subfamily under Cold and Osmotic Stresses in Ammopiptanthus Nanus . Plants 9, 455. 10.3390/plants9040455 Cao S. Wang Y. Li Z. Shi W. Gao F. Zhou Y. (2019). Genome-Wide Identification and Expression Analyses of the Chitinases under Cold and Osmotic Stress in Ammopiptanthus Nanus . Genes 10, 472. 10.3390/genes10060472 Chang H.-C. Tsai M.-C. Wu S.-S. Chang I.-F. (2019). Regulation of ABI5 Expression by ABF3 during Salt Stress Responses in Arabidopsis Thaliana . Bot. Stud. 60, 16. 10.1186/s40529-019-0264-z Chen C. Chen H. Zhang Y. Thomas H. R. Frank M. H. He Y. (2020). TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol. Plant 13, 1194–1202. 10.1016/j.molp.2020.06.009 Chen Y. Chen Z. Kang J. Kang D. Gu H. Qin G. (2013). AtMYB14 Regulates Cold Tolerance in Arabidopsis. Plant Mol. Biol. Rep. 31, 87–97. 10.1007/s11105-012-0481-z Clément T. Salone V. Rederstorff M. (2015). Dual Luciferase Gene Reporter Assays to Study miRNA Function. Methods Mol. Biol. 1296, 187–198. 10.1007/978-1-4939-2547-6_17 Dai X. Zhuang Z. Zhao P. X. (2018). psRNATarget: A Plant Small RNA Target Analysis Server (2017 Release). Nucleic Acids Res. 46, W49–W54. 10.1093/nar/gky316 Dewey C. N. (2010). RNA-Seq Gene Expression Estimation with Read Mapping Uncertainty. Bioinformatics 26, 493–500. 10.1093/bioinformatics/btp692 Ding Y. Shi Y. Yang S. (2020). Molecular Regulation of Plant Responses to Environmental Temperatures. Mol. Plant 13, 544–564. 10.1016/j.molp.2020.02.004 Dubos C. Stracke R. Grotewold E. Weisshaar B. Martin C. Lepiniec L. (2010). MYB Transcription Factors in Arabidopsis . Trends Plant Sci. 15, 573–581. 10.1016/j.tplants.2010.06.005 Formey D. Iñiguez L. P. Peláez P. Li Y.-F. Sunkar R. Sánchez F. (2015). Genome-Wide Identification of the Phaseolus Vulgaris sRNAome Using Small RNA and Degradome Sequencing. BMC Genomics 16, 423. 10.1186/s12864-015-1639-5 Foyer C. H. Kyndt T. Hancock R. D. (2020). Vitamin C in Plants: Novel Concepts, New Perspectives, and Outstanding Issues. Antioxid. Redox Signaling 32, 463–485. 10.1089/ars.2019.7819 Gao F. Wang X. Li X. Xu M. Li H. Abla M. (2018). Long-Read Sequencing and De Novo Genome Assembly of Ammopiptanthus Nanus, a Desert Shrub. Gigascience 7, 074. 10.1093/gigascience/giy074 Gao F. Wang N. Li H. Liu J. Fu C. Xiao Z. (2016). Identification of Drought-Responsive microRNAs and Their Targets in Ammopiptanthus Mongolicus by Using High-Throughput Sequencing. Sci. Rep. 6, 34601. 10.1038/srep34601 Hajyzadeh M. Turktas M. Khawar K. M. Unver T. (2015). miR408 Overexpression Causes Increased Drought Tolerance in Chickpea. Gene 555, 186–193. 10.1016/j.gene.2014.11.002 Jang J.-C. (2016). Arginine-Rich Motif-Tandem CCCH Zinc finger Proteins in Plant Stress Responses and Post-Transcriptional Regulation of Gene Expression. Plant Sci. 252, 118–124. 10.1016/j.plantsci.2016.06.014 Jones-Rhoades M. W. Bartel D. P. (2004). Computational Identification of Plant MicroRNAs and Their Targets, Including a Stress-Induced miRNA. Mol. Cel 14, 787–799. 10.1016/j.molcel.2004.05.027 Kang Y.-J. Yang D.-C. Kong L. Hou M. Meng Y.-Q. Wei L. (2017). CPC2: A Fast and Accurate Coding Potential Calculator Based on Sequence Intrinsic Features. Nucleic Acids Res. 45, W12–W16. 10.1093/nar/gkx428 Karpinska B. Zhang K. Rasool B. Pastok D. Morris J. Verrall S. R. (2018). The Redox State of the Apoplast Influences the Acclimation of Photosynthesis and Leaf Metabolism to Changing Irradiance. Plant Cel Environ 41, 1083–1097. 10.1111/pce.12960 Khan T. A. Fariduddin Q. Yusuf M. (2017). Low-temperature Stress: Is Phytohormones Application a Remedy? Environ. Sci. Pollut. Res. 24, 21574–21590. 10.1007/s11356-017-9948-7 Komori H. Higuchi Y. (2015). Structural Insights into the O2 Reduction Mechanism of Multicopper Oxidase. J. Biochem. 158, 293–298. 10.1093/jb/mvv079 Kozomara A. Birgaoanu M. Griffiths-Jones S. (2019). miRBase: From microRNA Sequences to Function. Nucleic Acids Res. 47, D155–D162. 10.1093/nar/gky1141 Kozomara A. Griffiths-Jones S. (2014). miRBase: Annotating High Confidence microRNAs Using Deep Sequencing Data. Nucl. Acids Res. 42, D68–D73. 10.1093/nar/gkt1181 Krützfeldt J. Stoffel M. (2006). MicroRNAs: A New Class of Regulatory Genes Affecting Metabolism. Cel Metab. 4, 9–12. 10.1016/j.cmet.2006.05.009 Lei P. Han B. Wang Y. Zhu X. Xuan Y. Liu X. (2019). Identification of microRNAs that Respond to Soybean Cyst Nematode Infection in Early Stages in Resistant and Susceptible Soybean Cultivars. Int. J. Mol. Sci. 20, 5634. 10.3390/ijms20225634 Lescot M. Dehais P. Thijs G. Marchal K. Moreau Y. Van De Peer Y. (2002). PlantCARE, a Database of Plant Cis-Acting Regulatory Elements and a portal to Tools for In Silico Analysis of Promoter Sequences. Nucleic Acids Res. 30, 325–327. 10.1093/nar/30.1.325 Liu C. T. Wang W. Mao B. G. Chu C. C. (2018). Cold Stress Tolerance in Rice: Physiological Changes, Molecular Mechanism, and Future Prospects. Yi Chuan 40, 171–185. 10.16288/j.yczz.18-007 Liu H.-H. Tian X. Li Y.-J. Wu C.-A. Zheng C.-C. (2008). Microarray-Based Analysis of Stress-Regulated microRNAs in Arabidopsis Thaliana . Rna 14, 836–843. 10.1261/rna.895308 Liu Y. Dang P. Liu L. He C. (2019). Cold Acclimation by the CBF-COR Pathway in a Changing Climate: Lessons from Arabidopsis thaliana . Plant Cel Rep 38, 511–519. 10.1007/s00299-019-02376-3 Lv D.-K. Bai X. Li Y. Ding X.-D. Ge Y. Cai H. (2010). Profiling of Cold-Stress-Responsive miRNAs in rice by Microarrays. Gene 459, 39–47. 10.1016/j.gene.2010.03.011 Lv L.-M. Zuo D.-Y. Wang X.-F. Cheng H.-L. Zhang Y.-P. Wang Q.-L. (2020). Genome-Wide Identification of the Expansin Gene Family Reveals that Expansin Genes Are Involved in Fibre Cell Growth in Cotton. BMC Plant Biol. 20, 223. 10.1186/s12870-020-02362-y Lyzenga W. J. Stone S. L. (2012). Abiotic Stress Tolerance Mediated by Protein Ubiquitination. J. Exp. Bot. 63, 599–616. 10.1093/jxb/err310 Ma C. Burd S. Lers A. (2015a). miR408is Involved in Abiotic Stress Responses in Arabidopsis. Plant J. 84, 169–187. 10.1111/tpj.12999 Ma F. Liu Z. Huang J. Li Y. Kang Y. Liu X. (2019). High-Throughput Sequencing Reveals microRNAs in Response to Heat Stress in the Head Kidney of Rainbow trout (Oncorhynchus Mykiss). Funct. Integr. Genomics 19, 775–786. 10.1007/s10142-019-00682-3 Ma X. Tang Z. Qin J. Meng Y. (2015b). The Use of High-Throughput Sequencing Methods for Plant microRNA Research. RNA Biol. 12, 709–719. 10.1080/15476286.2015.1053686 Megha S. Basu U. Kav N. N. V. (2018). Regulation of Low Temperature Stress in Plants by microRNAs. Plant Cel Environ 41, 1–15. 10.1111/pce.12956 Millar A. A. Lohe A. Wong G. (2019). Biology and Function of miR159 in Plants. Plants 8, 255. 10.3390/plants8080255 Morin R. D. Aksay G. Dolgosheina E. Ebhardt H. A. Magrini V. Mardis E. R. (2008). Comparative Analysis of the Small RNA Transcriptomes of Pinus Contorta and Oryza Sativa . Genome Res. 18, 571–584. 10.1101/gr.6897308 Nagatoshi Y. Fujita M. Fujita Y. (2018). Casein Kinase 2 α and β Subunits Inversely Modulate ABA Signal Output in Arabidopsis Protoplasts. Planta 248, 571–578. 10.1007/s00425-018-2919-5 Nithin C. Thomas A. Basak J. Bahadur R. P. (2017). Genome-Wide Identification of miRNAs and lncRNAs in Cajanus Cajan . BMC Genomics 18, 878. 10.1186/s12864-017-4232-2 O'leary B. M. Rico A. Mccraw S. Fones H. N. Preston G. M. (2014). The Infiltration-Centrifugation Technique for Extraction of Apoplastic Fluid from Plant Leaves Using Phaseolus Vulgaris as an Example. J. Vis. Exp., 52113. 10.3791/52113 Rajesh K. S. Bharath B. R. Rao C. V. Bhat K. I. Bhat K. S. C. Bhat P. (2017). Neutralization of Naja naja Venom Induced Lethality, Edema and Myonecrosis by Ethanolic Root Extract of Coix Lacryma-Jobi . Toxicol. Rep. 4, 637–645. 10.1016/j.toxrep.2017.11.005 Ranocha P. Mcdougall G. Hawkins S. Sterjiades R. Borderies G. Stewart D. (1999). Biochemical Characterization, Molecular Cloning and Expression of Laccases - a Divergent Gene Family - in Poplar. Eur. J. Biochem. 259, 485–495. 10.1046/j.1432-1327.1999.00061.x Sailaja B. Voleti S. R. Subrahmanyam D. Sarla N. Prasanth V. V. Bhadana V. P. (2014). Prediction and Expression Analysis of miRNAs Associated with Heat Stress in Oryza Sativa . Rice Sci. 21, 3–12. 10.1016/s1672-6308(13)60164-x Sanmartin M. Pateraki I. Chatzopoulou F. Kanellis A. K. (2007). Differential Expression of the Ascorbate Oxidase Multigene Family during Fruit Development and in Response to Stress. Planta 225, 873–885. 10.1007/s00425-006-0399-5 Song C. Wang C. Zhang C. Korir N. K. Yu H. Ma Z. (2010). Deep Sequencing Discovery of Novel and Conserved microRNAs in Trifoliate orange (Citrus Trifoliata). BMC Genomics 11, 431. 10.1186/1471-2164-11-431 Song X. Li Y. Cao X. Qi Y. (2019). MicroRNAs and Their Regulatory Roles in Plant-Environment Interactions. Annu. Rev. Plant Biol. 70, 489–525. 10.1146/annurev-arplant-050718-100334 Sun L. Chen Z. J. (2004). The Novel Functions of Ubiquitination in Signaling. Curr. Opin. Cel Biol. 16, 119–126. 10.1016/j.ceb.2004.02.005 Sunkar R. Zhu J.-K. (2004). Novel and Stress-Regulated MicroRNAs and Other Small RNAs from Arabidopsis[W]. Plant Cell 16, 2001–2019. 10.1105/tpc.104.022830 Szittya G. Moxon S. Santos D. M. Jing R. Fevereiro M. P. Moulton V. (2008). High-Throughput Sequencing of Medicago Truncatula Short RNAs Identifies Eight New miRNA Families. BMC Genomics 9, 593. 10.1186/1471-2164-9-593 Tang W. Thompson W. A. (2019). OsmiR528 Enhances Cold Stress Tolerance by Repressing Expression of Stress Response-Related Transcription Factor Genes in Plant Cells. Curr. Genomics 20, 100–114. 10.2174/1389202920666190129145439 Varkonyi-Gasic E. Wu R. Wood M. Walton E. F. Hellens R. P. (2007). Protocol: A Highly Sensitive RT-PCR Method for Detection and Quantification of microRNAs. Plant Methods 3, 12. 10.1186/1746-4811-3-12 Wang B. Sun Y.-F. Song N. Wei J.-P. Wang X.-J. Feng H. (2014). MicroRNAs Involving in Cold, Wounding and Salt Stresses in Triticum Aestivum L. Plant Physiol. Biochem. 80, 90–96. 10.1016/j.plaphy.2014.03.020 Wang M. Guo W. Li J. Pan X. Pan L. Zhao J. (2021). The miR528-AO Module Confers Enhanced Salt Tolerance in rice by Modulating the Ascorbic Acid and Abscisic Acid Metabolism and ROS Scavenging. J. Agric. Food Chem. 69, 8634–8648. 10.1021/acs.jafc.1c01096 Wu J. Yang R. Yang Z. Yao S. Zhao S. Wang Y. (2017). ROS Accumulation and Antiviral Defence Control by microRNA528 in rice. Nat. Plants 3, 16203. 10.1038/nplants.2016.203 Xu S. Ha C. H. Wang W. Xu X. Yin M. Jin F. Q. (2016). PECAM1 Regulates Flow-Mediated Gab1 Tyrosine Phosphorylation and Signaling. Cell Signal. 28, 117–124. 10.1016/j.cellsig.2015.12.007 Yang A. Dai X. Zhang W.-H. (2012). A R2R3-Type MYB Gene, OsMYB2, Is Involved in Salt, Cold, and Dehydration Tolerance in rice. J. Exp. Bot. 63, 2541–2556. 10.1093/jxb/err431 Yang G. Rhodes D. Joly R. (1996). Effects of High Temperature on Membrane Stability and Chlorophyll Fluorescence in Glycinebetaine-Deficient and Glycinebetaine-Containing maize Lines. Funct. Plant Biol. 23, 437. 10.1071/pp9960437 Yang G. Zou H. Wu Y. Liu H. Yuan Y. (2011). Identification and Characterisation of Candidate Genes Involved in Chilling Responses in Maize (Zea mays L.). Plant Cel Tiss Organ. Cult 106, 127–141. 10.1007/s11240-010-9900-8 Yasar F. Ellialtioglu S. Yildiz K. (2008). Effect of Salt Stress on Antioxidant Defense Systems, Lipid Peroxidation, and Chlorophyll Content in green Bean. Russ. J. Plant Physiol. 55, 782–786. 10.1134/s1021443708060071 Yin Z. Wang N. Duan W. Pi L. Shen D. Dou D. (2021). Phytophthora Capsici CBM1‐Containing Protein CBP3 Is an Apoplastic Effector with Plant Immunity‐Inducing Activity. Mol. Plant Pathol. 22, 1358–1369. 10.1111/mpp.13116 Yoo S.-D. Cho Y.-H. Sheen J. (2007). Arabidopsis Mesophyll Protoplasts: A Versatile Cell System for Transient Gene Expression Analysis. Nat. Protoc. 2, 1565–1572. 10.1038/nprot.2007.199 Yu Y. Zhang Y. Chen X. Chen Y. (2019). Plant Noncoding RNAs: Hidden Players in Development and Stress Responses. Annu. Rev. Cel Dev. Biol. 35, 407–431. 10.1146/annurev-cellbio-100818-125218 Zeng Q.-Y. Yang C.-Y. Ma Q.-B. Li X.-P. Dong W.-W. Nian H. (2012). Identification of Wild Soybean miRNAs and Their Target Genes Responsive to Aluminum Stress. BMC Plant Biol. 12, 182. 10.1186/1471-2229-12-182 Zhang R. Wang Y.-Q. Su B. (2008). Molecular Evolution of a Primate-Specific microRNA Family. Mol. Biol. Evol. 25, 1493–1502. 10.1093/molbev/msn094 Zhang W. Ruan J. Ho T.-H. D. You Y. Yu T. Quatrano R. S. (2005). Cis-Regulatory Element Based Targeted Gene Finding: Genome-Wide Identification of Abscisic Acid- and Abiotic Stress-Responsive Genes in Arabidopsis Thaliana . Bioinformatics 21, 3074–3081. 10.1093/bioinformatics/bti490 Zhang Y. Cao Y. Zheng H. Feng W. Qu J. Fu F. (2021). Ectopic Expression of Antifreeze Protein Gene from Ammopiptanthus Nanus Confers Chilling Tolerance in maize. Crop J. 9, 924–933. 10.1016/j.cj.2020.08.011 Zhou M. Tang W. (2019). MicroRNA156 Amplifies Transcription Factor-Associated Cold Stress Tolerance in Plant Cells. Mol. Genet. Genomics 294, 379–393. 10.1007/s00438-018-1516-4