Integration of miRNAs, Degradome, and Transcriptome Omics Uncovers a Complex Regulatory Network and Provides Insights Into Lipid and Fatty Acid Synthesis During Sesame Seed Development

Sesame (Sesamum indicum L.) has always been known as a health-promoting oilseed crop because of its nutrient-rich oil. In recent years, studies have focused on lipid and fatty acid (FA) biosynthesis in various plants by high-throughput sequencing. Here, we integrated transcriptomics, small RNAs, and the degradome to establish a comprehensive reserve intensive on key regulatory micro RNA (miRNA)-targeting circuits to better understand the transcriptional and translational regulation of the oil biosynthesis mechanism in sesame seed development. Deep sequencing was performed to differentially express 220 miRNAs, including 65 novel miRNAs, in different developmental periods of seeds. GO and integrated KEGG analysis revealed 32 pairs of miRNA targets with negatively correlated expression profiles, of which 12 miRNA-target pairs were further confirmed by RT-PCR. In addition, a regulatory co-expression network was constructed based on the differentially expressed gene (DEG) profiles. The FAD2, LOC10515945, LOC105161564, and LOC105162196 genes were clustered into groups that regulate the accumulation of unsaturated fatty acid (UFA) biosynthesis. The results provide a unique advanced molecular platform for the study of lipid and FA biosynthesis, and this study may serve as a new theoretical reference to obtain increased levels of UFA from higher-quality sesame seed cultivars and other plants.


