Identification of regulatory pathways of miRna and phasiRna related to flax lignan synthesis

ABSTRACT Lignan is an important functional health care substance. At present, research on the molecular mechanism of flax lignan biosynthesis is primarily focused at the transcriptional level, with no research on the regulatory mechanism of small RNAs, particularly phased small interfering RNAs (phasiRNAs) related to lignan content. In this study, we used flax seeds at different developmental stages as materials to construct small RNA libraries and degradation libraries. Eighteen known and 4 new miRNA families containing 40 and 4 miRNA members, respectively, were identified. A total of 104 target genes were identified. Using small RNA data for a genome-wide search, a total of twenty-two 21-nt PHAS loci and twenty-eight 24-nt PHAS loci were obtained. Five miRNAs capable of triggering target genes to produce phasiRNAs were identified. These 5 miRNAs triggered their target genes to produce 39 phasiRNAs. Three complete miRNA-PHAS loci-phasiRNA-target gene functional regulatory pathways were identified. Furthermore, qRT‒PCR confirmed an obvious regulatory relationship among the regulatory pathways. In this study, the molecular regulatory network among flax miRNAs, phasiRNAs and target genes was described from the perspective of lignan synthesis, providing an important scientific basis for the study of flax molecular breeding for high lignan content.


Introduction
Flax (Linum usitatissimum L.) is an important cash crop in China. Research shows that linseed is rich in a variety of functional nutrients such as α-linolenic acid, gum, protein, and lignan (Ezzat et al. 2018). Among these nutrients, the content of lignan in linseed is very high, dozens or even hundreds of times higher than that of other plants (Blitz, Murphy, and Au 2007;Kiely et al. 2003;Meagher and Beecher 2000). Therefore, flax is also known as the "king of lignan." Lignan is a plant polyphenol that is also called phytoestrogen because its structure is similar to human estrogen. The development and utilization of lignan as a natural plant extract has been highly valued. Studies have shown that lignan has a significant anticancer effect and slows a series of chronic diseases, making it the focus of an active research area in natural nutritional and health care substances (Dobrowolska and Regulska-Ilow 2021). Therefore, the development and comprehensive utilization of linseed phenolics is indispensable to the sustainable development of linseed and can greatly improve the added value of linseed products.
Noncoding RNA refers to RNA that is not translated into proteins in the transcriptome and includes housekeeping noncoding RNA and regulatory noncoding RNA (Vaucheret et al. 2004). The regulation of short noncoding RNA (small noncoding RNA) is a recent research hotspot that includes microRNA (miRNA) and phased small interfering RNAs (phasiRNAs). In addition to directly silencing target genes and regulating plant growth and development, there is a class of 22-nt miRNAs that can target genes and trigger cutting under the joint action of a variety of enzymes and producing phase-aligned secondary phasiRNAs. These 21-nt phasiRNAs have regulatory functions similar to those of miRNAs. The phasiRNA originally reported in Arabidopsis was named trans-acting siRNA (tasiRNA) due to its trans actions on downstream target genes. There are four types of tasiRNAs triggered by miR173, miR390 and miR828 in Arabidopsis, each of which have their own functions (Fei, Xia, and Meyers 2013). tasiRNA produced by TAS1 targets and negatively regulates two genes, HTT1 and HTT2, which improve heat resistance and participate in the response to heat stress in Arabidopsis (Li et al. 2014). tasiRNA produced by TAS3 and targeting the auxin response factor is involved in regulating the development time and morphology of Arabidopsis, particularly leaf morphogenesis (Adenot et al. 2006;Fahlgren et al. 2006). phasiRNA, miRNA, and endogenous target mic (ETM) coregulate lipid synthesis in soybean (Ye et al. 2014). tasiRNA produced by TAS4 targets the MYB transcription factor to participate in the biosynthesis of grape nutrients (Rock 2013). Thus, these phasiRNAs play important regulatory roles in plant stress resistance, growth and development, and nutrient synthesis (lipids, anthocyanins, etc.).
Although some progress has been made in research on the molecular mechanism of flax lignan synthesis, it has primarily been focused at the mRNA level with little research on the molecular interactions with small molecule RNAs (miRNAs and phasiRNAs). Therefore, this study combined high-throughput sequencing technology and bioinformatics analysis methods to analyze the regulatory mechanism and network of miRNAs and phasiRNAs related to flax lignan synthesis at the genome-wide level. The transcriptome and sRNA sequencing data of the flax varieties "SY 4" and "NEW" were analyzed at 10, 20 and 30 d after anthesis. The phasiRNAs mediated by flax seed lignan synthesis-related miRNAs were screened. Based on the differential expression of the screened miRNAs, the triggered phasiRNAs and target genes in the two flax varieties exhibited significant differences in lignan content. This study provides an important theoretical basis for the molecular breeding of high-lignan flax varieties by preliminarily clarifying the regulatory role of sRNAs and the functional regulatory pathway of the synthesis of lignan in flax seed.

