\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Multi-omics analysis reveals a link between Brassica -specific miR1885
and rapeseed tolerance to low temperature}



\author[1]{Xiang Yu}%
\author[1]{Pengfei Xu}%
\author[2]{Wenting Zhang}%
\author[3]{Xuan Wang}%
\author[4]{Yantao Zhu}%
\author[1]{Wanqi Liang}%
\author[2]{Yuke He}%
\affil[1]{Shanghai Jiao Tong University School of Life Sciences and Biotechnology}%
\affil[2]{Institute of Plant Physiology and Ecology Shanghai Institute for Biological Sciences}%
\affil[3]{Shanghai Normal University}%
\affil[4]{Hybrid Rapeseed Research Center of Shaanxi Province}%


\vspace{-1em}



  
  \date{June 05, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
\emph{Brassica} crops include various edible vegetable and plant oil
crops, and their production is limited by low temperature beyond their
tolerant capability. The key regulators of low-temperature resistance in
\emph{Brassica} remain largely unexplored. To identify
post-transcriptional regulators of plant response to low temperature, we
performed small RNA profiling, and found that 16 known miRNAs were
responsive to cold treatment in \emph{Brassica rapa}. The cold response
of seven of those miRNAs were further confirmed by qRT-PCR and/or
northern blotting analyses. In parallel, a genome-wide association study
of 220 accessions of \emph{Brassica napus} identified four candidate
MIRNA genes, all of which were cold-responsive, at the loci associated
with low temperature resistance. Specifically, these large-scale data
analyses revealed a link between miR1885 and the plant response to low
temperature in both \emph{B. rapa} and \emph{B. napus}. Using 5' rapid
amplification of cDNA ends approach, we validated that miR1885 can
cleave its putative target transcripts, \emph{Bn.TIR.A09} and
\emph{Bn.TNL.A03,} in \emph{B. napus}. Furthermore, overexpression of
miR1885 in Semi-winter type \emph{B. napus} decreased the mRNA abundance
of \emph{Bn.TIR.A09} and \emph{Bn.TNL.A03}, resulting in increased
sensitivity to low temperature. Knocking down of miR1885 in Spring type
\emph{B. napus} led to increased mRNA abundance of its targets and
improved rapeseed tolerance to low temperature. Together, our results
suggested that the loci of miR1885 and its targets could be potential
candidates for the molecular breeding of low temperature-tolerant Spring
type \emph{Brassica} crops.%
\end{abstract}%



\sloppy


\textbf{INTRODUCTION}

\emph{Brassica} crop species include \emph{Brassica napus}
(rapeseed),\emph{Brassica rapa} , \emph{Brassica oleracea} ,
\emph{Brassica juncea} , and \emph{Brassica nigra} . Among them,
\emph{B. napus} arose approximately 7500 years ago by
allopolyploidization between its ancestors, \emph{B. rapa} and \emph{B.
oleracea} (Chalhoub et al., 2014). \emph{B. rapa} and \emph{Arabidopsis
thaliana} share a common ancestor, but \emph{B. rapa} has a larger
genome than that of \emph{A. thaliana} due to genome triplication (Wang
et al., 2011). Both \emph{B. napus} and \emph{B. rapa} are important
food sources. \emph{B. napus}mainly provides edible plant oil and
\emph{B. rapa} is a valuable oil and vegetable crop. However,
\emph{Brassica} crops are negatively affected by low temperature,
especially those that encounter frigid environment beyond their
cold-tolerant capability. Worldwide, extensive economic losses of
\emph{Brassica} crops are attributed to low temperature, often in
combination with other abiotic or biotic stresses (Sanghera et al.,
2011). Brassica crops mainly have three ecotypes: winter, spring, and
semi-winter types. The winter and semi-winter type Brassica require
vernalization, which are more tolerant to low temperature as compared to
spring type Brassica (Gomez-Campo and Prakash, 1999). Low temperature
includes coldness (\textgreater{} 0  and \textless{} 10 ) and freeze
(\textless{} 0 ) (Ding et al., 2019). At the early seeding stage,
long-term low temperature stress (4 \selectlanguage{ngerman}°C) would restrict the growth of
spring Brassica crops and even cause plant damage and yield loss (O'
neill et al., 2019). It is important to identify key regulators
contributing to plant response to cold stress for improving the cold
resistance of spring type Brassica crops.

The vernalization of semi-winter type Brassica crops require moderate
low temperature. However, long time exposure to low temperature in
winter leads to plant damage at vegetative stage and even lead to yield
loss (Liao and Guan, 2001; Zhang et al., 2015). The Brassica crops are
usually sown in autumn after the harvest of rice in Yangtze River basin,
China (Cong et al., 2019). In recent years, the delay of rice harvest
usually led to the postpone of rapeseed sowing resulting in poor
Brassica crop seedings establishment due to low temperature. Lacking
enough biomass of Brassica crops in winter cause the plants more
susceptible to cold or freezing stress (Luo et al., 2019; Zhang et al.,
2012). Moreover, some Spring type accessions would die due to long-term
low temperature and decrease of seed yield (Ozer, 2003).

Previous studies have identi\selectlanguage{english}fied a number of genes involved in the plant
response to low temperature (Ding et., 2019). Cold stress induced the
expression of the transcription factor gene \emph{CBF} (encoding
C-repeat binding factor), and CBF triggered the expression of
\emph{COR}(COLD RESPONSIVE) genes (Shi et al., 2018). CBF/COR-dependent
factors allowed plants to withstand subsequent freezing stress (McClung
and Davis, 2010). In addition, many important genes, which contributed
to low-temperature response in \emph{CBF} -independent manner, were also
found (Bolt et al., 2017; Liu et al., 2018). Although breeding for
low-temperature resistance is possible, the genetic improvement of crops
is still hindered by the lack of elite genetic resources related to
low-temperature resistance.

MicroRNAs (miRNAs) are regulators of plant adaptation to abiotic or
biotic stress (Biggar and Storey, 2015). They are approximately 20--24
nucleotides (nt) in length and originate from pre-miRNAs with stem-loop
structures (Song et al., 2019). The mature miRNAs are loaded into the
RNA-induced silencing complex that regulate target transcripts involved
in various aspects of plant growth, development, and stress responses,
including response to cold or freezing stress (Chen, 2009; Anjali and
Sabu, 2020). For instance, overexpression of miR156 increased cold
tolerance in rice, \emph{A. thaliana,} and pine (Zhou and Tang, 2019).
In rice, OsmiR319b down-regulated \emph{OsPCF6} and \emph{OsTCP21} ,
resulting in enhanced cold tolerance (Wang et al., 2014), and
overexpressing OsmiR1320 improved cold tolerance, although OsmiR1320 was
repressed under cold stress (Sun et al., 2022). In \emph{A.
thaliana,}plants knocking-down miR165/166 displayed a drought- and
cold-resistant phenotype (Yan et al., 2016). Other known cold-responsive
miRNAs included miR172, miR396, miR845, and miR168 (Gupta et al., 2014;
Zeng et al., 2018, Chen et al., 2019). However, miRNAs in response to
low temperature stress have not been well studied in \emph{Brassica}
crops.

Plant \emph{R} (RESISTANCE) gene encodes immune protein involved in the
response to pathogen infection (Jones and Dangl, 2006). Most of the
identified \emph{R} genes contain a nucleotide-binding site-leucine-rich
repeat (NBS-LRR) domain with a coiled-coil (CC) or a Toll/Interleukin-1
receptor (TIR) domain at the amino terminus (Dangl and Jones, 2001).
Emerging evidence has shown that cold stress and pathogen attack induced
the expression of some common genes, such as
\emph{PR}(PATHOGENESIS-RELATED) genes, which played roles in the
response to pathogen infection and low temperature stress (Snider et
al., 2000; Seo et al., 2008). Some R genes negatively regulated plants
tolerance to low temperature. In \emph{A. thaliana} , the \emph{rpp4-1d}
mutant plant with a gain-of-function mutation in \emph{RPP4} (TIR-NB-LRR
protein PERONOSPORA PARASITICA 4) was chilling-sensitive (Huang et al.,
2010). Mutations in \emph{HSP90} suppressed the chilling-sensitive
phenotype of\emph{rpp4-1d} mutant plant, and HSP90 interacted with the
RPP4 protein in planta (Bao et al., 2014). Interestingly, a significant
increase in\emph{HSP90} mRNA was detected in \emph{B. napus} exposed to
low temperature (Krishna et al., 1995). In addition, a gain-of-function
mutant of \emph{CHS3} , encoding a TIR-NB-LRR-LIM protein, resulted in
enhanced defense responses and chilling-sensitivity phenotype (Yang et
al., 2010). Other studies found that some \emph{R} genes played positive
roles in plants tolerance to low temperature. A point missense mutation
of \emph{CHS1} , encoding a TIR-NB protein, resulted in a
chilling-sensitive phenotype (Wang et al., 2013; Zbierzak et al.,
2013).\emph{R} genes were also regulated by genes that contributed to
cold tolerance. In tobacco, overexpression of \emph{SgRVE6}
(\emph{REVEILLE} ) improved plant tolerance to low temperature and
up-regulated the expression of genes encoding NB-LRR proteins (Chen et
al., 2020).

Although many studies have investigated the molecular mechanisms
underlying low-temperature resistance, the roles of miRNAs in
low-temperature resistance in \emph{Brassica} crops are largely unknown.
In this study, we aimed to identify cold-responsive miRNAs in winter or
semi-winter Brassica and applied them to improving cold resistance of
Spring type Brassica crops. To do this, we constructed and sequenced
small RNA libraries from \emph{B. rapa} under cold stress and performed
a genome-wide association study (GWAS) related to the phenotype of cold
resistance using a \emph{B. napus} population. Combing the results from
both small RNA sequencing and GWAS, we found that\emph{Brassica}
-specific miR1885 was responsive to cold stress in both\emph{B. rapa}
and \emph{B. napus} . Previous studies have revealed that heat, as well
as \emph{Turnip mosaic virus} (\emph{TuMV} ) and\emph{Plasmodiophora
brassicae} infection, induced miR1885 accumulation (He et al., 2008; Yu
et al., 2011, Paul et al., 2021). Here, we further demonstrated that
miR1885 was also involved in the response to low-temperature stress
through its target genes and functioned as a negative regulator of
\emph{B. napus} in low temperature tolerance, providing potential
genetic candidates for cold-resistant rapeseed breeding.

\textbf{RESULTS}

\textbf{Genome-wide analysis of small RNAs responsive to cold stress
in\emph{B. rapa}}

To study the cold response of \emph{B. rapa} , the degree of cold
tolerance was estimated by evaluating morphological traits under cold
stress. After 7 days of cold treatment (4 degC), the plants stopped
growing and their leaves drooped (Figure 1A). Therefore, to identify
genome-wide cold-responsive miRNAs, we performed cold treatment for 7
days on 3-week-old seedlings of \emph{B. rapa} in a growth chamber, and
plants grown at 22 degC was used as control for small RNA sequencing.

Totally, 17.96 million and 17.17 million raw reads were obtained from
these two small RNA sequencing samples treated with 22 degC and 4 degC,
respectively (Figure 1B). Among them, most of small RNA sequences were
shared between these two libraries, while 2.82 million small RNA
sequences were specific in the seedling under normal temperature, and
2.23 million small RNA sequences were specific in the seedling under 4
degC (Figure 1B). Small RNA length distribution analysis showed that 24-nt
small RNAs were most abundant, and 21-nt small RNAs were second most
abundant (Figure 1C).

To learn the genomic loci that produced these small RNAs, 11.46 million
(22 degC) and 11.36 million (4 degC) clean reads were mapped to the \emph{B.
rapa} reference genome (http://brassicadb.cn/\#/). We found that the
small RNAs were generated from exons, introns, miRNAs, rRNAs, snRNAs,
snoRNAs, tRNAs, and repeat regions (Table S1). Notably, focusing on
small RNAs from non-coding gene regions, most abundant small RNAs were
produced from ribosomal RNAs (rRNAs) (22 degC, 44.67\% and 4 degC, 48.10\%),
unknown RNAs (22 degC, 27.65\% and 4 degC, 25.22\%), and miRNA precursors
(22 degC, 8.50\% and 4 degC, 8.71\%) (Figure S1).

\selectlanguage{english}\textbf{Identification of \emph{B. rapa} miRNAs responsive to cold
stress}

To identify known miRNAs in our sequence libraries, we aligned the small
RNA sequences with mature miRNAs in miRBase and the \emph{B.
rapa}reference genome, allowing no more than two mismatches. In total,
77 mature miRNAs related to 184 conserved miRNA precursors were found in
our libraries. Among them, 52 miRNAs were conserved with those
in\emph{Arabidopsis} , while 16 miRNAs were \emph{Brassica} -specific
miRNAs (Table S2). To further identify novel miRNA candidates, the
flanking sequences of unknown small RNAs in genome were extracted for
RNA secondary structure prediction using RNAfold. Based on the plant
miRNA criteria (Axtell et al., 2018, see the detail in the method), we
identified 11 highly credible novel miRNAs (Table S3). Among them,
novel-miR1 and novel-miR11 were identified with more than 100 reads.

Among the 77 known miRNAs, 16 miRNAs were differentially expressed
between the sample under cold treatment and the control
(\textgreater{}1.5-fold) (Table 1). Among these mature miRNAs,
bra-miR157a\textasciitilde{}c, bra-miR161, bra-miR393f, bra-miR396b,
bra-miR319c, bra-miR156k, bra-miR1140, bra-miR158a and
bra-miR156a\textasciitilde{}f were down-regulated; and bra-miR168d,
bra-miR403, bra-miR390a\textasciitilde{}b,
bra-miR398b\textasciitilde{}c, bra-miR408, and bra-miR1885 were
up-regulated under cold treatment. Among them, the\emph{Brassica-}
specific bra-miR1885 was induced more than four-fold by cold stress.
Among the novel miRNAs, only novel-miR11 was cold-responsive. Taken
together, the small RNA profiling provided potential candidates of
cold-responsive miRNAs in \emph{B. rapa} .

\textbf{Validation of cold-responsive mature miRNAs in \emph{B. rapa}}

To validate the cold-responsive miRNAs detected by small RNA sequencing,
we further analyzed the cold responsiveness of mature bra-miR1885,
bra-miR398, bra-miR408, bra-miR161, bra-miR157, bra-miR393, and
bra-miR168 using qRT-PCR. We found that the abundance of bra-miR1885,
bra-miR398 and bra-miR168 were significantly up-regulated by cold
stress, whereas bra-miR161, bra-miR157 and bra-miR393 were significantly
down-regulated under cold stress, as compared to normal temperature at
the same growth stage. These results were consistent with the small RNA
sequencing data (Figure 2A). Among these miRNAs, bra-miR398 and
bra-miR168 were conserved miRNAs that were induced more than three-fold
by cold treatment (Figure 2A). We then performed northern blotting
analyses to further con\selectlanguage{english}firm the cold responsiveness of bra-miR398,
bra-miR168 and bra-miR1140. Consistently, the accumulation of bra-miR398
and bra-miR168 was much higher in 4 degC than in 22 degC, and\emph{Brassica}
-specific miR1140 was repressed after cold treatment (Figure 2B). In
addition, the results of qRT-PCR analyses further confirmed the
abundance of mature bra-miR1885 was increased in leaf under cold stress
as compared to normal condition (Figure 2C).

Next, to identify the potential cold-regulated \emph{cis} -elements in
the \emph{bra-MIR1885} promoter, we extracted the promoter sequence
of\emph{bra-MIR1885} using Promoter 2.0 Prediction Server
(https://services.healthtech.dtu.dk/service.php?Promoter-2.0), and
analyzed the \emph{cis-} elements in the putative
\emph{bra-MIR1885}promoter sequence using PlantCare
(http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). These
analyses revealed that the promoter contained light-responsive,
defense-responsive, stress-responsive, and low temperature-responsive
elements (LTR: CCGAAA) (Figure S2). A previous study has revealed that
LTR motif (CCGAAA) within the promoter region of \emph{CsINV5} was the
core \emph{cis} element in response to low temperature (Qian et al.,
2018). This result suggested that the \emph{bra-MIR1885} promoter
contained a low temperature-responsive element that potentially
regulated the expression of \emph{bra-MIR1885} under cold stress.

\textbf{Identification of low temperature resistance-associated loci
containing MIRNA genes in \emph{B. napus} using Genome-wide association
study}

\emph{B. napus} is closely related to \emph{B. rapa} and both
are\emph{Brassica} crops. To identify MIRNA genes associated with low
temperature stress in \emph{Brassica} crops, we constructed a \emph{B.
napus} population with 220 accessions. The phenotype of low-temperature
resistance was divided into 5 levels (Grade 1 indicates most
cold-resistant, while Grade 5 indicates most cold-sensitive) and was
measured in the field in January of 2017 and 2018, and the phenotypic
data of rapeseed population showed a nearly skewed distribution (Xu et
al., 2021). The genome re-sequencing of these 220 \emph{B.
napus}accessions generated 2.13 Tb clean reads, and single nucleotide
polymorphism (SNP) genotyping of \emph{B. napus} population identified
3.80 million highly quality SNPs (MAF \textgreater{} 0.05, geno
\textless{} 0.2) (Xu et al., 2021). Overall, the SNPs were not evenly
distributed across the whole genome. The SNP density in the A sub-genome
(5.25 SNP/Kb) was higher than that in the C sub-genome (4.74 SNP/Kb)
(Figure S3A). A principal component analysis (PCA) was used to assess
the genetic relationship in the association panel. Based on the PCA,
the\emph{B. napus} population was roughly divided into three groups
corresponding to winter (W), semi-winter (SW), and spring (S) ecotypes
(Figure S3B). We further analyzed the linkage disequilibrium (LD)
throughout the genome. The LD was estimated
as\emph{r\textsuperscript{2}} (the squared Pearson correlation
coef\selectlanguage{english}ficient) between all pairs of SNP markers. The average distance over
which LD decayed to \emph{r\textsuperscript{2}} =0.2 of its maximum
value was 100 kb, indicating that SNPs located in the 200-kb genomic
region around each peak SNP represented a GWAS-quantitative trait locus
(QTL) (Figure S3C).

The GWAS was performed using a mixed linear model, and the SNPs
significantly associated with low-temperature resistance were identified
at a threshold of \emph{p} \textless{} 10\textsuperscript{-6}. As a
result, we identified 32 highly credible GWAS-QTLs (Figure 3A). To
discover candidate MIRNA genes within these QTLs, we collected 163
pre-miRNA sequences specific to the Cruciferae family from four species
(\emph{B. napus} , \emph{B. rapa} , \emph{B. oleracea,}
and\emph{Arabidopsis} ) according to previous reports (Jones-Rhoades and
Bartel, 2004; Alves-Junior et al., 2009; Yu et al., 2011; He et al.,
2018; Zhang et al., 2018; Li et al., 2019). All the pre-miRNA sequences
were mapped to the \emph{B. napus} reference to identify their location
(Table S4). We found that four candidate MIRNA genes,\emph{bna-MIR166C}
, \emph{bna-MIR1885} , \emph{bna-MIR168A} and\emph{bna-MIR845A,} were in
the QTLs associated with low-temperature resistance.

A previous study proved that knock-down of miR165/166 expression
conferred a cold-resistant phenotype in \emph{A. thaliana} (Yan et al.,
2016). Conserved miR845a was known to be differentially expressed in
leaves of winter turnip rape under cold stress (Zeng et al., 2018). To
further determine whether these four miRNAs were responsive to low
temperature, we measured the abundance of mature miRNAs in \emph{B.
napus} under cold stress by qRT-PCR. From this analysis, we found that
all of them were cold-responsive. bna-miR1885 and bna-miR168a were
up-regulated under cold stress compared to under normal condition at the
same growth time points, while bna-miR845a and bna-miR166c were
down-regulated (Figure 3B and Figure 4A).

To identify if any genetic variation in cis-regulatory element
contribute to the accumulation of mature miRNAs in the population of
rapeseed, we studied the LD of SNPs in \emph{bna-MIR1885, bna-MIR166C,
bna-MIR168A and bna-MIR845A} gene locus, respectively. We found that the
SNPs in their corresponding genes exhibited significant linkage (Figure
S3D, F, G, H). Importantly, we identified 22 SNPs and 7 INDELs
in\emph{bna-MIR1885} promoter region, which were found to be significant
associated with cold resistance. We further investigated the five
variants located in the 200-bp genomic region around the low
temperature-responsive (LTR: CCGAAA) element identified above.
The\emph{B. napus} accessions were classified into three haplotype
groups based on the genotypes of these variants (Figure 3C). We found
that 179 accessions contained the H1 haplotype, which variants were
identical with the genome reference, and this haplotype was associated
with low temperature-sensitive phenotype. 24 accessions carried H2
haplotype, which contained mutated LTR cis-element and was associated
with low temperature-resistant phenotype. Consistently, H3 haplotype had
an average resistance degree, which indicated more resistant than H1
haplotype, and less resistant than H2. Next, we randomly selected 4
accessions with low temperature-resistant variant alleles and low
temperature-sensitive variant alleles, respectively, to measure their
abundance of miR1885 after cold treatment by qRT-PCR. The qRT-PCR data
indicated that the abundance of miR1885 was significantly higher in low
temperature-sensitive accessions with normal LTR in the promoter (H1
haplotype), as compared to low temperature-resistant accessions with
mutated LTR in their promoter (H2 haplotype) (Figure S3E). In total,
these results suggested that genetic variation in the LTR region
at\emph{MIR1885} promoter in the rapeseed population potentially
affected the miR1885 expression and were associated with cold
resistance.

\textbf{Identification of miR1885 target genes and their cold
responsivity in \emph{B. napus}}

Because bra-miR1885 and bna-miR1885 shared the same mature sequence
(Figure 4B), we used miR1885 to represent the two miRNAs
in\emph{Brassica} crops hereafter. We \selectlanguage{english}first predicted the candidate
target genes of miR1885 in \emph{B. napus} using a plant small RNA
target analysis tool (psRNATarget). Two \emph{R} genes,
\emph{Bn.TIR.A09}(BnaA09g14980D) and \emph{Bn.TNL.A03} (BnaA03g56180D),
were identified as highly credible targets, and the miR1885 target sites
were within the TIR domain. To verify the miR1885 cleavages of the
predicted target genes, we performed a 5\selectlanguage{english}' rapid amplification of cDNA
ends (5\selectlanguage{english}'-RACE) assay in \emph{B. napus} and detected miR1885 cleavage
sites frequently within \emph{Bn.TIR.A09} and \emph{Bn.TNL.A03} (Figure
4C and Figure S4), indicating that they were targeted by miR1885 in
\emph{B. napus} .

Next, we examined the transcript levels of these two target genes under
cold stress, using three low temperature-responsive marker genes
(\emph{Bn.CBF1.C03} , \emph{Bn.HSP70.A01} and \emph{Bn.ERF105.C09} ) as
positive controls (Wang and Hua, 2009; Bolt et al., 2017; Liu et al.,
2018). The target genes were down-regulated under cold stress in\emph{B.
napus} (Figure 4D), while the expression of \emph{Bn.CBF1.C03}and
\emph{Bn.HSP70.A01} were induced and \emph{Bn.ERF105.C09} were repressed
in samples under cold treatment as compared to those under normal
condition, consistent with those reported in previous studies (Figure
4E).

\textbf{Overexpression of miR1885 in \emph{B. napus} led to
down-regulation of its target genes \emph{Bn.TIR.A09}
and\emph{Bn.TNL.A03}}

Genomic synteny analysis using TBtools
(https://github.com/CJ-Chen/TBtools) suggested that bna-miR1885
in\emph{B. napus} is inherited from the genome of \emph{B. rapa} through
genome merging (Figure S5), and bna-miR1885 and bra-miR1885 shared the
same mature sequence. To further verify the function of miR1885, we
constructed a vector with the \emph{bra-MIR1885} gene under the control
of the AA6 promoter (Patent WO 2007/069894) and transformed it
into\emph{B. napus} cv. K407 (Figure 5A). The strategy using \emph{B.
napus}system for functional analysis of \emph{B. rapa} genes have been
used in several studies (Liu et a li., 2014; Tian et al., 2018; Cui et
al., 2020). In total, seven positive transgenic lines
(T\textsubscript{0}) were screened by PCR (Figure 5B and C). The
expression level of pri-miR1885 was con\selectlanguage{english}firmed by qRT-PCR. We found that
pri-miR1885 was overexpressed in the leaves of X2, X8, X11, X12, and X15
lines (Figure 5D). Compared with WT, there were no obvious morphological
abnormalities among miR1885-OE lines at seeding stage (Figure 5B). We
then further checked the abundance of mature miR1885 in the transformed
lines by qRT-PCR and northern blotting analyses. As expected, mature
miR1885 was found to be over-accumulated in the X2, X8, X11, X12, and
X15 lines (Figure 5E and 5F). And then, we conducted qRT-PCR analyses to
validate the in\selectlanguage{english}fluence of the miR1885 on the expression level of the two
target genes in the transgenic lines. As expected, transcript levels
of\emph{Bn.TIR.A09} and \emph{Bn.TNL.A03} were lower in the transgenic
lines than in wild type (Figure 5G and 5H).

Next, we checked the RNA levels of miR1885 and its targets in miR1885-OE
lines and WT with or without cold treatment by qRT-PCR. As expected, the
abundance of miR1885 was increased in miR1885-OE lines and WT after cold
treatment (Figure 6A). In addition, the RNA levels of
\emph{Bn.TIR.A09}were further decreased in miR1885-OE under cold stress
(Figure 6B), but those of \emph{Bn.TNL.A03} were unchanged (Figure 6C),
suggesting that there was additional cold-triggered regulation of
\emph{Bn.TNL.A03} at the transcriptional level, while the cold response
of \emph{Bn.TIR.A09}was majorly contributed by miR1885-mediated
post-transcriptional regulation.

Given that CBF/CORs played critical roles in plant response to low
temperature, we investigated whether miR1885 can affect the cold
induction of CBF/COR factors. To do this, we analyzed the expression
levels of \emph{Bn.CBF1.CO3} , \emph{Bn.CBF2.Ann} and
\emph{Bn.COR15.A03}in the miR1885-OE under cold treatment. However,
there were no significant differences at the expression levels of
\emph{Bn.CBF1.CO3} ,\emph{Bn.CBF2.Ann} and \emph{Bn.COR15.A03} after
cold treatment between WT and miR1885-OE lines (Figure 6D-F). This
result suggested that miR1885 regulated low temperature tolerance in
\emph{Brassica} through CBF/CORs-independent pathway.

\textbf{miR1885 negatively regulated low temperature tolerance
in\emph{Brassica}}

A previous study reported that knockdown of \emph{CHS1} , encoding a
TIR-NB protein, led to a chilling-sensitive phenotype in \emph{A.
thaliana} (Zbierzak et al., 2013). We found that two TIR-NB transcripts
were the targets of miR1885 and both miR1885 and these two TIR-NB genes
were cold-responsive. Therefore, we further tested the effects of
miR1885 and its targets on low-temperature resistance in \emph{B. napus}
. Semi-winter type WT K407 and miR1885-OE plants were grown for 6 weeks
at 22 degC, and then were transferred to 4 degC for 30 days or 2 degC for 7
days. We found that the miR1885-OE plants were more sensitive as
compared to WT under 2 degC for 7 days (Figure 7A) but showed no apparent
phenotypic difference under 4 degC (Figure S6A). This result indicated
miR1885 functioned as negative regulator in rapeseed response to low
temperature. We further performed phenotypic studies using WT K407 and
miR1885-OE seeds during post-germination under normal and cold
environment. The seeds were harvested at the same season and the cold
environment was 6 degC. The hypocotyl length (NHL) of WT K407 and
miR1885-OE showed no difference after growing for 2 days on the
germination bed without cold treatment. However, after 7 days cold
treatment, the miR1885-OE exhibited shorter NHL compared with WT K407
(Figure S6B-D), indicating that miR1885 also played critical role in
cold response in the early developmental stage, when no significant
phenotypic difference between WT and miR1885-OE was observed under
normal condition.

To test that whether the cold-responsive miR1885 can be applied to
improve low temperature tolerance of Spring type \emph{Brassica} , the
miR1885 knockdown plants (\emph{STTM1885} ) generated in Sping
type\emph{Brassica napus cv.} `Westar' were used for further study (Cui
et al., 2020). The qRT-PCR results showed that the accumulation of
miR1885 was reduced in \emph{STTM1885} plants. Compared with those in
wildtype Westar, the expression levels of \emph{Bn.TNL.A03} and
\emph{Bn.TIR.A09}were increased in \emph{STTM1885} plants (Figure 7B).
To study morphological responses to cold stress, we then moved the
Westar and\emph{STTM1885} plants into 2 degC environment. After 14 days
cold treatment, a cold-tolerant phenotype was observed in
\emph{STTM1885}plants (Figure 7C). This finding further suggested that
attenuating miR1885 resulted in enhanced cold tolerance of rapeseed
potentially via their targets.

We further tested whether overexpression of miR1885 affected plant
growth at freezing temperature compared to its background winter ecotype
K407, while the spring ecotype Westar was used as negative control
(Figure S7A). The abundance of mature miR1885 in wild-type \emph{B.
napus} under freezing stress was examined by real-time PCR. Based on
previous studies, we set the freezing temperature at \selectlanguage{english}-10 degC (Xu and Cai,
2019; Ljubej et al., 2021). We found that miR1885 was up-regulated under
freezing stress, and its expression level was approximately 2.4-fold
than in wild type at 2.5 h (Figure S6B). The wild type and miR1885-OE
lines were subjected to a \selectlanguage{english}-10 degC treatment for 2.5 h, and then allowed
to recover at 22 degC for 7 days. The survival rate of wild type K407 was
78.00\%, and those of miR1885-OE were much lower, ranging from 22.33\%
to 55.33\%, whereas 100\% of spring ecotype Westar died (Figure S7A and
Figure S7C). Accumulation of free proline (Pro) during cold or freezing
stress is thought to protect plants (Xin and Browse, 1998). To test
whether miR1885-OE had altered Pro content, we determined the Pro
content in wild type and miR1885-OE with or without freezing treatment
(Figure S7D). As we expected, the Pro content was higher in wild type
and miR1885-OE grown at \selectlanguage{english}-10 degC than in those grown at 22 degC. In plants
grown at \selectlanguage{english}-10 degC, more Pro accumulated in the wild type plants than in
the miR1885-OE. This \selectlanguage{english}finding was consistent with the freezing-sensitive
phenotype of transgenic lines subjected to a \selectlanguage{english}-10 degC treatment. These
results demonstrated that overexpression of miR1885 increases plant
sensitivity to freezing at the seedling stage.

\textbf{DISCUSSION}

At present, many MIRNA genes functioning in stress response have been
identified by small RNA profiling and RNA-seq (Ahmed et al., 2020; Ma et
al., 2022; Zhu et al., 2022). However, the MIRNA genes that function in
the response to cold stress in crops remain largely unexplored. In this
study, we used small RNA profiling and GWAS to identify cold-responsive
candidate miRNAs or MIRNA genes associated with low-temperature
resistance. The results from these two methods provide mutual support
for the roles of miR1885 in the cold response of \emph{Brassica} crops.

We identified up-regulated or down-regulated miRNAs in small RNA
sequencing libraries of \emph{B. rapa} and candidate MIRNA genes from
GWAS loci associated with low-temperature resistance in \emph{B. napus}
. The small RNA profiling suggested that miR1885 was drastically
triggered by cold stress. Using GWAS, we found \emph{MIR1885} gene was
in the loci associated with low-temperature resistance. We further
confirmed that miR1885 was induced by cold stress in both \emph{B. rapa}
and \emph{B. napus} . We used the newest version of \emph{B. rapa} and
\emph{B. napus}genome reference to analyze their genomic synteny. As
expected, the\emph{B. napus} miR1885 (ChrA06) is inherited from the
genome of\emph{B. rapa} (ChrA06 genome) through genome merging (Fig.
S5). Additionally, For the cold responsive miRNAs in this study, we
found some miRNAs have been reported to regulate flowering, such as
miR156/7, miR168 and miR158. This suggested that some the
cold-responsive miRNAs also play roles in vernalization.

Mature miR1885 is a \emph{Brassica} -specific miRNA that targets
the\emph{R} gene family. In \emph{A. thaliana} , the TIR-NB-LRR protein
RPP4 and CHS1 were involved in plant resistance to low temperature
stress (Huang et al., 2010; Zbierzak et al., 2013). Here, two \emph{R}
genes,\emph{Bn.TIR.A09} and \emph{Bn.TNL.A03} , were predicted as
targets for cleavage by miR1885 in \emph{B. napus} . These \emph{R}
genes belonged to two different subgroups: \emph{Bn.TIR.A09} encoded a
protein with the TIR motif lacking NBS-LRR domain; while
\emph{Bn.TNL.A03} encoded a protein with the entire TIR-NBS-LRR
structure. We found that these two targeted genes were down-regulated in
wild-type \emph{B. napus} under cold treatment, opposite to the
expression pattern of miR1885 (Figure 6). These results demonstrated
that miR1885 negatively regulated these\emph{R} genes under cold
temperatures.

Overexpression of miR1885 in semi-winter type \emph{B. napus} attenuated
rapeseed sensitivity to low temperature at the seedling stage. While
knocking down of miR1885 in Spring type \emph{B. napus}
(\emph{STTM1885}lines) improved plant tolerance to low temperature
(Figure 7). Additionally, we did not find obvious differences between
knock-down (STTM1885) lines and control lines at the early seeding
stage. But\emph{STTM1885} did delayed the flowering time and increased
the plant height at a later stage. The RNA levels of miR1885-cleaved R
genes\emph{Bn.TIR.A09} and \emph{Bn.TNL.A03,} which increased
in\emph{STTM1885} rapeseed and decreased in miR1885-OE lines, may
contribute to rapeseed response to cold stress. Although a link between
miR1885 and cold response has been revealed, the hiding mechanism is
still uncovered. In addition, we have checked some CBF-dependent genes
in miR1885-OE lines, including \emph{CBFs} and \emph{CORs} , and found
that the transcript levels of these genes in miR1885-OE in response to
cold stress were no different with that of wild type, suggesting that
miR1885-targets regulated plant resistance to low temperature
independent on CBFs pathway.

In this study, we revealed a link between miR1885 and cold response in
Brassica through its targeted R genes. However, how miR1885 and these R
genes involve in CBF-independent cold response in rapeseed requires
further examination. In Arabidopsis, 22-nt miRNAs can trigger the
production of 21-nt phased secondary siRNAs (Cuperus et al., 2010; Li et
al., 2012). 22-nt miR482/2118, which widely present in seed plants,
triggered phasiRNA biogenesis from its targeted R genes (NBS-LRR),
forming another layer of regulation to reinforce its silencing effect
(Fei et al., 2013; Liu et al.,2020). Therefore, 22-nt miR1885 in
rapeseed possibly triggers phasiRNAs biogenesis from its targeted R
genes, and these phasiRNAs may target unknown downstream genes to
regulate plant response to low temperature stress. Additionally,
co-immunoprecipitation (co-IP) followed by~mass spectrometry can be
performed to identify proteins interacting with miR1885-regulated R
proteins under cold stress, which may play role in cold signal transfer
and resistance to cold stress.

In rice, OsmiR1320 was repressed by cold stress. However, overexpressing
OsmiR1320 enhanced plant tolerance to cold stress, indicating that
OsmiR1320 played positive role in cold response (Sun et al., 2022). On
the contrary, we found that, as a negative regulator of cold response,
miR1885 was induced under cold stress potentially via its
low-temperature response (LTR) cis-element in \emph{MIR1885} promoter.
Therefore, \emph{MIRNA} transcription could be regulated incoherently by
cold stress. The genetic mutation in and close to this LTR in the
rapeseed population may contribute to their diversity of cold response
(Figure 3C and Figure S3E). Genetic editing in this LTR could be used
for detaching the cold response of \emph{MIR1885} , and therefore
enhance rapeseed resistance to cold stress. Previous works showed that
over-expression of miR1885 affected development traits and biotic
response, including flowering time and disease resistance. The abundance
of miR1885 may regulated by different signals through potentially
different promoter cis elements and played a very important role in
regulating both plant growth and development.

In summary, using multi-omics data, we identified miR1885 involved in
plant response to low temperature stress by small RNA sequencing and
GWAS. Overexpression of miR1885 in \emph{B. napus} not only increased
plant sensitivity to low temperature, while knockdown of miR1885
improved plant tolerance to low temperature. Our findings suggested
that\emph{MIR1885} and its target genes can be used as genetic resource
for improving resistance to low temperature in the breeding
of\emph{Brassica} crops.

\textbf{METHODS AND MATERIALS}

\textbf{Plants materials, cold treatment, and freezing tolerance assay}

The wild-type \emph{B. rapa} variety `Bre' was grown under a 16-h
light/8-h dark (22 degC) photoperiod for 3 weeks. To analyze the miRNAs
abundance in `Bre' under cold stress, some `Bre' plants were moved to 4
degC condition for 7 days cold treatment and some were grown at 22 degC for
additional 7 days as control. The aboveground parts of treated and
untreated were then harvested, frozen in liquid nitrogen, and stored at
\selectlanguage{english}-80 degC until extraction of total RNA. Small RNA sequencing (one
biological replicate) was performed using the Illumina GAII sequencer at
BIG (ShenZhen, China).

The wild-type \emph{B.} \emph{napus} variety `K407', a Semi-winter
ecotype, was grown~under a 16-h light/8-h dark (22 degC) photoperiod with
light intensity of 250\selectlanguage{english} \selectlanguage{greek}μ\selectlanguage{english}mol m\textsuperscript{-2}s\textsuperscript{-1}.
For the cold resistance assay, the plants were grown under standard
conditions (22 \selectlanguage{ngerman}°C) for 6 weeks before transferring to a cold environment
(4 °C or 2 °C, 16-h/8-h light/dark photoperiod). The similar methods for
studying cold stress in Brassica crops have been used in previous
reports (Raza et al., 2021; Hussain et al., 2022). After cold treatment,
the cold resistance of plants was determined. The freezing tolerance
assay was performed in a freezing chamber, after cold acclimation for 24
h, in which wild-type and transgenic lines were subjected to a \selectlanguage{english}-10 degC
treatment for 2.5 h. The plants were then allowed to recover at 22 degC
for 7 days before determining the survival rate. The Pro content was
measured immediately after the freezing treatment using a commercial
assay kit (A107-1-1, Nanjing JianCheng Bioengineering Institute,
Nanjing, China), according to the manufacturer's instructions.

The wild-type \emph{B. rapa} cultivar ``Bre'', \emph{B. napus} cultivar
`K407' are all inbred lines.

\textbf{Identification of known and novel miRNAs in \emph{B. rapa}}

We aligned the clean reads with known miRNAs in the miRBase
(https://www.mirbase.org/), allowing two mismatches using TopHat
(http://ccb.jhu.edu/software/tophat/index.shtml). The potential miRNAs
were mapped to \emph{B. rapa} reference genome
(http://brassicadb.cn/\#/) to extract the 600-bp flanking sequence. The
flanking sequence of miRNAs were also aligned with known pre-miRNAs to
confirm their identities of miRNA family. The flanking sequence of
unknown small RNAs was also used to predict secondary structures using
RNAfold software
(http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). To
identify novel miRNA, Small RNAs that showed less than 2 mismatches with
known miRNAs in miRbase were filtered. The potential miRNAs that had
fewer than 10 hits in the reference genome were retained. 600-bp
flanking sequences of small RNAs were used to predict their secondary
structures using RNAfold software with default parameters. Novel miRNAs
were identified according to the plant miRNA criteria (Axtell et al.,
2018, Yu et al., 2011). In brief, a key judgement was made by
determining whether the potential miRNAs and their miRNA* were on the
stem with fewer than four bulges, and the miRNA candidate were most
abundant among small RNAs generated from given miRNA precursors. After
Blast these novel miRNA precursors in NCBI, we found them were only
detected in Brassica, therefore, they were defined as Brassica-specific
miRNAs.

\textbf{Real time PCR}

Total RNA was isolated with TRIzol reagent (Invitrogen, Carlsbad, CA,
USA) according to the manufacturer's instruction. High-quality total RNA
(A260/280 \textgreater{} 2.00) was reverse-transcribed into cDNA. Small
RNA qRT-PCR was conducted using a miRNA First Strand cDNA Synthesis kit
(stem-loop method) (Sangon Biotech, Shanghai, China). A total of 3 \selectlanguage{greek}μ\selectlanguage{english}g
total RNA was used for first-strand cDNA synthesis with the stem-loop
primer and \emph{U6} reverse primer in a 20-\selectlanguage{greek}μ\selectlanguage{english}L reaction volume.
The\emph{U6} snoRNA was used as an internal control. For qRT-PCR
analysis of coding genes, 2 \selectlanguage{greek}µ\selectlanguage{english}g total RNA treated with DNase I was
reserve-transcribed into cDNA using M-MuLV (Takara, Otsu, Japan) with
oligos (dT). The internal control was \emph{ACT7} . qRT-PCR was
performed using SYBR mix (YEASEN, China) on a Bio-CFX96 instrument.
Three biological replicates and three technological replicates were
analyzed. Real-time PCR and reverse transcriptase (RT)-PCR were
performed using gene-specific primer pairs (Table S5).

\textbf{Northern blot analyses}

A 30--45 \selectlanguage{greek}µ\selectlanguage{english}g aliquot of high-quality total RNA (A260/280
\textgreater{}2.00) was separated on a 19\% polyacrylamide denaturing
gel, and then transferred to a Hybond membrane (GE healthcare) for 2 h
at 200 mA at 4 \selectlanguage{ngerman}°C. In brief, after crosslinking for 5 min with
ultraviolet irradiation, the membrane was hybridized overnight at 45 °C
with DNA probes (Yu et al., 2011). Washing and autoradiography of the
Hybond membrane were performed according to the instructions of the
North2South\selectlanguage{english} Chemiluminescent Hybridization and Detection kit (No.17097,
Pierce, Rockford, IL).

\textbf{Genome-wide association analysis (GWAS)}

A diversity panel, consisting of 220 \emph{B. napus} inbred lines was
used for the GWAS in this study. The accessions were grown using a
randomized complete block design with two replicates. Each accession was
grown in plots with four rows. For the phenotypic data about low
temperature resistance of rapeseed population showed a nearly skewed
distribution. For more details about field experiment design and the
distribution of phenotype were described in our previous reports (Xu et
al., 2021). More specifically, the population included winter,
semi-winter and spring oilseed rape. The accessions originated from
Asia, Europe, and North America. The \emph{B. napus} population was
grown at the Farm Station in Yangling, Shaanxi province, China
(107deg59\selectlanguage{english}'--108deg08\selectlanguage{english}'E; 34deg14\selectlanguage{english}'--34deg20\selectlanguage{english}'N), where the local temperature ranged
from \selectlanguage{english}-15 to 10 degC in winter
(\emph{https://lishi.tianqi.com/xianyang/202101.html}). The reference
genome of \emph{B. napus} ``Darmor'' (version 4.1) was used to call
variants, which was downloaded from BRAD (http://brassicadb.cn/\#/)
(Chalhoub et al., 2014). Low-temperature resistance was measured in the
field in January of 2017 and 2018. The grade of low-temperature
resistance was evaluated with criteria including grade 1 (the plants are
normal, and the leaves are green), grade 2 (a small part of leaves are
necrotic), grade 3 (a half of leaves are necrotic), grade 4 (most of
leaves are necrotic), and grade 5 (all leaves are necrotic). The
association analysis for low-temperature resistance was carried out
using a MLM with the software package TASSEL 5
(https://www.softpedia.com/get/Science-CAD/TASSEL.shtml). The population
structure was calculated using PCA. The software PopLDdecay
(https://github.com/BGI-shenzhen/PopLDdecay) was used to calculate LD.
The parameter \emph{r\textsuperscript{2}} =0.2 was read as LD length.
The threshold of the GWAS was set to \emph{P} -value \textless{}
10\emph{\textsuperscript{--6}} .

\textbf{5'-RACE}

The 5'modified RACE was performed using the FirstChoice RLM-RACE Kit
(Invitrogen). In brief, total RNA was isolated from young leaves
of\emph{B. napus} . Then, 8--10 \selectlanguage{greek}µ\selectlanguage{english}g high-quality total RNA
(A260/280[?]2.00) was directly ligated to the 5'-RACE adapter without calf
intestinal phosphatase and tobacco acid pyrophosphatase treatment. The
first-strand cDNA was synthesized using M-MLV reverse transcriptase
according to the manufacturer's instructions. The PCR amplification was
performed using the 5' outer primer and gene-specific outer primer, and
0.1 \selectlanguage{greek}µ\selectlanguage{english}L PCR product was used as the template for nested PCR. Nested PCR
amplification was performed using the 5' inner primer and gene-specific
inner primer. The PCR products were gel-purified and then cloned and
sequenced.

\textbf{Plasmid construction, genetic transformation, and phenotypic
analysis}

The plant overexpression vector was generated using pCAMBIA1301
harboring the AA6 promoter and the tAA6 terminator.
The\emph{bra-MIR1885} genomic DNA fragment was isolated from `Bre'
seedlings and cloned into the pCAMBIA1301 binary vector. The binary
construct was introduced into \emph{Agrobacterium tumefaciens} strain
GV3101. Genetic transformation of the \emph{bra-MIR1885} gene construct
into \emph{B. napus} was performed as described elsewhere (Bhalla and
Singh, 2008). Seedlings exhibiting resistance to hygromycin were
transplanted and grown in a greenhouse at 22 \selectlanguage{ngerman}°C under a 16-h/8-h
light/dark photoperiod. The genomic DNA extracted from
T\textsubscript{0} plants was used for PCR amplification with specific
primers. The seeds (T\textsubscript{1}) were harvested separately from
positive plants for further analyses.

The wild-type and transgenic lines (T2\textasciitilde{}T3 generations
homozygote lines) were grown in the field, and the genomic DNA extracted
from plants was used for PCR amplification with specific primers. Each
transgenic line was grown in a row of 10 plants, with 30-cm spacing
between plants and 40 cm spacing between rows. The phenotypes of
transgenic and wild-type plants were measured in the field.

\textbf{AUTHOR CONTRIBUTIONS}

P.X., W.Z., X.Y., and Y.H. conceived the project; P.X., W.Z. and Y.H.
designed the experiments; P.X., X.W. and Y.Z. performed the experiments;
P.X. and W.Z. analyzed and interpreted the data; P.X., W.L., X.Y., and
Y.H. wrote the manuscript.

\textbf{ACKNOWLEDGEMENTS}

We would thank Dr. Hui-Shan Guo, the Institute of Microbiology, Chinese
Academy of Sciences for providing the STTM1885 rapeseed material. This
work was supported by China Postdoctoral Science Foundation
(2022M72211), Shanghai Pujiang Program (21PJ1407600), National Natural
Science Foundation of China (Grant No. 32170581, 31471883, 31571261),
National Programs for Science and Technology Development of China (Grant
No. 2016YFD0101900), and Science and Technology Program of Shaanxi
Province (2019NY-038).

\textbf{CONFLICT OF INTEREST}

The authors declare no conflict of interest.

\textbf{DATA AVAILABILITY STATEMENT}

All the data and materials that support the findings of this study are
available upon request from the corresponding author. The sequencing
data generated in this study have been deposited in the National Center
for Biotechnology Information Sequence Read Archive database under the
accession PRJNA748869 and PRJNA745269.

\textbf{TABLES}

\textbf{Table 1. The miRNAs that were down- or up-regulated by more than
1.5-fold change under cold treatment (} 4 °C\textbf{) as compared to
normal temperature (} 22 °C\textbf{).}

RPM: Reads per million.

\textbf{FIGURE LEGENDS}

\textbf{Figure 1. Phenotypes of \emph{B. rapa} after cold treatment and
global view of small RNA reads in sequencing libraries.}

\textbf{(A)} Phenotypes of wild-type \emph{B. rapa} plants grown in soil
at 22 °C (left) and 4 °C (right) for 7 days. Red scale bar represents 5
cm. \textbf{(B)} Number of total reads in small RNA library of \emph{B.
rapa} . Shared small RNAs were presented in both libraries. Specific
small RNAs were those detected only in normal temperature \textbf{(} 22
°C\textbf{)} or low temperature \textbf{(} 4 °C\textbf{)}
library.\textbf{(C)} Length distribution of unique and raw reads in the
two small RNA libraries.

\textbf{Figure 2. Abundance of differentially expressed miRNAs between
cold-stress and normal conditions.}

\textbf{(A)} qRT-PCR analysis of seven differentially expressed miRNAs
in libraries with normal temperature and low temperature.
\textbf{(B)}Relative abundance of three differentially expressed miRNAs
under normal condition and cold condition detected by northern blotting.
\textbf{(C)}Relative abundance of miR1885 of \emph{B. rapa} plants under
22 °C and 4 °C condition for 0, 7, 14, 21 and 28 days, respectively.
SAM: shoot apical meristem. ``D'' represents days. Error bar represents
mean ± SE;\emph{n} = 3. \emph{U6} was used as the endogenous reference
gene in expression analyses. * indicates \emph{P} \textless{} 0.05; **
indicates\emph{P} \textless{} 0.01; Student's \emph{t} -test.

\textbf{Figure 3. Genome-wide association study of low-temperature
resistance in a panel of \emph{B. napus} accessions.}

\textbf{(A)} Manhattan plots and Q-Q plots for genome-wide association
study of low temperature resistance phenotype measured in Year 2017
(top) and Year 2018 (bottom), respectively. \textbf{(B)} qRT-PCR
analysis of the abundance of three candidate miRNAs under 22 °C and 4 °C
condition for 0, 10, 20 and 30 days, respectively. Error bar represents
mean ± SE; \emph{n} = 3. \emph{U6} was used as the endogenous reference
gene. * indicates \emph{P} \textless{} 0.05; ** indicates
\emph{P}\textless{} 0.01; Student's t-test. (\textbf{C} ) Haplotypes of
candidate SNPs within \emph{MIR1885} promoter for low-temperature
resistance in\emph{B. napus} population. LTR represents low-temperature
element (CCGAAA). Five variants, including four SNPs and the one-bp
deletion, which were significantly associated with low-temperature
resistance, were marked with dark bold vertical lines. The five variants
located in the 200-bp genomic region around LTR were used for analysis.
`Hap.' denotes the haplotypes of \emph{MIR1885} ; `Number' denotes the
number of lines belonging to each haplotype group. `Degree' denotes the
average grade of low-temperature resistance.

\textbf{Figure 4. qRT-PCR analysis of the abundance of miR1885 and its
targets under cold stress in \emph{B. napus}.}

\textbf{(A)} qRT-PCR analysis of the abundance of miR1885 in \emph{B.
napus} under 22 °C and 4 °C condition. \emph{U6} was used as endogenous
reference gene. \textbf{(B)} Alignment of miR1885 mature sequence
between \emph{B. rapa} and \emph{B. napus} . Short black line represents
the perfect matches. \textbf{(C)} Validation of miR1885 target genes and
cleavage sites using 5'-RACE PCR. Predicted TIR, NBS, and LRR domains
encoded by target genes are labeled (Left). Reverse complementary
matches of miR1885 and their target sites were showed and vertical
arrows indicated 5 termini of miRNA-guided cleavage products, as
identi\selectlanguage{english}fied by 5'-RACE, with frequency of clones shown (Right).
(\textbf{D} ) The transcript levels of \emph{Bn.TIR.A09}
and\emph{Bn.TNL.A03} in \emph{B. napus} under cold treatment. ``D''
represents days. \textbf{(E)} Relative transcript levels of cold-related
genes in wild-type \emph{B. napus} under 22 degC and 4 degC
condition.\emph{ACT7} was used as endogenous reference gene. Error bar
represents mean +- SE; \emph{n} = 3. * indicates \emph{P} \textless{}
0.05; ** indicates \emph{P} \textless{} 0.01; Student's t-test.

\textbf{Figure 5. Transgenic \emph{B. napus} plants overexpressing
miR1885.}

\textbf{(A)} Construct with \emph{bra-MiR1885} under control of AA6
promoter. pCaMV35S, CaMV35S promoter; Hyg R, hygromycin-resistance gene;
pAA, AA6 promoter; LB, left border of T-DNA; t35S, CaMV35S terminator;
tAA6, AA6 terminator. \emph{Eco} RI, \emph{Nco} I and \emph{Bst} EII are
restriction sites. \textbf{(B)} Plants of T\textsubscript{1} transgenic
lines and wild type in the field. Red scale bar = 10 cm. X2--X16,
transgenic plants. \textbf{(C)} Detection of transgene in rapeseed
transgenic lines by PCR using specific primer spanning AA6 promoter
and\emph{MiR1885} fragment. M, DNA marker. \textbf{(D)} Real-time PCR
showing expression of pri-miR1885 in transgenic lines. \emph{ACT7} was
used as endogenous reference gene. \textbf{(E)} Relative abundance of
mature miR1885 in transgenic lines. \emph{U6} was used as endogenous
reference gene\emph{.} \textbf{(F)} Northern blot analysis of
accumulation of mature miR1885 in \emph{MIR1885} -OE lines. \emph{U6}
was used as endogenous reference gene. \textbf{(G-H)} Relative
transcript levels of \emph{Bn.TIR.A09} \textbf{(G)} and
\emph{Bn.TNL.A03}\textbf{(H)} in miR1885-OE lines as compared to WT.
\emph{ACT7} was used as endogenous reference gene. Error bar represents
mean +- SE; \emph{n} = 3. * indicates \emph{P} \textless{} 0.05; **
indicates \emph{P}\textless{} 0.01; Student's t-test.

\textbf{Figure 6. Relative transcript levels of miR1885 and its targets
in transgenic and wild-type plants with or without cold treatment.}

\textbf{(A-C)} qRT-PCR analysis of RNA levels of miR1885 (\textbf{A}
),\emph{Bn.TNL.A03} \textbf{(B)} and \emph{Bn.TIR.A09} \textbf{(C)} in
transgenic plants after 22degC and 4 degC cold treatment for 30 days. Error
bar represents mean +- SD; \emph{n} = 3. * indicates \emph{P} \textless{}
0.05 and ** indicates \emph{P} \textless{} 0.01 under normal condition;
\selectlanguage{english}+ indicates \emph{P} \textless{} 0.05 and \selectlanguage{english}+\selectlanguage{english}+ indicates
\emph{P}\textless{} 0.01 under 4 degC. Student's t-test. \emph{U6} and
\emph{ACT7}were used as endogenous reference gene for checking miR1885
and targets expression levels, respectively.

\textbf{Figure 7. Effects of overexpressing and knocking down of~miR1885
in \emph{B. napus} after cold treatment.} (\textbf{A} ) Overexpression
of miR1885 increased the sensitivity of cold in \emph{B. napus.} The
cold treatment was 2 degC for 7 days. (\textbf{B} ) Expression of miR1885
and its targets in miR1885 knock down plants (STTM1885). Error bar
represents mean +- SD; \emph{n} = 3. * indicates \emph{P} \textless{}
0.05; ** indicates \emph{P} \textless{} 0.01; Student's t-test.\emph{U6}
and \emph{ACT7} were used as endogenous reference gene for measuring
miR1885 and targets expression levels, respectively. (\textbf{C} )
Knocking down of miR1885 increased the resistance to cold in \emph{B.
napus} as compared to its background. The cold treatment was 2 degC for 14
days. Red scale bar was 5 cm.

\textbf{Supplementary datasets}

\textbf{Table S1.} Distribution of small RNAs among different categories
of genomic loci. rRNA: ribosomal RNA; snRNA: small nuclear RNA; snoRNA:
small nucleolar RNA; tRNA: transfer RNA; unann: small RNAs without any
specific annotation.

\textbf{Table S2.} Differential expression analysis of identified known
miRNAs in \emph{Brassica rapa} by small RNA sequencing. RPM: reads per
million.

\textbf{Table S3.} Putative target genes of novel miRNAs in \emph{B.
rapa} . The miRNA target genes were predicted by psRNAtarget.

\textbf{Table S4.} Location information of pre-miRNAs in \emph{Brassica
napus} reference genome.

\textbf{Table S5.} List of primers used in this study.

\textbf{Supplementary Figures}

\textbf{Figure S1.} Proportion of small RNAs generated from different
categories of genomic loci in the two small RNA libraries
of\emph{Brassica rapa} .

\textbf{Figure S2.} Cis-element analysis of \emph{bra-MIR1885} promoter.
G-box, GA-motif, P-box, TC, TCT, AAAC, ATCT, Box 4, I-box and LTR
represent \emph{cis} -elements. A06: chrA06. 24219516 and 24221516:
Location information. Red or green bar represent position.

\textbf{Figure S3.} Analysis of the single nucleotide polymorphism (SNP)
density distribution, principal component analysis (PCA), and linkage
disequilibrium (LD) of \emph{B. napus} population. (\textbf{A} ) SNP
density distribution (number of SNPs in 0.1 Mb sliding windows across
each chromosome). (\textbf{B} ) PCA analysis. Blue, green, and red
points represent winter (W), semi-winter (SW) and spring (S) ecotypes,
respectively. (\textbf{C} ) Genome-wide average LD decay estimated
from\emph{B. napus} . (\textbf{D} ) LD heatmap of SNPs in
\emph{bna-MIR1885}promoter. Numbers indicate physical position where
\emph{bna-MIR1885}promoter locates. Deeper red color indicates stronger
linkage relationship. (\textbf{E} ) Comparison of the miR1885 relative
expression abundance between the two groups with extremely different
phenotypes. ES: Low-temperature sensitive accessions; ER:
Low-temperature resistant accessions. \emph{U6} was used as an internal
control. Data are presented as boxplot. Significant differences were
determined by Student's t-test: \emph{*} indicates \emph{p \textless{}
0.05} . (\textbf{F-I} ) The local LD heatmap of SNPs in
\emph{bna-MIR166C} , \emph{bna-MIR168A} and\emph{bna-MIR845A} gene loci.

\textbf{Figure S4. 5'-RACE assays of \emph{Bn. TIR.A09} (BnaA09g14980D)
and \emph{Bn. TNL.A03} (BnaA03g56180D).} PCR products
of\emph{Bn.A09.TIR} and \emph{Bn.A03.TNL} were separated in 2\% agarose
gel. DNA ladders are labeled on the side. PCR products were cloned into
T vectors for Sanger sequencing. Red line represents target fragment;
black line represents inner primer.

\textbf{Figure S5. Genomic synteny analysis between \emph{B. rapa}
and\emph{B. napus} species.}

\textbf{Figure S6. Phenotypes of transgenic plants overexpressing
miR1885 after cold treatment.} (\textbf{A} ) Phenotypes of transgenic
lines grown at 22 degC or 4 degC. Plants were grown 22 degC for 6 weeks (upper
panel), and then transferred to 4 degC for an additional 30 days (middle
panel), and then allowed to recover at 22 degC for 7 days. (\textbf{B} )
Post-germination phenotypes of miR1885-OE and WT K407 after cold
treatment. Seeds were grown on the germination bed at 22 degC for 2 days
and then transferred to 6 degC for 7 days. (\textbf{C} ) Seeds were grown
on the germination bed at 22 degC for 7 days. (\textbf{D} ) Quantitative
analysis of net hypocotyl lengths from WT K407 and \emph{MIR1885} -OE
lines. \emph{n} = 10-15 plants. Similar results were obtained in three
independent experiments.

\textbf{Figure S7. Freezing sensitivity of miR1885 over-expressing
transgenic plants.}

\textbf{(A)} Phenotypes of transgenic lines grown at 22 degC and \selectlanguage{english}-10 degC.
Plants were grown at 22 degC for 6 weeks (upper panel), and then
transferred to \selectlanguage{english}-10 degC for 2.5 h (middle panel), and then allowed to
recover at 22 degC for 7 days. \textbf{(B)} Real-time PCR showing the
expression pattern of miR1885 under freezing treatment. \emph{U6} was
used as endogenous reference gene. Error bar represents mean +-
SD;\emph{n} = 3. \textbf{(C)} Survival rates of transgenic lines and
wild-type plants after freezing treatment (\selectlanguage{english}-10 degC in a growth chamber
for 2.5 h, \emph{n} = 15--20 plants). * indicates \emph{P} \textless{}
0.05; ** indicates \emph{P \textless{}} 0.01 (Student's \emph{t}
-test).\textbf{(D)} Proline content in transgenic plants and wild type
after freezing treatment. Error bar represents mean +- SD. \emph{*}
indicates\emph{P} \textless{} 0.05, significant difference from plants
subjected to \selectlanguage{english}-10 degC before freezing treatment (0 h). \selectlanguage{english}+ indicates
\emph{P}\textless{} 0.05, significant difference from plants subjected
to \selectlanguage{english}-10 degC for 2.5 h (Student's \emph{t} -test).

\textbf{REFERENCE}

\textbf{Ahmed, W., Li, R., Xia, Y., Bai, G., Siddique, K.H.M, Zhang, H.,
Zheng Y., Yang X. Guo P.} (2022). Comparative Analysis of miRNA
Expression Profiles Between Heat-Tolerant and Heat-Sensitive Genotypes
of Flowering Chinese Cabbage Under Heat Stress Using High-Throughput
Sequencing. Genes \textbf{11:} 264.

\textbf{Alves-Junior, L., Niemeier, S., Hauenschild, A., Rehmsmeier, M.,
and Merkle, T.} (2009). Comprehensive prediction of novel microRNA
targets in Arabidopsis thaliana. Nucleic Acids Research
\textbf{37:}4010-4021.

\textbf{Anjali, N.N., and Sabu, K.K.} (2020). Role of miRNAs in Abiotic
and Biotic Stress Management in Crop Plants. In: Sustainable Agriculture
in the Era of Climate Change--Roychowdhury, R., Choudhury, S.,
Hasanuzzaman, M., and Srivastava, S., eds. Cham: Springer International
Publishing. 513-532.

\textbf{Axtell, M., Meyers, B. (2018).} Revisiting criteria for plant
microRNA annotation in the era of big data. Plant Cell\textbf{30:}
272--284.

\textbf{Bao, F., Huang, X., Zhu, C., Zhang, X., Li, X., and Yang,
S.}(2014). Arabidopsis HSP90 protein modulates RPP4-mediated
temperature-dependent cell death and defense responses. New
Phytologist\textbf{202:} 1320-1334.

\textbf{Bhalla, P.L., and Singh, M.B.} (2008). Agrobacterium-mediated
transformation of Brassica napus and Brassica oleracea. Nature
Protocols\textbf{3:} 181-189.

\textbf{Biggar, K.K., and Storey, K.B.} (2015). Insight into
post-transcriptional gene regulation: stress-responsive microRNAs and
their role in the environmental stress survival of tolerant animals.
Journal of Experimental Biology \textbf{218:} 1281-1289.

\textbf{Bolt, S., Zuther, E., Zintl, S., Hincha, D.K., and Schmulling,
T.} (2017). ERF105 is a transcription factor gene of Arabidopsis
thaliana required for freezing tolerance and cold acclimation. Plant,
Cell \& Environment \textbf{40:} 108-120.

\textbf{Chalhoub, B., Denoeud, F., Liu, S., Parkin, I.A.P., Tang, H.,
Wang, X., Chiquet, J., Belcram, H., Tong, C., Samans, B., Correa, M., Da
Silva, C., Just, J., Falentin, C., Koh, C.S., Le Clainche, I., Bernard,
M., Bento, P., Noel, B., Labadie, K., Alberti, A., Charles, M., Arnaud,
D., Guo, H., Daviaud, C., Alamery, S., Jabbari, K., Zhao, M., Edger,
P.P., Chelaifa, H., Tack, D., Lassalle, G., Mestiri, I., Schnel, N., Le
Paslier, M.-C., Fan, G., Renault, V., Bayer, P.E., Golicz, A.A., Manoli,
S., Lee, T.-H., Thi, V.H.D., Chalabi, S., Hu, Q., Fan, C., Tollenaere,
R., Lu, Y., Battail, C., Shen, J., Sidebottom, C.H.D., Wang, X.,
Canaguier, A., Chauveau, A., Berard, A., Deniot, G., Guan, M., Liu, Z.,
Sun, F., Lim, Y.P., Lyons, E., Town, C.D., Bancroft, I., Wang, X., Meng,
J., Ma, J., Pires, J.C., King, G.J., Brunel, D., Delourme, R., Renard,
M., Aury, J.-M., Adams, K.L., Batley, J., Snowdon, R.J., Tost, J.,
Edwards, D., Zhou, Y., Hua, W., Sharpe, A.G., Paterson, A.H., Guan, C.,
and Wincker, P.} (2014). Early allopolyploid evolution in the
post-Neolithic Brassica napus oilseed genome. Science
\textbf{345:}950-953.\textbf{Chen, S., Huang, H.-A., Chen, J.-H., Fu,
C.-C., Zhan, P.-L., Ke, S.-W., Zhang, X.-Q., Zhong, T.-X., and Xie,
X.-M.} (2020). SgRVE6, a LHY-CCA1-Like Transcription Factor From
Fine-Stem Stylo, Upregulates NB-LRR Gene Expression and Enhances Cold
Tolerance in Tobacco. Front Plant Sci \textbf{11} . 1276

\textbf{Chen, X.} (2009). Small RNAs and Their Roles in Plant
Development. Annual Review of Cell and Developmental Biology\textbf{25:}
21-44.

\textbf{Chen, X., Jiang, L., Zheng, J., Chen, F., Wang, T., Wang, M.,
Tao, Y., Wang, H., Hong, Z., Huang., Y., and Huang, R.} (2019). A
missense mutation in Large Grain Size 1 increases grain size and
enhances cold tolerance in rice. Journal of Experimental
Botany\textbf{70} : 3851-3866.

\textbf{Cong, R., Zhang, Z., and Lu, J.} (2019). Climate impacts on
yield of winter oilseed rape in different growth regions of the Yangtze
River Basin. China Journal Oil Crop Science \textbf{41} :894-903.

\textbf{Cui, C., Wang, J., Zhao, J., Fang, Y., He, X., Guo, H., and
Duan, C.} (2020). A Brassica miRNA Regulates Plant Growth and Immunity
through Distinct Modes of Action. Molecular Plant \textbf{13:} 231-245.

\textbf{Cui, N., Sun, X., Sun, M., Jia B., Duanmu, H,. Lv, D., Duan, X.,
and Zhu Y.} (2015). Overexpression of OsmiR156k leads to reduced
tolerance to cold stress in rice (Oryza Sativa). Molecular
Breeding\textbf{35} : 214.

\textbf{Cuperus, J., Carbonell, A., Fahlgren, N., Garcia-Ruiz, H.,
Burke, R., Takeda, A., Sullivan, C., Gilbert. S., Montgomery. T., and
Carrington, J.} (2010). Unique functionality of 22-nt miRNAs in
triggering RDR6-dependent siRNA biogenesis from target transcripts in
Arabidopsis. Nature Structure Molecular Biology \textbf{17} :997-1003.

\textbf{Dangl, J.L., and Jones, J.D.G.} (2001). Plant pathogens and
integrated defence responses to infection. Nature \textbf{411:} 826-833.

\textbf{Ding, Y., Shi, Y., Yang, S.} (2019). Advances and challenges in
uncovering cold tolerance regulatory mechanisms in plants. New
Phytologist \textbf{222:} 1690-1704.

\textbf{Fei, Q., Xia, R., and Meyers, B.} (2013). Phased, secondary,
small interferingRNAs in posttranscriptional regulatory networks. Plant
Cell \textbf{25} :2400-241.

\textbf{Gomez-Campo, C., and Prakash, S.} (1999). Origin and
domestication. In: Gomez-Campo C (ed) Biology of BrassicaCoenospecies.
Developments in Plant Genetics and Breeding 4. Elsevier, Amsterdam.

\textbf{Gupta, O.P., Meena, N.L., Sharma, I., and Sharma, P.} (2014).
Differential regulation of microRNAs in response to osmotic, salt and
cold stresses in wheat. Molecular Biology Reports \textbf{41:}4623-4629.

\textbf{Hussain, M., Luo, D., Zeng, L., Ding, X., Cheng, Y., Zou, X.,
Lv, Y., and Lu G.} (2022). Genome-wide transcriptome profiling revealed
biological macromolecules respond to low temperature stress in Brassica
napus L. Fronters in Plant Science \textbf{13} :1050995.

\textbf{He, J., Xu, M., Willmann, M.R., McCormick, K., Hu, T., Yang, L.,
Starker, C.G., Voytas, D.F., Meyers, B.C., and Poethig, R.S.} (2018).
Threshold-dependent repression of SPL gene expression by miR156/miR157
controls vegetative phase change in Arabidopsis thaliana. Plos
Genetics\textbf{14:} e1007337-e1007337.

\textbf{He, X., Fang, Y., Feng, L., and Guo, H.} (2008).
Characterization of conserved and novel microRNAs and their targets,
including a TuMV-induced TIR--NBS--LRR class R gene-derived novel miRNA
in Brassica. FEBS Letters \textbf{582:} 2445-2452.

\textbf{Huang, X., Li, J., Bao, F., Zhang, X., and Yang, S.} (2010). A
Gain-of-Function Mutation in the Arabidopsis Disease Resistance Gene
RPP4 Confers Sensitivity to Low Temperature. Plant
Physiology\textbf{154:} 796-809.

\textbf{Jiang, L., Fu, Y., Sun, P., Tian, X., and Wang, G.} (2022).
Identification of microRNA158 from Anthurium andraeanum and Its Function
in Cold Stress Tolerance. Plants \textbf{11} :3371.

\textbf{Jones-Rhoades, M.W., and Bartel, D.P.} (2004). Computational
Identification of Plant MicroRNAs and Their Targets, Including a
Stress-Induced miRNA. Molecular Cell \textbf{14:} 787-799.

\textbf{Jones, J.D.G., and Dangl, J.L.} (2006). The plant immune system.
Nature \textbf{444:} 323-329.

\textbf{Krishna, P., Sacco, M., Cherutti, J.F., and Hill, S.} (1995).
Cold-Induced Accumulation of hsp90 Transcripts in Brassica napus. Plant
Physiology \textbf{107:} 915-923.

\textbf{Li, F., Orban, R., Baker, B.} (2012). SoMART: a web server for
plant miRNA, tasiRNA and target gene analysis. Plant Journal\textbf{70}
:891-901.

\textbf{Li, S., Cheng, Z., and Peng, M.} (2020). Genome-wide
identification of miRNAs targets involved in cold response in cassava.
Plant Omics \textbf{13} : 57-64.

\textbf{Li, X., Lian, H., Zhao, Q., and He, Y.} (2019). MicroRNA166
Monitors SPOROCYTELESS/NOZZLE for Building of the Anther Internal
Boundary1. Plant Physiology \textbf{181:} 208-220.

\textbf{Liu, F., Xiong, X., Wu, L., Fu, D., Hayward, A., Zeng, X., Cao,
Y., Wu, Y., Li, Y., and Wu, G.} (2014). BraLTP1, a Lipid Transfer
Protein Gene Involved in Epicuticular Wax Deposition, Cell Proliferation
and Flower Development in Brassica napus. PLoS ONE \textbf{9} : e11027.

\textbf{Liu, J., Pang, X., Cheng, Y., Yin, Y., Zhang, Q., Su, W., Hu,
B., Guo, Q., Ha, S., Zhang, J., and Wan, H.} (2018). The Hsp70 Gene
Family in Solanum tuberosum: Genome-Wide Identification, Phylogeny, and
Expression Patterns. Scientific Reports \textbf{8:} 16628.

\textbf{Liu,Y., Teng, C., Xia, R., and Meyers, B.} (2020). PhasiRNAs in
Plants: Their Biogenesis, Genic Sources, and Roles in Stress Responses,
Development, and Reproduction. Plant Cell \textbf{32} :3059-3080.

\selectlanguage{polish}\textbf{Ljubej, V., Radojčić Redovniković, I., Salopek-Sondi, B.,
Smolko, A., Roje, S., and Šamec, D.}\selectlanguage{english} (2021). Chilling and Freezing
Temperature Stress Differently Influence Glucosinolates Content in
Brassica oleracea var. acephala. Plants \textbf{10:} 1305.

\textbf{Luo, T., Xian, M., Zhang, C., Zhang, C., Hu, L., and Xu,
Z.}(2019). Associating transcriptional regulation for rapid germination
of rapeseed (Brassica napus L.) under low temperature stress through
weighted gene co-expression network analysis. Scientific
Report.\textbf{9} :55-70.

\textbf{Ma, X., Zhao, F., Zhou, B.} (2022). The Characters of Non-Coding
RNAs and Their Biological Roles in Plant Development and Abiotic Stress
Response. International Journal of Molecular Science \textbf{23:} 4124.

\textbf{McClung, C.R., and Davis, S.J.} (2010). Ambient Thermometers in
Plants: From Physiological Outputs towards Mechanisms of Thermal
Sensing. Current Biology \textbf{20:} R1086-R1092.

\textbf{O'neill, C., Lu, X., Calderwood, A., Tudor, E., and Penfield,
S.} (2019). Vernalization and floral transition in autumn drive winter
annual life history in oilseed rape. Current Biology\textbf{29}
:4300-4306.

\textbf{Ozer, H} . (2003). Sowing date and nitrogen rate effects on
growth, yield and yield components of two summer rapeseed cultivars.
European Journal of Agronomy \textbf{19} :453--463.

\textbf{Paul, P., Chhapekar, S.S., Rameneni, J.J., Oh, S.H., Dhandapani,
V., Subburaj, S., Shin, S.Y., Ramchiary, N., Shin, C., Choi, S.R., and
Lim, Y.P.} (2021). MiR1885 Regulates Disease Tolerance Genes in Brassica
rapa during Early Infection with Plasmodiophora brassicae. International
Journal of Molecular Sciences \textbf{22:} 9433.

\textbf{Qian, W., Xiao, B., Wang, L., Hao, X., Yue, C., Cao, H., Wang,
Y., Li, N., Yu, Y., Zeng, J., Yang Y., Wang, X.} (2018). CsINV5, a tea
vacuolar invertase gene enhances cold tolerance in transgenic
Arabidopsis. BMC Plant Biology \textbf{18:} 228.

\textbf{Raza, A., Su, W.,Hussain, M., Mehmood, S., Zhang, X., Cheng, Y.,
Zhang, X., Cheng, Y., Zou, X. and Lv, Y} . (2021). Integrated analysis
of metabolome and transcriptome reveals insights for cold tolerance in
rapeseed (Brassica napus l.).~Fronters in Plant Science.~\textbf{12} :
721681.

\textbf{Sanghera, G.S., Wani, S.H., Hussain, W., and Singh, N.B.}(2011).
Engineering cold stress tolerance in crop plants. Curr
Genomics\textbf{12:} 30-43.

\textbf{Seo, P.J., Lee, A.K., Xiang, F., and Park, C.M.} (2008).
Molecular and functional profiling of Arabidopsis pathogenesis-related
genes: insights into their roles in salt response of seed germination.
Plant \& cell physiology \textbf{49:} 334-344.

\textbf{Shi, Y., Ding, Y., and Yang, S.} (2018). Molecular Regulation of
CBF Signaling in Cold Acclimation. Trends in Plant Science
\textbf{23:}623-6.

\textbf{Snider, C.S., Hsiang, T., Zhao, G., and Griffith, M.} (2000).
Role of Ice Nucleation and Antifreeze Activities in Pathogenesis and
Growth of Snow Molds. Phytopathology \textbf{90:} 354-361.

\textbf{Song, X., Yan, L., Cao, X., Qi,Y.} (2019). MicroRNAs and Their
Regulatory Roles in Plant-Environment Interactions. Annual Review of
Plant Biology \textbf{70:} 489-525.

\textbf{Sun, M., Shen, Y., Chen, Y., Wang, Y., Cai, X., Yang, J., Jia,
B., Dong, W., Chen, X., Sun, X.} (2022). Osa-miR1320 targets the ERF
transcription factor OsERF096 to regulate cold tolerance via JA-mediated
signaling. Plant physiology \textbf{189} :2500-2516.

\textbf{Tian, N., Liu, F., Wang, P., Yan, X., Gao, H., Zeng, X., and Wu,
G} . (2018). Overexpression of BraLTP2, a Lipid Transfer Protein of
Brassica napus, Results in Increased Trichome Density and Altered
Concentration of Secondary Metabolites. International Journal of
Molecular Sciences. \textbf{19} :1733.

\textbf{Thiebaut, F., Rojas c., A., Almeida, K., L., Grativol, C.,
Domiciano G., C., Lamb, C., R., Engler, J.,} \textbf{Hemerly A., S., and
Ferreira, P} . (2012). Regulation of miR319 during cold stress in
sugarcane. Plant cell and Environment \textbf{35} :502-512.

\textbf{Wang J., Czech B., and Weigel D} . (2009). miR156-regulated SPL
transcription factors define an endogenous flowering pathway in
Arabidopsis thaliana. Cell \textbf{138} :738-49.

\textbf{Wang, S., Sun, X., Hoshino, Y., Yu, Y., Jia, B., Sun, Z., Sun,
M., Duan, X., and Zhu, Y.} (2014). MicroRNA319 Positively Regulates Cold
Tolerance by Targeting OsPCF6 and OsTCP21 in Rice (Oryza sativa L.).
Plos One \textbf{9} . e91357

\textbf{Wang, X., and Wang, H., and Wang, J., and Sun, R., and Wu, J.,
and Liu, S., and Bai, Y., and Mun, J.-H., and Bancroft, I., and Cheng,
F., and Huang, S., and Li, X., and Hua, W., and Wang, J., and Wang, X.,
and Freeling, M., and Pires, J.C., and Paterson, A.H., and Chalhoub, B.,
and Wang, B., and Hayward, A., and Sharpe, A.G., and Park, B.S., and
Weisshaar, B., and Liu, B., and Li, B., and Liu, B., and Tong, C., and
Song, C., and Duran, C., and Peng, C., and Geng, C., and Koh, C., and
Lin, C., and Edwards, D., and Mu, D., and Shen, D., and Soumpourou, E.,
and Li, F., and Fraser, F., and Conant, G., and Lassalle, G., and King,
G.J., and Bonnema, G., and Tang, H., and Wang, H., and Belcram, H., and
Zhou, H., and Hirakawa, H., and Abe, H., and Guo, H., and Wang, H., and
Jin, H., and Parkin, I.A.P., and Batley, J., and Kim, J.-S., and Just,
J., and Li, J., and Xu, J., and Deng, J., and Kim, J.A., and Li, J., and
Yu, J., and Meng, J., and Wang, J., and Min, J., and Poulain, J., and
Wang, J., and Hatakeyama, K., and Wu, K., and Wang, L., and Fang, L.,
and Trick, M., and Links, M.G., and Zhao, M., and Jin, M., and
Ramchiary, N., and Drou, N., and Berkman, P.J., and Cai, Q., and Huang,
Q., and Li, R., and Tabata, S., and Cheng, S., and Zhang, S., and Zhang,
S., and Huang, S., and Sato, S., and Sun, S., and Kwon, S.-J., and Choi,
S.-R., and Lee, T.-H., and Fan, W., and Zhao, X., and Tan, X., and Xu,
X., and Wang, Y., and Qiu, Y., and Yin, Y., and Li, Y., and Du, Y., and
Liao, Y., and Lim, Y., and Narusaka, Y., and Wang, Y., and Wang, Z., and
Li, Z., and Wang, Z., and Xiong, Z., and Zhang, Z., and The Brassica
rapa Genome Sequencing Project, C.} (2011). The genome of the
mesopolyploid crop species Brassica rapa. Nature Genetics
\textbf{43:}1035-1039.

\textbf{Wang, Y., Zhang, Y., Wang, Z., Zhang, X., and Yang, S.} (2013).
A missense mutation in CHS1, a TIR-NB protein, induces chilling
sensitivity in Arabidopsis. The Plant Journal \textbf{75:} 553-565.

\textbf{Xin, Z., and Browse, J.} (1998). eskimo1 mutants of Arabidopsis
are constitutively freezing-tolerant. Proceedings of the National
Academy of Sciences \textbf{95:} 7799.

\textbf{Xu, P., and Cai, W.} (2019). Function of Brassica napus BnABI3
in Arabidopsis gs1, an Allele of AtABI3, in Seed Development and Stress
Response. Fronters in Plant Science \textbf{10} . 00067

\textbf{Xu, P., Zhu, Y., Zhang, Y., Jiang, J., Yang, L., Mu, J., Yu, X.,
and He, Y.} (2021). Global Analysis of the Genetic Variations in
miRNA-Targeted Sites and Their Correlations with Agronomic Traits in
Rapeseed. Fronters in Genetics \textbf{12} . 741858.

\textbf{Yan, J., Zhao, C., Zhou, J., Yang, Y., Wang, P., Zhu, X., Tang,
G., Bressan, R.A., and Zhu, J.-K.} (2016). The miR165/166 Mediated
Regulatory Module Plays Critical Roles in ABA Homeostasis and Response
in Arabidopsis thaliana. Plos Genetics \textbf{12:} e1006416.

\textbf{Yang, H., Shi, Y., Liu, J., Guo, L., Zhang, X., and Yang,
S.}(2010). A mutant CHS3 protein with TIR-NB-LRR-LIM domains modulates
growth, cell death and freezing tolerance in a temperature-dependent
manner in Arabidopsis. The Plant Journal \textbf{63:} 283-296.

\textbf{Yu, X., Wang, H., Lu, Y., de Ruiter, M., Cariaso, M., Prins, M.,
van Tunen, A., and He, Y.} (2011). Identification of conserved and novel
microRNAs that are responsive to heat stress in Brassica rapa. Journal
of Experimental Botany \textbf{63:} 1025-1038.

\selectlanguage{ngerman}\textbf{Zbierzak, A.M., Porfirova, S., Griebel, T., Melzer, M., Parker,
J.E., and Dörmann, P.} (2013). A TIR--NBS protein encoded by Arabidopsis
Chilling Sensitive 1 (CHS1) limits chloroplast damage and cell death at
low temperature. The Plant Journal \textbf{75:} 539-552.

\textbf{Zeng, X., Xu, Y., Jiang, J., Zhang, F., Ma, L., Wu, D., Wang,
Y., and Sun, W.} (2018). Identification of cold stress responsive
microRNAs in two winter turnip rape (Brassica rapa L.) by high
throughput sequencing. BMC Plant Biology \textbf{18:} 52.

\textbf{Zhang, X., Li, X., and Liu, J.} (2014). Identification of
Conserved and Novel Cold-Responsive MicroRNAs in Trifoliate Orange
(Poncirus trifoliata (L.) Raf.) Using High-Throughput Sequencing. Plant
Molecular Biology Reporter \textbf{32} : 328-341.

\textbf{Zhang, X., Sun, J., You, Y., Song, J., and Yang, Z.} (2018).
Identification of Cd-responsive RNA helicase genes and expression of a
putative BnRH 24 mediated by miR158 in canola (Brassica napus).
Ecotoxicology and Environmental Safety \textbf{157:} 159-168.

\textbf{Zhou, M. and Tang, W.} (2019). MicroRNA156 amplifies
transcription factor-associated cold stress tolerance in plant cells.
Molecular Genetics and Genomics \textbf{294:} 379-393.

\textbf{Zhu, M., Wang, X., Zhou, Y., Tan, J., Zhou, Y., Gao, F.} (2022).
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.
Fronters in Genetics 13:870446.

\textbf{Zhang, S., Li, L., and Zhang, C.} (2012). Effects of sowing date
and planting density on the seed yield and oil content of winter oilseed
rape. Ying Yong Sheng Tai Xue Bao \textbf{23} :1326-1332.

\textbf{Table 1. The miRNAs that were down- or up-regulated by more than
1.5-fold change under cold treatment.}\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
RPM: Reads per million. & RPM: Reads per million. & RPM: Reads per
million. & RPM: Reads per million. & RPM: Reads per million. & RPM:
Reads per million.\tabularnewline
\textbf{miRNA} & \textbf{RPM\_22} & \textbf{RPM\_4 } & \textbf{Ratio}
& \textbf{Up/Down} & \textbf{Evidence}\tabularnewline
bra-miR157a\textasciitilde{}c & 2500.1 & 765.55 & 0.31 & Down &
qPCR\tabularnewline
bra-miR161 & 251.12 & 87.89 & 0.35 & Down & qPCR\tabularnewline
bra-miR393f & 186.74 & 68.64 & 0.37 & Down & qPCR\tabularnewline
bra-miR396b & 1431.99 & 726.84 & 0.51 & Down & Chen et al.,
2019\tabularnewline
bra-miR319c & 289.71 & 149.59 & 0.52 & Down & Thiebaut et al.,
2012\tabularnewline
bra-miR156k & 9852.92 & 5746.05 & 0.58 & Down & Cui et al.,
2015\tabularnewline
bra-miR1140 & 7314.42 & 4436.69 & 0.61 & Down & Northern
blot\tabularnewline
bra-miR158a & 168.42 & 103.83 & 0.62 & Down & Jiang et al.,
2022\tabularnewline
bra-miR156a\textasciitilde{}f & 2023.64 & 1322.56 & 0.65 & Down & Wang
et al., 2009\tabularnewline
bra-miR159a & 18327.95 & 12894.73 & 0.7 & Down & Li et al.,
2020\tabularnewline
bra-miR168d & 2498.35 & 3808.66 & 1.52 & Up & qPCR; Northern
blot\tabularnewline
bra-miR403 & 671.93 & 1004.9 & 1.5 & Up & Zhang et al.,
2004\tabularnewline
bra-miR390a\textasciitilde{}b & 251.32 & 417.09 & 1.66 & Up &
/\tabularnewline
bra-miR398b\textasciitilde{}c & 1000.04 & 2381.13 & 2.38 & Up & qPCR;
Northern blot\tabularnewline
bra-miR408 & 167.55 & 444.37 & 2.65 & Up & qPCR\tabularnewline
bra-miR1885 & 79.41 & 328.22 & 4.13 & Up & qPCR\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 1.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

Figure 2.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

Figure 3.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Figure 4.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

Figure 5.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

Figure 6.tiff\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

Figure 7.tiff

\selectlanguage{english}
\FloatBarrier
\end{document}