INTRODUCTION
Sesame (Sesamum indicum L.), as one of the oldest herbaceous diploid oilseed plants (2n = 26) in the Pedaliaceae family with a genome size of approximately 369 Mb, is mainly cultivated in tropical and subtropical regions of Asia, Africa, and South America Majdalawieh et al., 2020). Among the major cultivated oil crops worldwide, sesame as a commercial material contains the highest number of fatty acids (FAs) (Wacal et al., 2019;Zhang et al., 2020), and among the FA types, linoleic, oleic, palmitic, and stearic acids are the main components that make up the total FA configuration of sesame oil (Elleuch et al., 2007;Adebowale et al., 2011;Bhunia et al., 2015). In previous studies, unsaturated fatty acids (UFAs) are considered as a quality indicator of edible oils in the core area of research for the development of functional and industrial products Ahmed et al., 2020), which can prevent cardiovascular sclerosis and reduce blood pressure and blood lipids (Li et al., 2015). Therefore, increasing the UFAs content is a conventional way to improve the quality of sesame varieties and increase the seed oil content, which are also identified as the most important strategies to meet the current needs of consumers (Wang D.D. et al., 2018).
With the advent of high-throughput sequencing, genome breeding has developed rapidly, such as quantitative genetics, genetic breeding, genetic relatedness, and diversity Yan et al., 2016;Kumar et al., 2017). In recent years, although the whole transcriptome sequencing of sesame has been described, information on the role of micro RNAs (miRNAs) and target genes, especially miRNA expression profiles involved in lipid and FA biosynthesis of sesame oil, is still poorly understood. miRNAs as endogenous non-coding RNAs (18-26 nt) regulate gene expression through post-transcriptional translational inhibition or direct cleavage of transcripts based on complementarity between miRNA and its target gene (La Sala et al., 2018;Sevgi, 2018). Recent studies have demonstrated the role of miRNAs in almost all life activities, such as biosynthesis, cell reproduction, plant resistance to biotic and abiotic stresses, and metabolic pathways (Han et al., 2016;Zhao et al., 2018Zhao et al., , 2020Zheng et al., 2019). In banana fruits, 128 known miRNAs belonging to 42 miRNA families were determined, and 12 new miRNAs were identified; in olive, 135 conserved miRNAs were identified and 38 putative new miRNAs were discovered; in strawberry, a total of 88 known miRNAs and 103 targets cleaved by 19 known miRNAs families were discovered (Xu et al., 2013;Yanik et al., 2013;Dan et al., 2018). Notably, from 18 sesame libraries, 351 known and 91 novel miRNAs, 116 miRNAs were reported to regulate salt stress response (Zhang et al., 2020). Similarly, major miRNA-target pairs and the complex regulatory network responsible for terpenoid biosynthesis in tea leaves (Zhao et al., 2018) and cadmium phytoremediation in hyperaccumulator Sedum alfredii were revealed by integrating small RNAs, degradome, and transcriptomics (Han et al., 2016). Subsequently, the lipid biosynthetic pathway and key periods of FA biogenesis and associated proteins were deciphered in the developing endosperm of Cocos nucifera (Reynolds et al., 2019) and tree peony seeds (Yin et al., 2018), respectively. Similarly, regulatory networks of lipid and FA metabolism were revealed in Carya cathayensis Sarg (Huang et al., 2016(Huang et al., , 2021 and oil palm (Zheng et al., 2019) by transcriptomic and miRNA sequencing, respectively.
From the previous studies, it can be speculated that the integration of transcriptomic and lipidomic studies can provide a comprehensive view of the dynamic processes regulating the seed development and oil biosynthesis in various plants (Dossa et al., 2017;Wang et al., 2019;Chen et al., 2020). However, no such data have been collected for sesame seeds to date. Although there is sparse information on the chemical compositions (lignan, sesamin, and sesamolin), amino acid types and contents, and FA properties of sesame seeds from different regions, so far no transcriptome and degradome analyzes have been performed (Wei L. B. et al., 2016;Gacek et al., 2017;Sui et al., 2018;Ahmed et al., 2020). The specific functional genes and key metabolic pathways that regulate the synthesis of UFA in sesame seeds in particular are still unknown and unexplored. Therefore, a more comprehensive study on the controlling roles of miRNAs and their targets in lipid and FA biosynthesis is needed to further elucidate the underlying mechanism in the differential seed development of sesame.
Herein, the FA compositions in developing sesame seeds were analyzed, and the first-hand data on the combination of transcriptomics, small RNAs, and degradome in sesame seeds harvested in different developmental periods were generated and integrated to develop a comprehensive resource. We anticipate that this study will provide valuable information on the mechanism of lipid and FA biosynthesis in sesame seeds in different developmental periods. Furthermore, this study would promote molecular breeding research to improve the valueadding effects of sesame oil and the resulting functional plant products. Moreover, the collected miRNA expression patterns and extensive miRNA-mRNA regulatory network data can also be used for the upcoming expansion of plant oils with enriched UFA content.

Chemicals, Plant Materials, and Total RNA Extraction
For this study, sesame seeds named Wan Zhi No.2 (WZ2) were collected from the experimental field of the Crop Research Institute of Anhui Academy of Agricultural Sciences, Hefei, China in 2019. The whole growth period of WZ2 was about 88 days, and its seed yield was recorded at 1,788 kg/ha with a crude fat content of 55% (Wang et al., 2011). It was bred from sexual hybridization and systematic breeding of a local variety "Xiaozibai" (female parent) and Yuzhi No.4 (male parent). Beginning at 7 days after flowering (DAF), three biological replicates of sesame seeds for each sampling development period (7, 14, 21, and 28 DAF in August 2019) were handpicked until fully matured and then stored in liquid nitrogen at −80 • C until further use.
For transcriptomic analysis, miRNA sequencing, and degradome analysis, total RNA of each sample was extracted from the collected seeds, including three biological replicates, using Trizol reagent (Invitrogen, Waltham, MA, United States) according to the previous method . a rotary evaporator (Dbrowski et al., 2020). The FA composition of the sesame oil was analyzed by gas chromatography-mass spectrometry (Agilent GC-MS 7890A, Agilent Technologies, Santa Clara, CA, United States) equipped with a DB-5MS column and NIST 08 (Gaithersburg, MD, United States) according to the procedure described (Zhu et al., 2021). Crude fat 50 mg was weighed into a 10-ml centrifuge tube followed by the addition of 2 ml 1% H 2 SO 4 -CH 3 OH reagent and heating for 30 min at 70 • C in a water bath. Then, 2 ml of N-hexane was added, and water was filled to the top of the centrifuge tube (10 ml). The whole mixture was shaken well and allowed to stand for 24 h. After the centrifugation, supernatant was analyzed by GC-MS. GC conditions were as follow: DB-5 column (30 mm × 0.32 mm, 0.25 µm), the carrier gas was high purity helium, the injection temperature was 250 • C, the column flow rate was 1.25 ml/min, and the micro injector was used to inject 1.0 µl. The temperature rise program was 50 • C (2 min), 5 • C/min to 270 • C. Mass spectrometry conditions were as follow: EI ion source, ion source temperature 220 • C, GC-MS interface temperature 250 • C, electron energy 70 eV, scan mass range 60-500 m/z.
According to the results of the GC-MS analysis, the percentage of each FA component was calculated by the peak area normalization method (i.e., the percentage of each peak area to the total peak area).

Transcriptome Sequencing and Assembly Analysis
A cDNA library constructed from the pooled RNA from the different samples of sesame seeds was sequenced using the Illumina 6000 sequencing platform according to Illumina pairedend RNA-seq. Prior to assembly, the low-quality reads were omitted, and overall cleaned paired-end reads were generated. The raw sequence data were submitted to the NCBI database, and the reads of the samples were mapped to the UCSC 1 reference genome 2 using the HISAT package. Then, all transcriptomes were generated using StringTie to examine the expression level for FPKM after the final transcriptome was constructed. Here, the differentially expressed genes (DEGs) and mRNAs with p value of less than 0.05 and a log 2 | fold change| of more than 1 were screened with the R package (Lyu et al., 2020).

Small RNA Sequencing and miRNA Identification
In this assay, the 12 small RNA libraries (three biological replicates) of the four different period samples were constructed and processed by screening with a high-throughput sequencing method (Illumina HiSeq 2500 platform) at LC-BIO (Hangzhou, China). Raw RNA reads were generated and estimated using Illumina FastQC to obtain Q30 data. The raw reads were subjected to an in-house program, ACGT101-miR (LC Sciences, Houston, TX, United States) to remove adapter dimers, junk data, low complexity reads, common RNA families (rRNA, tRNA, snRNA, and snoRNA) and repeats, and the valid reads were 1 http://genome.ucsc.edu/ 2 https://www.ncbi.nlm.nih.gov/genome/?term=Sesamum+indicum+L obtained. Four samples of small RNA sequencing were analyzed based on the previous bioinformatics method (Yuan et al., 2013;Bisgin et al., 2018). Read distribution was calculated for different libraries and normalized for miRNA prediction using global normalization methods (Lei et al., 2021;Zhang et al., 2021).
The differentially expressed miRNAs in different seed development periods were evaluated using the DESeq R package (1.8.3) and found to be significant at P < 0.05. Then, to identify the differentially expressed known and novel miRNAs in sesame seeds, unique sequences with length of 18-25 nt were selected to map the specific species precursors in miRBase 3 by BLAST search. The unique sequences that can be further compared with the genome were considered as known miRNAs. The unique sequences that can be compared with the genome but cannot be compared with the pre-miRNAs of the species selected in miRBase were assumed to be novel miRNA candidates.

Degradome Analysis and Target Identification
The RNA libraries containing three biological replicates from four samples of the developmental period of sesame seeds were mixed as one degradome library. Then, the four small RNA samples and one degradome composite sample were analyzed on the sequencing platform (IlluminaHiSeq 2500) at LC-Bio Co., Ltd., (Hangzhou, China), and the degradome libraries were constructed according to the description of the manufacturer. Furthermore, the degradome reads were aligned with the transcriptome data of the sesame samples. Finally, all the identified known and novel miRNA-target genes were mapped to the database (see text footnote 2). Analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was performed using Agrigo and the Perl script, respectively, and the non-redundant GO enrichment terms were visualized using Revigo. In this study, a miRNA target gene network was constructed using Cytoscape (version 3.4.0).

Analysis of Differentially Expressed Target Genes
To obtain the target genes of lipid and FA synthesis, 12 RNA libraries were constructed from the sesame seed samples collected in four different developmental periods. For each library, all sequences were mapped to the CDS, and the number of each gene was counted using standard parameters. All the clean reads were mapped to the assembled unigenes of S. indicum L. for annotation. The RPKM of each gene referred to the number of sequencing fragments mapped to the CDS of the gene and was calculated by one thousand transcripts per million sequencing bases. According to the similarity of the gene expression profiles of samples, the heat map of the clustering expression pattern was constructed using log 10 (FPKM + 1) to represent the miRNAs and gene expression. The standard for measuring the DEGs among the different samples was based on P-value < 0.05 and fold change greater than 2.

Co-expression Analysis
For further exploration of the relationship between each module and FA synthesis, DEGs and expression levels of each module were screened, and the correlation between each module and nine types of FA synthesis was tested. Weighted co-expression network analysis was performed using the three elements (node, edge line, and network) of Cytoscape network graph (LC-Bio Co., Ltd., Hangzhou, China). As a result of the data collection and import with the Cytoscape package, the basic framework of the network graph was created after confirmation. All the other parameters were used with modification and adjustment. The co-expression networks for selected miRNAs and target pairs were visualized using Cytoscape (version 3.4.0).

Analysis of Gene Expression and Validation by qRT-PCR
For data validation, we selected 12 differentially expressed miRNA-target pairs to perform qRT-PCR to verify the highthroughput sequencing results. SYBRPrimeScript TM miRNA RT-PCR Kit (TaKaRa, Dalian, China) was used for reverse transcription reactions for miRNAs, and PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) was used for reverse transcription of cDNAs. U6 and beta-tubulin (TUB) were selected for normalization of miRNAs and target genes, respectively (Zhang et al., 2020). The names of miRNAs and their target genes, sequences, and primers in RT-qPCR test experiment are given in Supplementary File 1. The 2 − Ct method was used to calculate the relative fold change of miRNA target genes, and the PCR reaction of each sample was repeated three times.

Lipid Accumulation in Sesame Seed Development
In this study, we harvested sesame seed samples of the four developmental periods to evaluate their FA composition by gas chromatography-mass spectrometry (GC-MS) (Figure 1). The results show five dominant components, namely, oleic acid (C18:1, 43.61% of total FAs at S1), linoleic acid (C18:2, 46.77% of total FAs at S4), palmitic acid (C16:0, 20.96% of total FAs at S3), cis-7-hexadecenoic acid (C16:1, 14.72% of total FAs at S2), and stearic acid (C18:0, 19.7% of total FAs at S3) as shown in Figure 1B. It is found that all these five FAs had higher percentage and dominated over the other FAs (more than 89% of the total FA) during the four consecutive periods and reached 98.8% of the total FA in the S4 periods. Based on this, the development of sesame seeds can be divided into three distinct periods; starting with a low content of oil and FAs in the initial stage (S1-S2), followed by a rapid accumulation of oil in the second and third stages (S3), and a gradual increase in oil and FAs until full maturation (S4) (Figures 1B,C). However, saturated fatty acids (SFAs) showed a slight decrease, and UFAs showed an increase at S4; in particular, oleic acid, linoleic acid, and hexadecenoic acid was the dominant factors to this significant increase ( Figure 1C). Overall, oleic acid and linoleic acid were maintained at relatively high levels during sesame seed development, and UFAs dominated total FAs with relatively high proportions ( Figure 1C).

Transcriptome Sequence Analysis in Sesame Seeds Under Different Growing Time Points
RNA pools were obtained from sesame seed samples harvested at four different developmental time points to construct the transcriptome library. After refining the raw sequence data, only 147,561,222 reads remained valid out of 213,772,171 raw reads. Table 1 shows the annotation of sequence alignments uniq_miRNAs with the NCBI non-redundant protein (Nr) database. In addition, the results of the BLASTX algorithm showed the presence of 6,292 uniq_miRNAs with an average length of 70b (Table 1). Twelve RNA sequencing data with three biological replicates from the four different developmental periods were accumulated for the identification of DEGs with the aim of gene expression profile.

Small RNA Sequencing Profile
Small RNA sequencing for 12 libraries was performed to decipher miRNA regulated post-transcriptional changes associated with lipid and FA biosynthesis. Figure 2A depicts the size distributions and similar patterns of unique and valid 18-to 25-nt-long sequences obtained after applying the different steps of data processing. Among the 12 libraries, the 24-nt (49.1%) RNAs were found to be most abundant and diverse, followed by 23% small RNAs with 21-nt-long sequences (Figure 2A and Supplementary File 2).
Based on the abundance in miRNA expression levels, 1,025 miRNAs were obtained and categorized into five different groups after database and sequence reading analysis (Figure 2C and Supplementary File 3). Among the unique mature miRNAs and pre-miRNAs, 354 unique mature miRNAs were categorized into groups 1, −2, and −3 (gp1, −2, and −3) (Supplementary File  3). In addition, the miRNAs from gp4 were not registered in miRBase, and 671 unique mature miRNAs corresponding to 664 pre-miRNAs were considered as new candidate miRNAs (Supplementary File 3). Interestingly, the pattern of the proportion of mature miRNAs in each sampling developmental period showed some similarity (Figure 2A). Moreover, a total of 220 miRNAs (114 known miRNAs and 106 novel miRNAs) were expressed in all four developmental periods, as shown in Figure 2B. In addition, the sum of the first nucleotide bases of miRNA was different between the known and novel miRNAs. Adenosine (A, 49.82%) was reported as the most abundant nucleotide, followed by uracil (U, 20.16%) and cytosine (C, 18.33%) (Supplementary File 4).

Differential Expression of miRNAs in Sesame Seeds Under Different Growing Time Points
We used the high-throughput sequencing method to analyze and compare the differential expression of miRNAs (DEmiRNAs) related to lipid and FA biosynthesis in four samples of  sesame seeds. According to the analysis, among the 220 DEmiRNAs, 126 miRNAs were up-regulated and 74 miRNAs were significantly down-regulated in S3 vs. S2 (Figure 3). Moreover, the distribution patterns showed significant upregulation of 113 miRNAs and down-regulation of 63 in S3 vs. S1 (Figure 3). These data suggest that stage 3 may be the critical control point that strongly influences miRNA expression levels during lipid and FA biosynthesis. According to the heat map of potential DEmiRNAs, 220 miRNAs were found to be significantly expressed across the four periods with similar distribution patterns of high or low expression levels verified after clustering analysis (Supplementary Figure 1).

Prediction of Known and Novel miRNAs
To further understand the possible function and the regulatory patterns of the identified miRNAs and their target genes responsible for lipid and FA biosynthesis during sesame seed  File 6). Based on the in-depth analysis, the miRNA (aly-MIR408-p3_2ss18CT19TG) was identified with 520 targets, which is the highest number of transcripts cleavages by the same miRNA (Supplementary File 6). To summarize, 12 miRNAs were associated with the regulation of lipid and FA biosynthesis during the developmental periods of sesame seed. Among the 12 miRNAs, 11 could cleave only one transcript