Plant materials
In this study, the flax varieties "SY 4" and "NEW" were sown in the Minzhu experimental field of Heilongjiang Academy of Agricultural Sciences (N 45° 51", E 126° 48"). Plants with good, uniform growth trends were selected for subsequent experiments. The seeds of flax were collected at three developmental stages, i.e., 10, 20 and 30 d postanthesis (DPA). Each stage included three biological replicates, each of which was represented by three individual flax plants, for a total of 18 samples. To ensure consistency, each sample was composed of 10-15 uniform flax seeds. Two subsamples for each developmental stage (10, 20 and 30 DPA) were collected, one for library construction and sequencing and one for lignan content determination. After a quick incubation in liquid nitrogen, they were stored at −80°C until use.

RNA-seq and small RNA library construction
Each sample was ground to fine powder in liquid nitrogen, and the total RNA from flax seed at each developmental stage was extracted using an Axyprep total RNA small preparation kit (Axygen Company, USA) following the instructions of the kit. The total RNA was enriched for small RNA and extracted for small RNA library construction. The small RNA fraction (18 to 30 nt in length) was gel-extracted, ligated to 5' and 3' adaptors and purified. These small RNAs were reverse transcribed by RT-PCR and finally amplified via PCR. For RNA-seq libraries, after total RNA extraction and DNase I treatment, mRNA was isolated using magnetic beads, and oligo (dT) was used for isolation. After mixing with fragment buffer, the mRNA was divided into short fragments. Complementary chains were synthesized using mRNA fragments as templates. Agarose gel electrophoresis was used to select suitable fragments as templates for PCR amplification. All microRNAs and RNA-seq libraries were sequenced on the Illumina HiSeq 2000 platform (Beijing, China).