GO and KEGG Pathway Analyzes of Targets
For the manifestation of the potential biological role of 2,296 miRNAs in sesame seed development, 2,186 target genes were categorized according to GO enrichment (Supplementary Figure 2). The enrichment data showed that most of the target genes were associated with biological processes involving "regulation of transcription, " "oxidation-reduction process, " and "protein phosphorylation" (Figure 5A). The GO terms "nucleus, " "cytoplasm, " and "plasma membrane" in the "cellular component" category accumulated the most frequent targets (Figure 5A). The GO terms "protein binding" and "molecular function" in the molecular function category had the highest number of targets compared with the other terms ( Figure 5A). There, interesting GO terms, such as "embryo development ending in seed dormancy" and "abscisic acid-activated signaling pathway" were revealed ( Figure 5A). For pathway analysis, KEGG annotations were performed in the top 20 pathways according to the KEGG annotated gene number. Most target genes were found enriched in "ribosome" and "protein processing in endoplasmic reticulum, " followed by "glycerolipid metabolism" and "glycerophospholipid metabolism" (Figure 5B). In addition, more genes were involved in "biosynthesis of secondary metabolites" consisting of "glycosylphosphatidylinositol (GPI)-anchor biosynthesis" and "biosynthesis of UFAs" (Figure 5B). Based on the enrichment data, "fatty acid degradation, " "beta-alanine, " and "inositol phosphate metabolism" were shown to be the most significantly enriched pathways (FDR 0.05) (Figure 5B).

miRNAs and Target Genes in the Context of Lipid and Fatty Acid Biosynthesis
The transcriptomics investigation revealed that a total of 220 DEmiRNAs corresponding to 3,705 target genes were also differentially expressed. After miRNA-target pair analysis, 32 miRNAs with 33 target genes were found to be associated  File 7). In this analysis, most DEmiRNAs and their targets showed a similar expression pattern as in high-throughput sequencing (Figure 6). According to the experimental qRT-PCR data, eight known miRNAs and their target pairs showed negative expression correlation (Figure 7). For example, miRNA (syl-miR166c-5p-1ss21AG) showed higher expression in S1 and S2, which decreased in S3 and S4, while target FAD2 showed lower expression in S1 and S2, which increased in S3 and S4. Moreover, four novel miRNA-target pairs (PC -3p-15838_568/LOC105162096, PC -5p-145800_91/LOC105178588, PC -3p-62_86681/LOC105162196, and PC -3p-14226_618/LOC105163373) also showed inverse association at the expression level, which control transcriptional repression by targets and their corresponding miRNAs (Figure 7). In addition, we detected three miRNAs-target pairs (sly-miR166c-5p_1ss21AG/FAD2, stu-MIR8005c-p3_1ss13AG/LOC105159459, and gma-MIR169o-p3_2ss16AC21CG/LOC105161564) that exhibited not only a relationship at the high expression level but also an apparent negative relationship. Interestingly, the miRNA bol-MIR9410-p3_2ss4TG17TA was found to have two target genes (LOC105176298 and LOC105167408). We hypothesize that the functional analysis and cloning of the three genes FAD2, LOC105159459, and LOC105161564 could potentially facilitate the deciphering of the mechanistic details associated with the biosynthesis of UFAs during sesame seed development.