Small RNA data analysis
Small RNA sequencing data were preprocessed by removing adapters and then mapped to the flax reference genome to obtain the distribution of tags within the genome (https://phytozome.jgi.doe.gov/pz/ portal.html) (Wang et al. 2012). In brief, small RNA sequencing data from different libraries were combined to increase the sequencing depth for PHAS locus identification. A phasing score of 25 was used as a stringent cutoff, followed by a manual check to remove loci producing highly abundant small RNAs of other sizes, which are most likely degradation products from t/rRNAs. The overall phasiRNA abundance for each PHAS locus was calculated by summing the normalized abundance of 21-nt or 24-nt small RNAs generated from each corresponding 21-PHAS and 24-PHAS locus. PHAS locus identification was performed using the same method previously described by Zhai et al. (2011). The naming convention for phasiRNA includes the reference sequence ID, identification added by analysis software and positive and negative chain information; for example, scaffold 200_4_756 (-) is the phasiRNA generated from 3' to 5' on the negative chain of the scaffold 200 reference sequence. Raw miRNA sequences were deposited in the National Center for Biotechnology Information (NCBI) database and can be accessed in the database (https://www.ncbi.nlm.nih.gov/) under accession number PRJNA478806.

Expression analyses of small RNAs and their target genes
To analyze the expression levels of miRNAs, phasiRNAs and target genes, stem-loop RT-PCR and qRT-PCR were both used in this study. Stem-loop RT-PCR was performed according to a previously published method (Xie et al. 2021). U6 (Phytozome gene ID: Lus10035151) was used as an internal control for stemloop RT-PCR (Kou et al. 2012). 18S (NCBI Accession: EU307117) was the internal control for qRT-PCR. Stem-loop RT primers and qPCR primers are listed in Table S1. The flax seed small RNA templates from 10 DPA, and the relative expression levels were calculated using 2 −∆∆CT , where ∆∆C T = (C T, target - (Livak and Schmittgen 2001). All experiments were performed with three biological replicates, each with three technical replicates under the same conditions.

Determination of lignan content
The seeds of the "SY 4" and "NEW" flax cultivars were obtained at 10, 20 and 30 DPA with three biological replicates prepared for each sample type. Fixation was performed at 105°C for 30 min, and the samples were then incubated at 80°C until they reached a constant weight. The dried fixed samples were then ground into a powder and passed through a #100 mesh. The resulting powder was stored in a desiccator until further use.
In this study, the high-performance liquid chromatography (HPLC) procedure of Charlet et al. (2002) (with minor modifications) was used to quantify the lignan content of the flax seeds. HPLC was performed using an Agilent 1260 HPLC and a diode array detector (DAD). The mobile phases were methanol (A) and ultrapure water (B). The chromatography column was an Agilent C18 (4.6 mm × 250 mm × 5 μm), and the temperature of the column was 35°C. The detection wavelength was 290 nm, and 20 μL of sample was injected into the column.

Accumulation of lignans during flax seed development
The accumulation of lignan at three stages of seed development (10, 20 and 30 DPA) in the flax varieties "SY 4" and "NEW" was studied. The lignan content in flaxseed was relatively low before 10 DPA (Ghose et al. 2014). The lignan content increased steadily from 10 DPA to 30 DPA in the two varieties. However, the lignan content in "SY 4" seeds was significantly higher than that in "NEW" seeds at all three stages. (Figure 1).

Mining of phasiRNA related to flax seed lignan synthesis
From the data on the sRNA and degradation groups, we obtained the degradation site information for 22-nt miRNAs and their corresponding target genes. Starting from the degradation site of the target gene, phasiRNA was formed by phase cutting units of 21 nt and 24 nt at the 3' end (the negative chain has two base extensions in the 3' direction).
Using the data from the small RNA library and degradation group library, flax miRNAs and target genes, a total of eighteen known (miR1511, miR156, miR157, miR159, miR162, miR164, miR165, miR166, miR167, miR168, miR171, miR172, miR319, miR390, miR396, miR397, miR399, miR408) and 4 new miRNA families (unconservative_scaffold43_4703, unconservative_scaffold373_19_1213(-), unconservative_scaf-fold1917_14_1162(+), unconservative_scaffold1199_40705) were obtained, which included 40 and 4 miRNA members, respectively. A total of 104 target genes were identified (Table 1). Based on the annotation of the nr database, the target genes were primarily cytochrome P450/ABC transporter, multicopper oxidase, indigoidine synthase a-like protein, MBD9, laccase and genes related to various metabolic pathways. In addition to the conserved target genes, the known miRNAs identified new target genes in flax, indicating that the known miRNAs not only maintained their basic function during evolution but also obtained new target genes to acquire new functions (Table 1). Among these miRNAs, there were five miRNAs capable of triggering one or more target genes to produce phasiRNA, namely aly miR390a-5p, gma-miR390e, cme- miR172a, scaffold373_19_1213(-) and scaffold 917_14_1162(+) ( Table S2). These five miRNAs jointly triggered their target genes to produce a total of 39 phasiRNAs, and the expression levels of these phasiRNAs at each stage of seed development was distinct for the two flax varieties (Table S3). Among them, the miR390 family exhibited the most prominent ability to trigger its target gene to produce phasiRNA, resulting in a total of 29 phasiRNAs. Among the phasiRNAs triggered by aly miR390a-5p and gma-miR390e, scaffold67_4_121(-), scaffold931_5_440(+), scaffold931_5_417(-), scaffold200_4_504(-), scaffold1686_4_315 (-), scaffold701_31_102 (+) and scaffold 406_1_671(+) were expressed at more than 10,000 RPM at different stages of seed development in the two different flax varieties. Scaffold200_4_800(+), scaffold701_31_102(+), scaffold416_16_270(+), scaffold416_16_226(-), triggered by aly miR390a-5p and gma-miR390e, and scaf-fold373_19_1213(-), triggered by scaffold38_1_294(-), had significantly higher expression in the seeds of "SY 4" at different developmental stages than those in the seeds of "NEW". These phasiRNAs and their target genes may be involved in the accumulation of lignan in flax seeds.

Identification of PHAS loci related to flax seed lignan synthesis
phasiRNA is a type of plant-specific secondary sRNA that depends on the cleavage of its transcript precursors by miRNA. To explore the PHAS locus involved in lignan synthesis in flax, 18 sRNA sequencing samples of flax were compared with the flax genome to predict the PHAS locus, and then the PHAS locus was screened in the phenylalanine biosynthesis and metabolism pathway, which is the pathway that begins lignan derivation. Using screening criteria of a phase score greater than 5 and a multiple test p value less than 0.05, twenty-two 21-nt PHAS loci (Table S4) were predicted of which 7 PHAS loci were determined to trigger phasiRNA/miRNA. Scaffold917_14 and scaffold373_19 were triggered by the corresponding phasiRNA. Scaffold701_31, scaffold1099_2, scaffold701_34 and scaffold 863_21 were triggered by multiple new miRNAs. Only scaffold416_16 was triggered by the known miRNAs aly miR390a-5p and gma-miR390e (Table 2). In this study, twenty-seven 21-nt phasiRNAs were predicted to target 47 genes. Among them, 11 phasiRNAs targeted only one gene, whereas 5 phasiRNAs targeted up to 4 genes (Table S5). Twenty-eight 24-nt PHAS loci were predicted (Table S6) of which only 2 PHAS loci were determined to trigger phasiRNAs/miRNAs. Scaffold695_1 was triggered by the corresponding phasiRNA. Scaffold695_2 was triggered simultaneously by the corresponding phasiRNA and a new miRNA (Table 2), but no PHAS locus was triggered by known or new miRNAs. In this study, three 24-nt phasiRNAs were predicted to target six genes. One gene was targeted by scaffold563_2_241, three genes were targeted by scaffold 947_11_324, and two genes were targeted by scaffold 947_ 11_ 516 (Table S7).

Functional regulatory pathways related to the synthesis of lignan by phasiRNA
In the process of phasiRNA mining of two varieties with significant differences in lignan content and samples at different stages of flax seed development, 21-nt reproducible phasiRNAs showed that there were 22 functional regulatory pathways of miRNA-PHAS loci-phasiRNA-target genes (Table S8). After screening, three complete regulatory pathways were obtained. The first regulatory pathway was the PHAS locus scaffold917_14 and could be targeted by scaffold 373_19_1213(-) and then generated scaffold 917_14_1162 (+) of phasiRNA; the target gene of scaffold 917_14_1162(+) was Lus10030865, which is annotated as methyl CpG binding domain 9 (MBD9). The second regulatory pathway was the PHAS locus scaffold 373_19, which could be targeted by scaffold 917_14_1162(+) and then generate scaffold 373_19_1213(-) of phasiRNA; the target gene of scaffold 373_19_1213(-) was Lus10030625, which is annotated as methyl CpG binding domain 9 (MBD9). The third regulatory pathway was the PHAS locus scaffold 416-16, which can be targeted by gma-miR390e and then generate scaffold 416_16_226(-) of phasiRNA; the target gene of scaffold 416_16_226(-) was Lus10026400, which is annotated as laccase (Table 3). A total of 28 functional regulatory pathways were found for 24-nt reproducible phasiRNAs (Table S9), but the complete miRNA-PHAS locus-phasiRNA-target gene regulatory pathway was not screened.

Expression analysis of small RNAs and target genes in flax varieties with different lignan contents
The expression patterns of miRNAs, phasiRNAs and target genes in the three complete regulatory pathways were analyzed. Using the seeds of two flax varieties, "SY 4" and "NEW," as experimental materials, we analyzed the miRNAs, phasiRNAs and target genes in the three complete regulatory pathways by qRT-PCR (Figure 2). Scaffold 917_ 14_1162(+) and scaffold 373_19_1213(-) targeted each other in the complete regulatory pathway related to the accumulation of lignan. Using qRT-PCR analysis, it was found that their expression levels were reduced in the seeds of the high lignan variety "SY 4" compared with those of the low lignan variety "NEW," whereas the expression of their corresponding target genes Lus1030865 and Lus1030625 were increased trend in the seeds of the high lignan variety "SY 4" compared with that of the low lignan variety "NEW" (Figure 2a, b). The scaffold 917_14_1162(+) for Lus10030865 and scaffold373_19_1213(-) for Lus10030625 showed significant negative regulation. Both target genes, Lus10030865 and Lus10030625, were annotated functionally as MBD9. Gma-miR390e targeted scaffold 416_16_226(-), and their expression levels were upregulated in the seeds of the high lignan variety "SY 4" compared with those of the low lignan variety "NEW," and expression of the corresponding target gene Lus10026400 was downregulated in the seeds of the high lignan variety "SY 4" compared with the low lignan variety "NEW" (Figure 2c). This indicates that scaffold 416_16_226(-) had an obvious negative regulatory effect on Lus10026400, and the function of the target gene Lus10026400 was annotated as laccase. Based on the above results, it is speculated that the three miRNAs in this study regulate the accumulation of flax lignan; scaffold 373_ 19_1213 (-) and scaffold 17_14_1162 (+) inhibited the accumulation of lignan, whereas gam-miR390e promoted the accumulation of lignan.

Discussion and conclusions
Using high-throughput sequencing data to screen miRNAs related to flax lignan synthesis, a total of 22 miRNA families were identified, including 18 known miRNA families and 4 new miRNA families. A total of 44 miRNAs in these miRNA families were significantly different between the high-lignan flax variety "SY 4" and the low-lignan flax variety "NEW" at different stages of flax seed development (Table 1). After studying their regulatory networks, it is suggested that these miRNAs may primarily be involved in the synthesis of flax lignan through negative regulation of cytochrome P450, multicopper oxidase, indigoidine synthase a like protein, MBD9, laccase and genes related to various metabolic pathways. The above target genes, such as cytochrome P450, laccase and multicopper oxidase, have been reported to be related to plant secondary metabolism (Kato and Esaka 2000;Mao et al. 2020;Tu et al. 2020;Zhou et al. 2021). These results demonstrate that the miRNAs identified here are likely to participate in the accumulation of flaxseed lignan by negative regulation of these target genes.
In addition to the screening of miRNAs related to flax lignan synthesis, we also showed that there is a complex phasiRNA regulatory network in flax, which involves a large number of coding and noncoding PHAS loci (Table S4, Table S6). A large amount of sRNA data was used for a detailed analysis of the PHAS loci related to flax lignan synthesis. We identified several phasiRNA regulatory pathways (Table S8, Table S9). However, compared with a large number of lignan content-related PHAS loci, the conserved phasiRNA pathway identified as miRNA-mediated represented only a small part of the phasiRNAs. Therefore, we inferred that the mechanism of phasiRNAs may be more complex than we currently understand. The tasiRNAs triggered by miR173, miR390 and miR828 in Arabidopsis have their own functions (Fei, Xia, and Meyers 2013). miR173 and miR828  may trigger the target gene to produce phasiRNA in Arabidopsis thaliana and were not detected in the sRNA sequencing of flax seeds in this study, whereas miR390 was significantly expressed in flax seed, triggering the production of phasiRNA during flax seed development. In other words, only miR390 in flax, which is homologous to Arabidopsis thaliana, can trigger its target gene to produce phasiRNA. It is currently unclear whether there are any other nonconserved miRNAs in flax that can trigger their target gene to produce phasiRNA. In this study, miRNAs in flax seeds during development were comprehensively identified using high-throughput sequencing technology, and phasiRNAs that negatively regulate target genes involved in flax lignan synthesis through translational inhibition and posttranscriptional cleavage were mined from the sequencing data. Three miRNA-PHAS loci-phasiRNA-target gene regulatory pathways were screened and identified, which may be involved in the synthesis of lignan in flax seeds (Table 3). Laccase has been reported to probably be involved in monolignol oxidation of lignin and lignan biosynthesis. However, Yonekura-Sakakibara et al. (2020) revealed that AtLAC5 was also essential for lignan biosynthesis. Laccase has been reported to be involved in the synthesis and metabolism of lignan in vitro (Mattinen et al. 2011;Tranchimand et al. 2006). The physiological and molecular genetic research data of Arabidopsis clearly showed that AtMBD9 positively regulates the lateral branch development and flowering time of Arabidopsis by participating in the double epigenetic pathways of DNA methylation and histone acetylation (Yaish, Peng, and Rothstein 2009). There are few studies on the function of the target gene MBD9, and there is no research on its participation in lignan synthesis at present. However, the target gene of the two regulatory chains related to lignan synthesis screened in this study was MBD9. It is speculated that this gene likely functions as a positive regulator of lignan synthesis. qRT-PCR was used to verify the expression in the seeds at different development stages of two flax materials with significant differences in lignan content. The expression patterns showed that miRNA positively regulated phasiRNA, whereas phasiRNA negatively regulated its target genes. Interestingly, scaffold 373_19_1213(-) and scaffold 917_14_1162(+) could trigger each other, simultaneously playing the roles of miRNA and phasiRNA. Flax lignan accumulation is an extremely complex physiological process. The phasiRNA mediated by miRNA related to lignan synthesis participates in many physiological and biochemical processes, including lignan synthesis. In addition to participating in the formation of secondary phasiRNAs, miRNA triggers are also directly involved in the accumulation of lignan via negative regulation. By screening and analyzing the miRNAs and PHAS loci related to flax lignan synthesis, this study preliminarily defined the specific functions of specific miRNA triggers and further explored the detailed regulatory process of phasiRNAs on flax lignan accumulation. These data could provide genetic resources for breeding high-lignan flax varieties through targeted molecular breeding.

Funding
This work was supported by the National Natural Science Foundation of China (31801409); Innovation and Entrepreneurship Training Program for College Students (202210304103Y); Social Livelihood Science and Technology Project of Nantong City (MS22021029).

Ethical approval
There are no ethical issues involved in this article.