Gene Co-expression Network Associated With Lipid and Fatty Acid Biosynthesis
To elucidate the function of miRNA-target pairs in the context of lipid and FA biosynthesis regulation in different developmental periods, we constructed a coexpression network with 33 genes from nine categories, mainly associated with glycerolipid metabolism, glycerophospholipid metabolism, UFAs, FA biosynthesis, FA elongation, FA degradation, steroid synthesis, arachidonic acid, and ether ester metabolism (Figure 8).
After further analysis of the functional assignments, seven and five genes were found to be responsible for the categories of glycerolipid metabolism and glycerophospholipid metabolism, respectively, among the 33 genes (Supplementary File 7). According to gene co-expression analysis, 13 genes were associated with FA biosynthesis and metabolism, of which three genes (FAD2, LOC105174578, and LOC105169658) regulated the biosynthesis of UFAs (Figure 8). In addition, four miRNAs (osa-miR167a-5p_2ss2GA21AG, osa-miR167a-5p_2ss14AG21A, osa-miR167a-5p_2ss2GA21AT, and osa-miR167a-5p) had the same target gene LOC10514358 involved in FA elongation (Figure 8). In the gene co-expression network, we found that miRNA (bol-MIR9410-p3_2ss4TG17TA) could simultaneously regulate three target genes (LOC105176298, LOC105167408, and LOC105158444) involved in glycerophospholipid metabolism, FA biosynthesis, and FA degradation, respectively (Figure 8). Similarly, some important miRNAs and their target genes involved in glyceride metabolism, ether ester metabolism, arachidonic acid, and steroid synthesis have been linked. Notably, the data suggest that the coexpression network of miRNA and target genes can be used in further studies to identify the function of known and novel genes related to lipid and FA biosynthesis in sesame seed development.

DISCUSSION
To improve the varietal characteristics of sesame, increasing the UFA content (especially oil yield and FA composition) is the first and most important priority for the development of sesame breeding (Wang D.D. et al., 2018;Ahmed et al., 2020). With the possibility of the presence of the sesame genome and the advent of high-throughput sequencing for small RNA analysis in sesame seeds, we hereby elaborated the scope of transcriptomic, miRNA, degradation group analysis data, and gene network data during sesame seed development (Zheng et al., 2019).
In this study, 32 miRNA and target pairs involved in lipid and FA biosynthesis were screened out. Importantly, most of these target pairs were related to glyceride metabolism and FA biosynthesis, and five miRNA target gene pairs related to glycerol phospholipid metabolism were also identified by degradome sequencing data. Interestingly, seven target genes involved in FA degradation were found, but three of them did not have miRNA pairs. We further analyzed and found that the expression levels of four miRNA genes changed significantly with the maturity of sesame, such as, stu-MIR8005c-p3_1ss13AG/LOC105159459 involved in glyceride metabolism, sly-miR166c-5p_1ss21AG/FAD2 involved in UFA metabolism, gma-MIR169o-p3_2ss16AC21CG/LOC105161564 involved in FA biosynthesis, and PC-3p-62_86681/LOC105162196 involved in steroid synthesis. Thus, these four miRNA-target pairs were identified as potential regulators of the lipid and FA biosynthesis process during sesame seed development. To the knowledge of the authors, this is the first study ever to report novel miRNAs FIGURE 4 | Some differentially expressed miRNA t-plots and cleavage sites identified by degradome sequencing.
targeting miRNA-target pairs that had not been previously discovered, for example, PC-3p-14226_618, PC-5p-228878_56, PC-5p-145800_91, and PC-3p-59282_203. The results suggest that these specific miRNA target pairs may be involved in the regulatory network of lipid and FA synthesis, and play a central role in controlling gene expression and metabolic pathways throughout development periods.
Previous studies have shown that oil synthesis regulates de novo FA biosynthesis in the plastid and TAG assembly in the endoplasmic reticulum (ER), and that oil body formation (Reynolds et al., 2019) in oil plants is represented (Zheng et al., 2019). This study signified the presence of 106 unigenes involved in oil biosynthesis pathways (Figure 9). In general, lipid and FA syntheses are initiated with the supply of acetyl-CoA substrate (Wang et al., 2019), which is first converted to malonyl-CoA (Chen et al., 2020), and carbon flux is initiated via FA synthase and most other enzymes associated with lipid and FA biosynthesis (Huang, 1996). In this study, the initiation and acyl chain elongation of de novo FA biosynthesis were reported to be possibly regulated by 50 unigenes, and the latter showed apparent down-regulation of at least one isoform in S4 compared with S1, suggesting their functional involvement in the initiation of FA biosynthesis (Figure 9). In this study, ACCase catalyzed the first reaction involving four genes (conversion of acetyl-CoA to malonyl-CoA) involved in the biosynthesis of FA (Figure 9). Furthermore, FATA is highly specific for 18:1-ACP, and FATB is specific for acyl-A, which are known to regulate chain termination during FA synthesis (Lu et al., 2009). Unigenes LOC105158149 and LOC105160370, which encode FATB, were up-regulated more than 10-fold during the S1 to S4 phase. In addition, four unigenes encoding LACS generating the acyl-CoA pool were detected, all of which were significantly down-regulated (Figure 9). The free long-chain FAs (16:0, 18:0, and 20:0) are esterified by LACS and output to  ER for further acyl processing (Bates et al., 2009). Pathways such as glycerolipid metabolism and ether-lipid metabolism are known to be linked by the compound, which converts to glycerone-3P (G3P) under the catalyzation of NAD + and FAD2 (Hu et al., 2009;Dehghan and Yarizade, 2014;Yasemin et al., 2018;Macovei et al., 2021). 49 unigenes related to the formation of TAG were identified, e.g., long-chain acyl-CoA and G3P, which could start the assemble of TAG; diacylglycerol acyltransferase (DGAT), which was the major enzymes that catalyze diacylglycerol (DAG); and phosphatidylcholine (an acyl donor), which could combine with oleosin to form an oil body (Kennedy, 1961;Peng and Weselake, 2011). Most of the unigenes encoding DAG were reported with high transcript levels, suggesting that they are an important signaling pathway during sesame seed development (Wang R.K. et al., 2018;Huai et al., 2020). On the other hand, relatively active unigenes, encoding acyl editing, PE-DAG, and PC-DAG interconversion were likely responsible for the incorporation of UFAs into TAG in sesame seeds (Peng and Weselake, 2011;Wang R.K. et al., 2018;Yang et al., 2018;Tyagi et al., 2019). The data emphasized the  Frontiers in Plant Science | www.frontiersin.org FIGURE 9 | Relation network of miRNAs and target genes with lipid and fatty acid biosynthesis in sesame seeds. The blue color represents the down-regulated gene expression pattern, while the red color represents the up-regulated gene expression pattern. essential roles of FAD2, FabF/KASII, LACS, ACO-E, LPCAT, and some unknown genes during the biosynthesis of FA and other UFAs during sesame seed development (Figure 9).
Overall, the regulatory co-expression network construction data showed that the nine kinds of oil synthesis related pathways were the most likely regulatory factors of lipid and FA biosynthesis; and so far, it can be stated that miRNA regulates the complex network of genes and their target genes involved in various physiological and biological functions such as regulation of sesame responses to oil synthesis. These target pairs gma-MIR1510a-p3_2ss12AT19CT/LOC105159459, sly-miR166c-5p_1ss21AG/FAD2, gma-MIR169o-p3_2ss16AC21CG/LOC10516 1564, and PC-3p-62_86681/LOC105162196 were found to be potentially relevant to oil synthesis in developing sesame seeds. Nevertheless, future investigation of these four pairs is warranted to reveal their precise molecular functions and regulatory mechanisms. Herein, complementary comparative multi-omics data were generated in four successive developmental periods to unravel the important pathways responsible for the high-oil yield and increased UFA content in sesame seeds. The results collected would also fill the knowledge gap on the mechanism of oil deposition and FA saturation transformation during embryonic development in other oil plant seeds.

CONCLUSION
In conclusion, this study provides a comprehensive framework for a better understanding of the lipid and FA biosynthesis mechanism of sesame under different growth periods. The1,006 target genes of 13,801 miRNAs were extracted from mixed degradome sequencing, which accounted for 65 novel miRNAs targeting 8,051 transcripts, and 3,705 target genes of 220 miRNAs were also differentially manifested. At least 12 potential lipidrelated miRNAs and target genes were analyzed by RNA-seq and qRT-PCR. Perfect harmonization of high level of FAD2 in the mature period promoted the increase of linoleic acid content in sesame. Thirty-three co-expressed genes formed a co-regulatory sub-network, which showed that some important miRNAs and their target genes were involved in the pathways of lipid and FA biosynthesis such as glyceride metabolism, ether ester metabolism, arachidonic acid, and steroid synthesis, and the pathways were also found to be interconnected. In addition, the biosynthetic and regulatory genes identified from the genomewide co-expression network may expand the understanding of lipid and FA biosynthesis in other plant seeds and increase the functional value of the resulting food products.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/bioproject/?term=PRJNA739094.

AUTHOR CONTRIBUTIONS
Y-PZ: conceptualization, methodology, software, and writingoriginal draft preparation. Y-YZ: data curation and validation. KT: data curation and investigation. FZ and J-GZ: methodology and software. FH: visualization and investigation. P-CW: visualization, supervision, and investigation. Z-JW: funding acquisition, supervision, and writing-reviewing and editing. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the Major Projects of Science and Technology in Anhui Province (201903a06020021, 201904a06020008, 202004a06020042, and 202004a06020052) and the National Natural Science Foundation of China (31850410476).