Science, Irving K. Barber Faculty of (Okanagan)
Non UBC
Biology, Department of (Okanagan)
DSpace
Gao, Yunfei
Wang, Le
Cao, Shanshan
Chen, Liangyu
Li, Xueying
Cong, Weixuan
Yang, Songnan
Zhang, Jian (Professor of biology)
Nie, Xiaojun
Zhang, Jun
Tiger nut (Cyperus esculentus L.) is recognized for its high oil and oleic acid content in underground tubers. However, the molecular mechanisms governing growth, development, and fatty acid accumulation in these tubers are not well understood. This study employed gas chromatography–mass spectrometry (GC–MS) and small RNA sequencing on tiger nut tubers across five developmental stages. The findings indicate that the critical period for accumulating dry matter and oils, particularly oleic acid, occurs between 35 and 75 days after tuber formation. A total of 183 microRNAs (miRNAs) were identified, comprising 31 known and 152 novel miRNAs. Approximately half of these miRNAs (such as ces-miR156b-3p and ces-miR166a-3p) exhibited differential expression during and around the key periods of metabolite synthesis. The predicted target genes of these miRNAs were significantly enriched in glycerate 3-phosphate metabolism and cell growth processes. Furthermore, 13 miRNA–mRNA interaction modules related to oil accumulation and tuber growth were identified, and these target genes’ expression levels showed significant differences during the tuber developmental stages. These findings advance the understanding of the molecular mechanisms underlying tuber development and oil accumulation in tiger nut.
Citation: Gao, Y.; Wang, L.; Cao, S.;Chen, L.; Li, X.; Cong, W.; Yang, S.;Zhang, J.; Nie, X.; Zhang, J.Identification of miRNAs Involved inLipid Metabolism and TuberDevelopment in Cyperus esculentus L.Plants 2024, 13, 3305. https://doi.org/10.3390/plants13233305Academic Editor: AlexandraS. DubrovinaReceived: 7 October 2024Revised: 13 November 2024Accepted: 20 November 2024Published: 25 November 2024Copyright: © 2024 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).ArticleIdentification of miRNAs Involved in Lipid Metabolism andTuber Development in Cyperus esculentus L.Yunfei Gao 1,†, Le Wang 2,† , Shanshan Cao 1, Liangyu Chen 1 , Xueying Li 1, Weixuan Cong 1 , Songnan Yang 1 ,Jian Zhang 1,3, Xiaojun Nie 4,* and Jun Zhang 1,5,*1 College of Agronomy, Jilin Agricultural University, Changchun 130118, China;20210055@mails.jlau.edu.cn (Y.G.); zhangjian@jlau.edu.cn (J.Z.)2 College of Plant Science, Jilin University, Changchun 130062, China; lewang27@jlu.edu.cn3 Department of Biology, University of British Columbia, Okanagan, Kelowna, BC V1V 1V7, Canada4 State Key Laboratory for Crop Stress Resistance and High-Efficiency Production, College of Agronomy,Northwest A&F University, Xianyang 712100, China5 National Crop Variety Approval and Characteristic Identification Station, Jilin Agricultural University,Changchun 130118, China* Correspondence: small@nwsuaf.edu.cn (X.N.); zhangjun@jlau.edu.cn (J.Z.)† These authors contributed equally to this work.Abstract: Tiger nut (Cyperus esculentus L.) is recognized for its high oil and oleic acid content inunderground tubers. However, the molecular mechanisms governing growth, development, and fattyacid accumulation in these tubers are not well understood. This study employed gas chromatography–mass spectrometry (GC–MS) and small RNA sequencing on tiger nut tubers across five developmentalstages. The findings indicate that the critical period for accumulating dry matter and oils, particularlyoleic acid, occurs between 35 and 75 days after tuber formation. A total of 183 microRNAs (miRNAs)were identified, comprising 31 known and 152 novel miRNAs. Approximately half of these miRNAs(such as ces-miR156b-3p and ces-miR166a-3p) exhibited differential expression during and around thekey periods of metabolite synthesis. The predicted target genes of these miRNAs were significantlyenriched in glycerate 3-phosphate metabolism and cell growth processes. Furthermore, 13 miRNA–mRNA interaction modules related to oil accumulation and tuber growth were identified, andthese target genes’ expression levels showed significant differences during the tuber developmentalstages. These findings advance the understanding of the molecular mechanisms underlying tuberdevelopment and oil accumulation in tiger nut.Keywords: Cyperus esculentus; miRNAs; GC-MS; oil synthesis; tuber development1. IntroductionTiger nut (Cyperus esculentus L.), a C4 plant in the Cyperaceae family, primarily yieldstubers as its main agricultural product [1,2]. These tubers accumulate substantial amountsof starch (23.2–29.9%) and oil (20.1–34.5%), and contain bioactive compounds includingorganic acids, alkaloids, and phenolics [3,4]. The tuber oil is predominantly composed ofmonounsaturated fatty acids, such as oleic acid, which impart exceptional oxidative stabilityand physicochemical properties [5,6]. This composition makes it suitable as a healthy edibleoil and a renewable raw material for biodiesel and other lipid-based chemical industries [7].The estimated tuber yield of tiger nut can reach up to 12 t/ha, which can produce anoil output of 1.5 t/ha, significantly surpassing conventional oil crops like soybean andpeanut [5]. Additionally, tiger nut demonstrates excellent soil adaptability, strong resistanceto poor soil conditions, low input costs, and resilience to pests and diseases [8]. Therefore,research on the development of tiger nut tubers and oil synthesis mechanisms will not onlyaid its genetic improvement but also elucidate the molecular mechanisms of oil productionin plants.Plants 2024, 13, 3305. https://doi.org/10.3390/plants13233305 https://www.mdpi.com/journal/plantsPlants 2024, 13, 3305 2 of 15Research indicates that tiger nut tubers exhibit rapid exponential growth in theirearly stages, transitioning to metabolite accumulation once they reach a certain size [9].While the timing of new tuber formation varies after sowing due to the different cultivars,the phase of rapid oil accumulation typically occurs 20 to 70 days post-formation [6,9,10].There is a significant accumulation of lipid droplets during tuber development, with lipidcomposition analysis identifying triacylglycerol (TAG) as the primary contributor to theincreased tuber oil content [6]. Transcriptomes across various developmental stages revealthat differentially expressed genes in the early stages are highly enriched in glycerolipid andlinoleic acid synthesis pathways [6]. Key candidate enzymes involved in the oil synthesisof encoding genes, such as acyl-ACP thioesterase (FATA) and diacylglycerol acyltransferase(DGAT), are activated early and sustain high expression levels until tuber maturation. Incontrast, stearoyl–ACP desaturase (SAD) and phosphatidate phosphohydrolase (PAH) areactivated and show expression during the mid- to late developmental stages. Comparativetranscriptomic studies with the related species Cyperus rotundus have revealed that theexpression levels of candidate genes involved in oil synthesis are significantly higher intiger nut; for example, FATA is 10.3 times higher [8]. Furthermore, the transcription factorsABSCISIC ACID INSENSITIVE3 and WRINKLED 1 are crucial in regulating the expressionof these oil synthesis-related genes. In addition, the introduction of exogenous genesCePAH1, CeDGAT2, and CeSAD into yeast and Arabidopsis significantly increases TAG andtotal oil content compared to the wild type [6,7,10,11]. However, the complete mechanismsunderlying tuber growth and fatty acid synthesis in tiger nut remain to be elucidated.Further investigation is needed to determine the potential existence of other key genesinvolved in oil synthesis and to understand how the expression levels of the identifiedgenes are regulated.MicroRNAs are short, endogenous, noncoding RNAs of 18–25 nucleotides (nt) inlength found in plants and animals. They regulate the expression levels of their targetgenes, which are involved in specific biological processes [12–15]. Numerous miRNAshave been identified in the seeds of oil crops via high-throughput small RNA [16–19].For example, 27 miRNAs regulate the gene expression related to seed development andoil accumulation in sea buckthorn [20]. In rapeseed, 32 miRNAs are directly involvedin fatty acid synthesis in seeds [16]. Moreover, some miRNAs have been identified asparticipating directly in organ development or oil synthesis. In rapeseed, miR166f regulatessilique length and harvest index [21]. Modified csa-miR159a increases oleic acid content inCamelina sativa [22]. In potatoes, the overexpression of miR156 can reduce tuber numberand size [14]. Despite the crucial roles of miRNAs in plant growth and oil accumulation,there are currently no reports on miRNAs in tiger nut.To elucidate the role of miRNAs in tuber development and fatty acid synthesis intiger nut, we employed GC-MS and small RNA sequencing to obtain dynamic informationon fatty acid composition and miRNA changes across five stages of tuber development.Then, bioinformatics tools were used to identify miRNA target genes and their biologicalannotations, aiming to associate them with tuber traits. Ultimately, this research offersnovel insights into the miRNAs and molecular mechanisms underlying tuber developmentand lipid accumulation in tiger nut.2. Results2.1. Changes in Tiger Nut Tuber Weight, Oil Content and Fatty Acid Composition at FiveDevelopment PeriodsIn tiger nut cultivar JINONGSHA6, tuber enlargement becomes distinctly observ-able 45 days after sowing (DAS). As growth progresses, the tuber color transitions fromwhite (S1) to dark brown (S5) (Figure 1A). Tuber weight and oil content both accumulateprogressively from S1 to S5 (Figure 1B). A quick increase in the dry weight of the tuberswas observed, with the peak growth rate occurring between stages S3 and S4 (85 DAS to105 DAS), at 0.02 g/day per tuber. The second highest growth rate was observed betweenstages S2 and S3 (65 DAS to 85 DAS), with an average increase of 0.014 g/day per tuberPlants 2024, 13, 3305 3 of 15(Figure 1B). Meanwhile, the fastest rate of oil accumulation occurred between stages S2and S3 at 6.06 mg/day per tuber, followed by a slow rate of 1.97 mg/day per tuber fromstages S3 to S4. These observations indicate that stages S2 to S4 are critical periods for theaccumulation of dry matter and oil content in tiger nut.Plants 2024, 13, x FOR PEER REVIEW 3 of 15 observed, with the peak growth rate occurring between stages S3 and S4 (85 DAS to 105 DAS), at 0.02 g/day per tuber. The second highest growth rate was observed between stages S2 and S3 (65 DAS to 85 DAS), with an average increase of 0.014 g/day per tuber (Figure 1B). Meanwhile, the fastest rate of oil accumulation occurred between stages S2 and S3 at 6.06 mg/day per tuber, followed by a slow rate of 1.97 mg/day per tuber from stages S3 to S4. These observations indicate that stages S2 to S4 are critical periods for the accumulation of dry matter and oil content in tiger nut. A total of 14 fatty acids were identified in tiger nut tubers using GC–MS (Table S1). Oleic acid, linoleic acid and palmitic acid were the three most abundant fatty acids at all developmental stages, collectively accounting for over 93% of the total fatty acid content, and are thus considered the primary fatty acids in tiger nut tubers (Figures 1C,D and S1). Furthermore, six unsaturated fatty acids (oleic acid, linoleic acid, palmitoleic acid, myristoleic acid, eicosapentaenoic acid, and vaccenic acid) constituted 75% of the total fatty acid content (Figure S1). Oleic acid comprised over 57% of the total fatty acid content except during the early developmental stages S1 and S2. Meanwhile, the composition ratios of fatty acids changed with developmental stages (Figure 1C). Consistent with the trends in tuber weight and total oil content, the content of most fatty acids accumulated continuously during tuber growth (S4 to S5). Notably, the stages from S2 to S4 were marked by the rapid accumulation of various fatty acids (Figure 1D). For example, the accumulation rates of oleic acid were 4.88 mg/day and 1.49 mg/day during the S2 to S3 and S3 to S4 stages, respectively. In summary, these results indicate that the fatty acids in tiger nut tubers are primarily composed of unsaturated fatty acids, with oleic acid being the most critical component. Additionally, t e period from S2 to S4 represents a crucial stage for fatty acid synth sis in the tubers. Figure 1. Growth, oil accumulation, and fatty acid composition of tiger nut tubers. (A) Five developmental stages of tiger nut tubers. (B) Curves of dry matter and oil accumulation. (C) Proportion of the six most abundant fatty acids. (D) Changes in the accumulation of various fatty acids over the five stages. Figure 1. Growth, oil accumulation, and fatty acid composition of tiger nut tubers. (A) Five develop-mental stages of tiger nut tubers. (B) Curves of dry matter and oil accumulation. (C) Proportion ofthe six most abundant fatty acids. (D) Changes in the accumulation of various fatty acids over thefive stages.A total of 14 fatty acids were identified in tiger nut tubers using GC–MS (Table S1).Oleic acid, linoleic acid and palmitic acid were the three most abundant fatty acids at alldevelopmental stages, collectively accounting for over 93% of the total fatty acid content,and are thus considered the primary fatty acids in tiger nut tubers (Figures 1C,D and S1).Furthermore, six unsaturated fatty acids (oleic acid, linoleic acid, palmitoleic acid, myris-toleic acid, eicosapentaenoic acid, and vaccenic acid) constituted 75% of the total fatty acidcontent (Figure S1). Oleic acid comprised over 57% of the total fatty acid content exceptduring the early developmental stages S1 and S2. Meanwhile, the composition ratios offatty acids changed with developmental stages (Figure 1C). Consistent with the trends intuber weight and total oil content, the content of most fatty acids accumulated continuouslyduring tuber growth (S4 to S5). Notably, the stages from S2 to S4 were marked by the rapidaccumulation of various fatty acids (Figure 1D). For example, the accumulation rates ofoleic acid were 4.88 mg/day and 1.49 mg/day during the S2 to S3 and S3 to S4 stages,respectively. In summary, these results indicate that the fatty acids in tiger nut tubersare primarily composed of unsaturated fatty acids, with oleic acid being the most criticalcomponent. Additionally, the period from S2 to S4 represents a crucial stage for fatty acidsynthesis in the tubers.2.2. Overview of sRNA Sequencing and the Identification of miRNAsFifteen small RNA libraries were constructed across five developmental stages toidentify the miRNAs involved in material accumulation in tiger nut tubers. A total of178,164,492 raw reads were obtained from all libraries through high-throughput sequencing,with Q30 values exceeding 95% for each library (Table S2). After the quality controlPlants 2024, 13, 3305 4 of 15described in the Materials and Methods Section, each library retained between 9,601,609and 13,391,806 clean reads. Finally, over 94% of the clean reads from each library weremapped to the tiger nut genome. These results demonstrate that the quality of the smallRNA sequencing data was sufficient for further analysis.According to the miRNA analysis pipeline described in the Materials and Methods,we identified 183 mature miRNAs in 15 libraries from five developmental stages of tubers(Table S3). Among these, 31 miRNAs are known plant miRNAs belonging to 11 miRNAfamilies. The miR166 and miR168 families have eight and five members, respectively(Figure 2A). The newly identified miRNA members were named in the form of their ces-miRn-number. In total, 152 novel mature miRNAs, ranging in length from 19 nt to 27 nt,were discovered. Of these, 24.34% have a length of 21 nt, while 41.44% have a length of24 nt (Figure 2B). In addition, there were differences in the preferential distribution of thefirst nucleotide between known and novel miRNAs (Figure 2C,D). Known mature miRNAstypically begin with 5′-G and 5′-C, whereas novel mature miRNAs often start with 5′-Uand 5′-G.Plants 2024, 13, x FOR PEER REVIEW 4 of 15 2.2. Overview of sRNA Sequencing and the Identification of miRNAs Fifteen small RNA libraries were constructed across five developmental stages to identify the miRNAs involved in material accumulation in tiger nut tubers. A total of 178,164,492 raw reads were obtained from all libraries through high-throughput sequencing, with Q30 values exceeding 95% for each library (Table S2). After the quality control described in the Materials and Methods Section, each library retained between 9,601,609 and 13,391,806 clean reads. Finally, over 94% of the clean reads from each library were mapped to the tiger nut genome. These results demonstrate that the quality of the small RNA sequencing data was sufficient for further analysis. According to the miRNA analysis pipeline described in the Materials and Methods, we identified 183 mature miRNAs in 15 libraries from five developmental stages of tubers (Table S3). Among these, 31 miRNAs are known plant miRNAs belonging to 11 miRNA families. The miR166 and miR168 families have eight and five members, respectively (Figure 2A). The newly identified miRNA members were named in the form of their ces-miRn-number. In total, 152 novel mature miRNAs, ranging in length from 19 nt to 27 nt, were discovered. Of these, 24.34% have a length of 21 nt, while 41.44% have a length of 24 nt (Figure 2B). In addition, there were differences in the preferential distribution of the first nucleotide between known and novel miRNAs (Figure 2C,D). Known mature miRNAs typically begin with 5′-G and 5′-C, whereas novel mature miRNAs often start with 5′-U and 5′-G. Figure 2. Overview of miRNAs: (A) Number of families of known miRNAs; (B) Length distribution of miRNAs; (C) First nucleotide preference of known miRNAs; (D) First nucleotide preference of novel miRNAs. 2.3. Differential Expression and Functional Analysis of miRNAs at Different Growth Stages Differences in miRNA expression profiles across the developmental stages of tiger nut tubers may be associated with growth and metabolite accumulation. Expression levels were normalized using transcripts per million (TPM), with a threshold for reliable expression set at TPM > 1. The numbers of expressed miRNAs across the five stages were 146, 138, 141, 145, and 139, respectively, resulting in 169 unique miRNAs when combined (Table S4). For differential expression analysis, miRNAs reliably expressed in any two stages were compared, with differentially expressed miRNAs (DEmiRs) defined by an absolute value of fold change (FC) ≥ 1.5 and p < 0.05. The analysis revealed that there were fewer DEmiRs between adjacent stages than those between non-adjacent stages (Figure 3A). In comparisons of adjacent developmental stages, the highest number of DEmiRs was Figure 2. Overview of miRNAs: (A) Number of families of known miRNAs; (B) Length distributionof iR As; (C) First nucleotide preference of kno n iR As; ( ) First nucleotide preference ofnovel iR As.2.3. Differential Expression and Functional Analysis of miRNAs at Different Growth StagesDifferences in miRNA expression profiles across the developmental stages of tigernut tubers may be associated with growth and metabolite accumulation. Expressionlevels were normalized using transcripts per million (TPM), with a threshold for reliableexpression set at TPM > 1. The numbers of expressed miRNAs across the five stages were146, 138, 141, 145, and 139, respectively, resulting in 169 unique miRNAs when combined(Table S4). For differential expression analysis, miRNAs reliably expressed in any twostages were compared, with differentially expressed miRNAs (DEmiRs) defined by anabsolute value of fold change (FC) ≥ 1.5 and p < 0.05. The analysis revealed that there werefewer DEmiRs between adjacent stages than those between non-adjacent stages (Figure 3A).In comparisons of adjacent developmental stages, the highest number of DEmiRs wasobserved between S2 and S3, with a total of 40, comprising 23 that were upregulated and17 that were downregulated. This was followed by 26 DEmiRs between S4 and S5. Incontrast, fewer than 10 DEmiRs were identified between S1 and S2, as well as between S3and S4. Based on miRNA expression levels and differential expression analysis, 25 miRNAsexhibited differential expression between key periods and non-key periods of metaboliteaccumulation, and their expression level either increased or decreased in accordance withPlants 2024, 13, 3305 5 of 15developmental stages (Figure 3B). It is noteworthy that five known miRNAs from themiR166 family showed lower expression levels from S3 to S5 compared to S1 and S2, whiletwo members from the miR156 family displayed an opposite expression pattern to thatobserved in the miR166 family. Additionally, a novel miRNA, ces-miRn100, showed a veryreduced expression (TPM < 1) in S3 and S4, but its expression level in S5 was higher thanin S1 and S2, suggesting a potential close association with tuber metabolite synthesis.Plants 2024, 13, x FOR PEER REVIEW 5 of 15 observed between S2 and S3, with a total of 40, comprising 23 that were upregulated and 17 that were downregulated. This was followed by 26 DEmiRs between S4 and S5. In contrast, fewer than 10 DEmiRs were identified between S1 and S2, as well as between S3 and S4. Based on miRNA expression levels and differential expression analysis, 25 miRNAs exhibited differential expression between key periods and non-key periods of metabolite accumulation, and their expression level either increased or decreased in accordance with developmental stages (Figure 3B). It is noteworthy that five known miRNAs from the miR166 family showed l wer expression levels from S3 to S5 compared to S and S2, whil two m mb r from the miR156 family displayed opposite expression pattern to that observed in the miR166 family. Additionally, a novel miRNA, ces-miRn100, showed a very reduced expressio (TPM < 1) in S3 and S4, but its expression level in S5 was higher than in S1 and S2, suggesting a potential close association with tuber metabolite synthesis. Figure 3. Identification and functional analysis of DEmiRs. (A) Number of DEmiRs between any two stages. (B) miRNAs with expression changes before and after key stages of metabolite accumulation. (C) GO enrichment results of DEmiRs. The target genes of miRNAs were predicted to assess their potential functions, resulting in the identification of 618 target genes, 102 of which were annotated as transcription factors (Table S5). Gene ontology (GO) and KEGG enrichment analyses were conducted to elucidate the biological processes in which target genes of DemiRs were Figure . I tific ti functional analysis of DEmiRs. (A) Number of DEmiRs between a y twostages. (B) miRNAs with expression changes before and after key stages of metabolite accumulation.(C) GO enrichment results of DEmiRs.The genes of miRNAs were predicted to asses their potential functions, result-ing in the identification of 618 target genes, 102 of which were annotated as transcriptionf ctors (Table S5). Gene ontology (GO) and KEGG e richment analyses were conductedt el idate the b ological pr cesses in which target genes of DemiR were involved atvarious growth stages. GO annotations indicated that significantly enriched biologicalprocesses were related to cell growth and glycerol phosphate metabolism (Figure 3C).Remarkably, most of the previously mentioned target genes of 25 DEmiRs were involved inthese pathways. For instance, pathways associated with glycerol-3-phosphate metabolismwere enriched during both early and later developmental stages. In these pathways,two glycerol-3-phosphate dehydrogenases (GPDH) were regulated by ces-miR156b-3p,which showed low expression levels during early development but consistently highPlants 2024, 13, 3305 6 of 15expression from stages S3 to S5 (Figure 3B). Processes such as multicellular organismdevelopment, actin cytoskeleton organization, and base-excision repair were significantlyenriched at key early and late developmental stages (Figure 3C). These processes involveeight APETALA2/ethylene-responsive (AP2/ERF) transcription factors, two suppressorsof cAMP receptor (SCAR-LIKE) proteins, and two endonuclease III enzymes. Specifically,ces-miR172b-3p, ces-miR172d-3p, and ces-miRn2 target eight AP2, exhibiting high expres-sion only at stage S5 (Table S4). The two SCAR-LIKE proteins are regulated by ces-miRn108,which shows progressively increased expression levels during tuber growth. While noexpression differences are observed from stages S1 to S3, distinct differences emerge atstages S4 and S5 compared to earlier stages (Figure 3B). Endonuclease III is targeted byfive miRNAs: ces-miRn27, ces-miRn31, ces-miR166l-3p_1, miR166l-3p_2, and miR166k-3p.These miRNAs maintain high expression levels at stages S1 and S2, but their expressionis reduced to one-fourth of their previous abundance after S3 (Figure 3B and Table S4).The KEGG enrichment results were close to the GO analysis (Figure S2). Differing fromthese results, the pathway “starch and sucrose metabolism” was significantly enrichedaccording to the KEGG analysis. The highly expressed ces-miR168a-5p_4 in S4 and S5regulated three isoamylase 1 enzymes in this pathway (Figure 3B). Additionally, the folatebiosynthesis pathway was only not significantly enriched in the S1vsS2 and S4vsS5 compar-isons. This pathway involved three aminodeoxychorismate synthases which were targetedby ces-miRn25. This novel miRNA, like ces-miR168a-5p_4, showed a gradual increase inexpression during tuber growth (Figure 3B). In conclusion, the miRNA expression profilesat various growth stages corresponded to the curves of tuber growth and metabolites accu-mulation, and their target genes were mainly involved in cell growth, glycerophosphate,and sugar synthesis.2.4. miRNA–mRNA Regulatory Network Involved in Metabolite Sythesis and DevelopmentProgress in Tiger Nut TubersTo further understand the regulatory mechanisms of miRNAs on tiger nut tubergrowth and metabolite synthesis, a miRNA–mRNA regulatory network was constructed us-ing DEmiRs information and function annotations of target genes in Cytoscape v3.10.2 [23].A total of eight miRNAs were identified to be associated with the synthesis of fattyacids, amino acids, and carbohydrates (Figure 4A), and these miRNAs generally exhib-ited higher expression levels starting from S3 (Figure 3B and Table S4). In addition toces-miR156b-3p-Cyp0115470.1/Cyp0164100.1 (GPDH), the interaction pair related to fattyacid synthesis included ces-miRn92-Cyp0186910.1 [Lipin/Ned1/Smp2 (LNS2)] and ces-miRn117-Cyp0194050.1/Cyp0339760.1/Cyp0347520.1/Cyp0404730.1 [squamosa promoter-binding-like protein (SPL)]. Notably, the expression level of ces-miRn92 decreased onlyat S5 (Table S4). The amino acid metabolism involved two miRNA–mRNA interactionpairs: ces-miR396b-5p_2/ces-miRn24-Cyp0205270.1/Cyp0311610.1/Cyp0322150.1 [aminoacid permeases (AAP)] and ces-miRn25-Cyp0353620.1 [amino acid permeases (AAP)]. TwomiR168as are involved in carbohydrate metabolism. miR168a-5p_4 interacted with threeisoamylase 1 genes: Cyp0449600.1, Cyp0463730.1, and Cyp0523270.1. miR168a-3p regulatedthree alkaline neutral invertase CINV2-like proteins (CINV): Cyp0159510.1, Cyp0197230.1,and Cyp0205110.1. This miRNA only exhibited higher expression levels at S4 (Table S4).On the other hand, the miRNA–mRNA regulatory network related to growth anddevelopment processes included 22 miRNAs (Figure 4B); of these, 13 miRNAs exhib-ited lower expression levels after S3 (including S3) compared to earlier developmen-tal stages (S1 and S2) (Figure 3B and Table S4). The largest module consisted of fivemiR166 family members (ces-miR166a-3p, ces-miR166c-3p, ces-miR166k-3p, ces-miR166l-3p_1, and ces-miR166l-3p_2) and seven novel miRNAs (ces-miRn27, ces-miRn31, ces-miRn69, ces-miRn70, ces-miRn71, ces-miRn94, and ces-miRn98); their target genes includedfive HD-ZIP (Cyp0091390.1/Cyp0446160.1/Cyp0465030.1/Cyp0469270.1/Cyp0537450.1),two endonuclease III (Cyp0237730.1 and Cyp0414850.1), one pentatricopeptide repeat-containing protein/PPR (Cyp0381820.1), and one tetratricopeptide repeat protein/TPRPlants 2024, 13, 3305 7 of 15(Cyp0281570.1). The remaining miRNA–mRNA interactions only involved a single genefamily. Among these, ces-miR159f-5p_1, ces-miR159f-5p_2, and ces-miRn77 collectivelyacted on three PPR genes: Cyp0259970.1, Cyp0270890.1, and Cyp0284730.1. Meanwhile,ces-miR172b-3p and ces-miR172d-3p regulated six AP2 genes: Cyp0074770.1, Cyp0081000.1,Cyp0111500.1, Cyp0450180.1, Cyp0463160.1, and Cyp0522690.1. Notably, the expressionlevels of these two miRNAs were higher in S4 compared to the first three stages (S1, S2,and S3) (Table S4). In another interaction pair, ces-miRn2-Cyp0251950.1/Cyp0276990.1,ces-miRn2 also acted on two AP2 genes and its expression level in S5 was more thantwice as high as that in the other four stages (Table S4). Ces-miRn17, which was highlyexpressed in the later stages of development (Figure 3B), regulated five argonaute/AGOproteins (Cyp0381320.1/Cyp0418970.1/Cyp0488250.1/Cyp0506070.1/Cyp0529560.1). Ces-miRn15 participated in cell division by regulating five cell division cycle proteins/CDC(Cyp0104710.1/Cyp0140030.1/Cyp0240780.1/Cyp0263450.1/Cyp0295900.1), and its expres-sion level was only low in S1 (Table S4).Plants 2024, 13, x FOR PEER REVIEW 7 of 15 Figure 4. Hypothetical miRNA–mRNA interaction networks. (A) Metabolite synthesis of tubers. (B) Growth and development of tubers. On the other hand, the miRNA–mRNA regulatory network related to growth and development processes included 22 miRNAs (Figure 4B); of these, 13 miRNAs exhibited lower expression levels after S3 (including S3) compared to earlier developmental stages (S1 and S2) (Figure 3B and Table S4). The largest module consisted of five miR166 family members (ces-miR166a-3p, ces-miR166c-3p, ces-miR166k-3p, ces-miR166l-3p_1, and ces-miR166l-3p_2) and seven novel miRNAs (ces-miRn27, ces-miRn31, ces-miRn69, ces-miRn70, ces-miRn71, ces-miRn94, and ces-miRn98); their target genes included five HD-ZIP (Cyp0091390.1/Cyp0446160.1/Cyp0465030.1/Cyp0469270.1/Cyp0537450.1), two endonuclease III (Cyp0237730.1 and Cyp0414850.1), one pentatricopeptide repeat-containing protein/PPR (Cyp0381820.1), and one tetratricopeptide repeat protein/TPR (Cyp0281570.1). The remaining miRNA–mRNA interactions only involved a single gene family. Among these, ces-miR159f-5p_1, ces-miR159f-5p_2, and ces-miRn77 collectively acted on three PPR genes: Cyp0259970.1, Cyp0270890.1, and Cyp0284730.1. Meanwhile, ces-miR172b-3p and ces-miR172d-3p regulated six AP2 genes: Cyp0074770.1, Cyp0081000.1, Cyp0111500.1, Cyp0450180.1, Cyp0463160.1, and Cyp0522690.1. Notably, the expression levels of these two miRNAs were higher in S4 compared to the first three stages (S1, S2, and S3) (Table S4). In another interaction pair, ces-miRn2-Cyp0251950.1/Cyp0276990.1, ces-miRn2 also acted on two AP2 genes and its expression level in S5 was more than twice as high as that in the other four stages (Table S4). Ces-miRn17, which was highly expressed in the later stages of development (Figure 3B), regulated five argonaute/AGO proteins (Cyp0381320.1/Cyp0418970.1/Cyp0488250.1/Cyp0506070.1/Cyp0529560.1). Ces-miRn15 participated in cell division by regulating five cell division cycle proteins/CDC (Cyp0104710.1/Cyp0140030.1/Cyp0240780.1/Cyp0263450.1/Cyp0295900.1), and its expression level was only low in S1 (Table S4). 2.5. Quantitative Real-Time PCR (qRT-PCR) Validation of miRNA and Their Target Genes To assess the validity of small RNA sequencing, six miRNAs were randomly selected from two regulatory networks, and qRT-PCR was conducted on samples from five developmental stages of tiger nut tubers, with CeU6 serving as the housekeeping gene. The results showed that the TPM values and relative expression levels of these six miRNAs displayed consistent trends across all five developmental stages (Figure 5). These findings confirm the robustness of high-throughput sequencing. i re 4. Hypothetical miRNA–mRNA interactio networks. (A) Metabolite synthesis of tubers.(B) Growth and developmen of tubers.2.5. Quantitative Real-Time PCR (qRT-PCR) Validation of miRNA and Their Target GenesTo assess the validity of small RNA sequencing, six miRNAs were randomly selectedfrom two regulatory networks, and qRT-PCR was conducted on samples from five develop-mental stages of tiger nut tubers, with CeU6 serving as the housekeeping gene. The resultsshowed that the TPM values and relative expression levels of these six miRNAs displayedconsistent trends across all five developmental stages (Figure 5). These findings confirmthe robustness of high-throughput sequencing.Furthermore, to validate the miRNA–mRNA regulatory relationships, 10 target geneswere randomly selected from the predicted interaction network (Figure 4). The relativeexpression levels of these mRNAs varied according to tuber growth (Figure 6), with notabledifferences observed between critical periods of metabolite accumulation or development(one-way analysis of variance followed by Duncan’s multiple range tests, p < 0.05). Forinstance, Cyp0074770.1 showed significantly higher expression during S2 and S3 comparedto other stages, while being barely detectable in later developmental stages. The expressionanalysis revealed that the majority of miRNAs exhibited highly significant (Pearson’s corre-lation, p < 0.01) or significant (p < 0.05) negative correlations with their target genes (Figure 6and Table S7). However, three mRNA–miRNA pairs showed no significant correlations:Cyp0469270.1 and Cyp0091390.1 with ces-miRn69 and ces-miRn71, and Cyp0115470.1 withces-miR156b-3p. Nevertheless, both Cyp0469270.1 and Cyp0091390.1 were identified aspotential targets of multiple miRNAs (Figure 4), and all these relationships demonstratedhighly significant negative correlations (Table S7).Plants 2024, 13, 3305 8 of 15Plants 2024, 13, x FOR PEER REVIEW 8 of 15 Figure 5. qRT-PCR validation of the six DEmiRs at five development stages (bar plot shows miRNA expression from small RNA sequencing and line plot shows relative expression from qRT-PCR). Furthermore, to validate the miRNA–mRNA regulatory relationships, 10 target genes were randomly selected from the predicted interaction network (Figure 4). The relative expression levels of these mRNAs varied according to tuber growth (Figure 6), with notable differences observed between critical periods of metabolite accumulation or development (one-way analysis of variance followed by Duncan’s multiple range tests, p < 0.05). For instance, Cyp0074770.1 showed significantly higher expression during S2 and S3 compared to other stages, while being barely detectable in later developmental stages. The expression analysis revealed that the majority of miRNAs exhibited highly significant (Pearson’s correlation, p < 0.01) or significant (p < 0.05) negative correlations with their target genes (Figure 6 and Table S7). However, three mRNA–miRNA pairs showed no significant correlations: Cyp0469270.1 and Cyp0091390.1 with ces-miRn69 and ces-miRn71, and Cyp0115470.1 with ces-miR156b-3p. Nevertheless, both Cyp0469270.1 and Cyp0091390.1 were identified as potential targets of multiple miRNAs (Figure 4), and all these relationships demonstrated highly significant negative correlations (Table S7). Figure 5. qRT-PCR validation of the six DEmiRs at five development stages (bar plot shows miRNAexpression from small RNA sequencing and line plot shows relative xpression from qRT-PCR).Plants 2024, 13, x FOR PEER REVIEW 8 of 15 Figure 5. qRT-PCR validation of the six DEmiRs at five development stages (bar plot shows miRNA expression from small RNA sequencing and line plot shows relative expression from qRT-PCR). Furthermore, to validate the miRNA–mRNA regulatory relationships, 10 target genes were randomly selected from the predicted interaction network (Figure 4). The relative expression levels of these mRNAs varied according to tuber growth (Figure 6), with notable differences observed between critical periods of metabolite accumulation or development (one-way analysis of variance followed by Duncan’s multiple range tests, p < 0.05). For instance, Cyp0074770.1 showed significantly higher expression during S2 and S3 compared to other stages, while being barely detectable in later developmental stages. The expression analysis revealed that the majority of miRNAs exhibited highly significant (Pearson’s correlation, p < 0.01) or significant (p < 0.05) negative correlations with their target genes (Figure 6 and Table S7). However, three mRNA–miRNA pairs showed no significant correlations: Cyp0469270.1 and Cyp0091390.1 with ces-miRn69 and ces-miRn71, and Cyp0115470.1 with ces-miR156b-3p. Nevertheless, both Cyp0469270.1 and Cyp0091390.1 were identified as potential targets of multipl miRNAs (Figure 4), and all th s relationships demonstrated highly significant negative correlati s (Table S7). Figure 6. Expression analysis between target genes and their corresponding miRNA. The line plotsrepresent the relative expression levels of mRNAs, with data shown as means ± standard deviationfrom three independent experiments. Different lowercase letters indicate significant differences(p < 0.05, Duncan’s multiple range tests). The bar charts show miRNA expression levels, and thecorrelations between mRNAs and miRNAs were analyzed using Pearson’s correlation coefficient(p < 0.05 indicates significant correlation and p < 0.01 indicates highly significant correlation). Fortarget genes regulated by multiple miRNAs, only the correlation with the lowest p value is shown inthe figure.3. DiscussionTiger nut is a plant known for its ability to accumulate substantial amounts of oil inits tuber [24], making it essential for understanding the patterns and mechanisms of oilproduction. This study utilized GC–MS to determine the accumulation and changes in theoil content and fatty acid composition of tubers at five developmental stages. ConsistentPlants 2024, 13, 3305 9 of 15with previous findings, the tubers are rich in various unsaturated fatty acids, with oleicacid constituting more than 50% of the total oil content, and thus being the predominantcomponent [6,9]. Moreover, both dry matter weight and total oil content increase rapidlyduring the mid-developmental stages and remain elevated for over 40 days (Figure 1B) [6,8,9].Among the 14 detected fatty acids, this trend is particularly evident for oleic acid [6].Notably, in the later developmental stages, the levels of major fatty acids such as oleicacid, linoleic acid, and palmitic acid decrease, while those of minor fatty acids increase(Figures 1C,D and S1). This may reflect shifts in the synthesis and consumption of fattyacid components. Overall, our results suggest that 65 DAS to 105 DAS (from 35 to 75 daysafter tuber formation) are critical for metabolite accumulation in tiger nut tubers.miRNAs have been demonstrated to influence the oil content and composition ofplant seeds and are crucial in potato tuber development [12,14,15,25]. In this study, weidentified 183 miRNAs using high-throughput sequencing and mapped their expressionprofiles (Tables S3 and S4). Among these, miR166, miR168, and miR156 were identifiedin 11 miRNA families within the tubers (Figure 2A). These miRNAs are also found in oil-rich plants such as rapeseed [16], sesame [26], and peony [17], indicating their conservednature. Nonetheless, a significant number of unidentified miRNA types were observedin tiger nut (Table S3). Some miRNAs known to be involved in plant fatty acid andoil production, such as miR169 and miR390, were not detected, which may be due tointerspecies differences [16,26]. These results underscore the need for small RNA studiesacross various tissues and cultivars to enhance our understanding of miRNA types andfunctions in diverse plants.The analysis of differential miRNA expression and the functions of their target genesat various developmental stages in plants provide insights into how miRNAs mediatebiological functions [17,18]. For instance, the rapid oil accumulation stage in sesame seedscoincides with the most significant changes in miRNA expression, with these miRNAs’target genes being involved in numerous biological processes related to lipid synthesis [26].Similarly, in tiger nut tubers, the expression changes of many miRNAs across five develop-mental stages closely align with the rate of metabolite accumulation (Figure 3B). Duringthe stages of the most rapid synthesis of dry matter, total oil, and fatty acids, there is anincreased number of DEmiRs (Figure 3A). These miRNAs directly target genes involved inthe glycerol-3-phosphate catabolic process and cell growth (Figure 3C and Table S5). Theces-miR156-3p directly targets two GPDH genes (Figure 4A), with its expression signifi-cantly increasing at 85 DAS and remaining at high levels thereafter (Figure 3B). As a keyenzyme in the Kennedy pathway, GPDH catalyzes the synthesis of glycerol-3-phosphate.In soybeans, GmGPDH overexpression can directly increase diacylglycerol content andpromote oleic acid synthesis, leading to simultaneous increases in both soybean oil andoleic acid content [27]. Unexpectedly, although the GPDH gene Cyp0115470.1 shows higherrelative expression levels during S3 and S4 compared to other stages, it exhibits only a weakpositive correlation with ces-miR156-3p expression (Figure 6). These findings suggest thatthe interaction between GPDH, a component of the oil biosynthesis pathway in tiger nuttubers, and ces-miR156-3p might be more complex than anticipated. Ces-miRn117, whichshares a similar expression pattern with ces-miR156-3p, targets five SPL genes (Figure 4A).In rapeseed, miR156 interacts with SPL2/10 and SPL10/11 to regulate oil production andseed maturation, respectively [28,29]. StSPL9 is involved in potato tuber formation throughtis regulation of the expression levels of miR172. Cyp0347520.1, identified as an SPL gene,exhibited high expression levels during the S2 stage and was progressively suppressed asthe tuber matured (Figure 6). This expression pattern suggests that it might play a rolein regulating oil accumulation during the rapid oil biosynthesis phase (Figure 6). Ces-miR172b-3p and ces-miR172d-3p are key modules in tuber growth, targeting six AP2 genes(Figure 4B). Two reports indicate that AP2 genes not only control the size and number ofArabidopsis seeds but also influence dry mass and oil accumulation [30,31]. Cyp0074770.1,an AP2 gene, as a target gene of two Ces-miR172 family members, was found to be tran-scriptionally repressed (Figure 6). This gene exhibits higher expression levels during the S2Plants 2024, 13, 3305 10 of 15to S3 stages compared to other developmental stages, with minimal expression detected inthe later developmental phases. Furthermore, the miR172 family has been identified as akey factor mediating fat formation in oilseed crop seeds [16,17,26]. Thus, ces-miRn117-SPLand ces-miR172-AP2 may interact to jointly regulate the oil synthesis and development oftuber.In addition, ces-miR168a-5p_4, which regulates four genes encoding isoamylase,involved in starch accumulation in seeds or tubers [32,33], exhibited a rapid increasein expression to twice its previous level at 85 DAS and further increased as the tubersmatured (Figure 3B and Table S4). These results suggest that miR168a-5p_4 is involvedin the tuber starch synthesis pathway. Notably, several studies have demonstrated thatDGAT [10,11] and SAD [7] play a crucial role in oil accumulation in tiger nut tubers.However, the miRNAs identified in this study did not act on DGAT or SAD, highlightingthe complexity of oil production mechanisms in tiger nut tubers. Thus, future researchshould conduct transcriptomic and other omics studies across different varieties of tigernut, particularly among those with significant differences in oil content, to elucidate theregulatory mechanisms of miRNAs in tuber oil accumulation.miR166 family members are highly conserved and commonly found in oilseed crops,representing a significant group of miRNAs [17,19]. The miR166-HD-ZIP module is akey regulator of plant organ growth [13,15,21]. In this study, miR166, one of the mostabundant conserved miRNA families, showed that the expression levels of five membersgradually decreased during development or were repressed after 85 DAS, maintaining theirexpression at less than 1.3-fold of the previous level or becoming undetectable (Figure 3Band Table S4). Two HD-ZIP family target genes, Cyp0091390.1 and Cyp0469270.1, whichwere negatively regulated by miR166, exhibited progressively increased expression levelsduring tuber maturation. The relative expression levels of these genes at the S5 stage weremore than three-fold higher compared to those at the S1 stage (Figure 6 and Table S7). Re-markably, ces-miRn17 not only regulates five AGO proteins but also displays an expressionpattern that is opposite to that of the five miR166 members (Figures 3B and 4B and Table S4).Moreover, one of the AGO genes, Cyp0506070.1, exhibited high expression levels exclusivelyduring the S1 and S2 stages (Figure 6). Previous studies have demonstrated that AtAGO1regulates vegetative organ development by activating miR166, which in turn suppressesHD-ZIP gene expression [15]. Hence, the ces-miRn17-AGO-ces-miR166-HD-ZIP modulemay form a similar regulatory network to influence the growth and development of tigernut tubers. Meanwhile, these miR166 members, together with ces-miRn27, ces-miRn31,and others, form an extensive regulatory network targeting HD-ZIP, endonuclease III, PPR,and TPR genes (Figure 4B). The loss of the AtTPR in Arabidopsis can cause abnormalcell proliferation [34]. In addition, the CDC48 gene directly regulates the cell divisionprocess [35,36]. Ces-miR15 targets five CDC48 genes, with its expression level increas-ing twofold after 65 DAS and remaining elevated until tuber maturation (Figure 4B andTable S4). In contrast, the expression of the CDC48 gene (Cyp0263450.1) showed an inversepattern to Ces-miR15, with its expression level sharply declining after S1 and becomingbarely detectable in subsequent stages (Figure 6). Therefore, it can be inferred that celldivision in tiger nut tubers primarily occurs during the first 35 days of tuber formation,with ces-miR15 playing a crucial role. In conclusion, the results demonstrate that miRNAsare closely associated with oil accumulation, fatty acid composition, and growth duringthe developmental stages of tiger nut tubers. However, further validation is required toelucidate how miRNA–mRNA interactions influence these traits.4. Materials and Methods4.1. Plant MaterialsThe tiger nut cultivar, JINONGSHA6, was planted in the experimental field of JilinAgricultural University (Jilin Province, China). In 2022, tiger nut tubers began to form30 DAS and reached maturity at 125 DAS. Consequently, five developmental stages wereselected for this study: 45 DAS (S1), 65 DAS (S2), 85 DAS (S3), 105 DAS (S4), and 125 DASPlants 2024, 13, 3305 11 of 15(S5). Uniform-sized tubers were obtained from three plants exhibiting similar growth,washed with distilled water to remove soil, and dried with absorbent paper. They werethen immediately frozen in liquid nitrogen and subsequently stored in a −80 ◦C freezer forfurther analysis4.2. Analysis of Oil Content and Fatty Acid CompositionA total of 200 mg of tiger nut tuber powder was mixed with 5 mL of 5% hydrochloricacid methanol solution, 10 mL of dichloromethane, and 15 mL of methanol. The mixturewas heated in an 80 ◦C water bath for 1 h and then cooled to 25 ◦C. The solution wasadjusted to a final volume of 20 mL with methanol to obtain the test solution. The samplesolution was prepared by mixing 1 mL of the test solution with an internal standard solution(methyl nonadecanoate) and 5% aqueous ammonia. The oil content and fatty acid profile ofthe samples were determined by GC–MS using a Trace 1310 ISQ (Thermo Fisher Scientific,Waltham, MA, USA). Separation was performed on a TG-5MS capillary column (30 m ×0.25 mm × 0.25 µm) with helium as the carrier gas at a flow rate of 1.2 mL/min. A splitinjection mode with a split ratio of 1:50 was used, and the injector temperature was set to290 ◦C. The oven temperature was programmed to start at 80 ◦C for 1 min, ramp to 200 ◦Cat a rate of 10 ◦C/min, hold at 250 ◦C for 5 min, then increase to 270 ◦C at a rate of 2 ◦C/min,and hold at 270 ◦C for 3 min. The temperatures of the transfer line and ion source wereset to 280 ◦C. The mass spectrometry scanning range was 30–400 amu. All experimentswere performed in triplicate. A statistical analysis of tuber traits was conducted using the Rv4.3.1 (https://www.r-project.org/ (accessed on 6 October 2024)). For comparisons amongmultiple groups, one-way analysis of variance (ANOVA) and Duncan’s multiple range testwere conducted using the R package agricolae [37].4.3. RNA Isolation and Small RNA SequencingTotal RNA from 15 samples at five tuber developmental stages was extracted usingthe RN40-EASYspin reagent (Aidlab, Beijing, China) according to the manufacturer’sinstructions. The quality and quantity of the total RNA from each sample were analyzedusing a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).Small RNA was separated from high-quality RNA and used to construct small RNAlibraries with the NEBNext Ultra Small RNA Library Prep Kit for Illumina (NEB, Ipswich,MA, USA). Sequencing was performed on the Illumina NovaSeq 6000 platform in single-end 50 bp (Biomarker Technology, Beijing, China).Raw reads were processed by removing the adaptor sequence using Trimmomaticv0.39 with the “SE” module [38]. Low-quality reads, poly-A tails, shorter than 18bp, andlonger than 34 bp were filtered using Cutadapt v4.1 [39] with the parameter setting “-aA{20} -m 18 -M 34”. Clean reads from all samples were combined to identify known andnovel miRNAs. Combined reads were converted into unique tags using the FASTX-Toolkit0.0.14 (https://www.encodeproject.org/software/fastx_toolkit/ (accessed on 6 October2024)). Ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), smallnucleolar RNA (snoRNA), and repeat sequences were filtered by aligning unique tagsto the Rfam database [40] and Repbase database [41]. Known miRNAs in tiger nut wereidentified based on known plant miRNAs in the miRbase database v22 [42]. To accuratelyidentify novel miRNAs, the plant-specific miRNA analysis tool miRDeep-P2 [43] was usedto identify novel miRNAs with the minor modification [44]. The chromosome-level genomeassembly sequence of tiger nut (unpublished) was available to detect novel miRNA. Atotal of 68,504 coding sequences from the gene annotation were applied to predict miRNAtargets using psRNATarget [45]. The eggNOG-mapper v2 [46] and plantTFDB v5.0 [47]were used for the functional annotation of genes and identification of transcription factorsbased on gene protein sequences.Plants 2024, 13, 3305 12 of 154.4. Differential Expression Analysis of miRNAsAfter the identification of known and novel miRNAs, the expression of miRNAs ineach sample was quantified by aligning clean reads to the final set of miRNA sequencesusing Bowtie v1.3.1 with default parameters [48]. The read counts for each miRNA ineach sample were generated using SAMtools [49]. The expression level of miRNA wasnormalized to TPM. The differential expression of miRNAs was analyzed based on readcounts using DESeq2 across different developmental stages [50,51]. Differentially expressedmiRNAs were identified using an absolute value of FC ≥ 1.5 and a p-value < 0.05. GO andKEGG enrichment analyses were performed for the target genes of differentially expressedmiRNAs between different stages using the clusterProfiler 4.0 [52]. All figures were plottedusing ggplot2.4.5. Expression Analysis of miRNA and Their Target Genes Using qRT-PCRThe qRT-PCR analysis of miRNAs and their target genes was performed usinga QuantStudioTM 5 Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA,USA). For miRNAs, residual total RNA from small RNA library construction was reverse-transcribed to cDNA using TransScript miRNA First-Strand cDNA Synthesis SuperMix(AT351, TransGen, Beijing, China) according to the manufacturer’s protocol. qRT-PCRwas conducted using TransStart Top Green qPCR SuperMix (AQ132-11, TransGen, Beijing,China), with CeU6 as the internal reference control. The qRT-PCR cycling conditions wereas follows: initial denaturation at 95 ◦C for 2 min, followed by 45 cycles of 95 ◦C for 15 s and60 ◦C for 30 s. For mRNA analysis, total RNA was reverse-transcribed using StarScript IIIRT Kit (A232, GeneStar, Beijing, China), with CeUCE2 as the reference gene [53]. Subsequentsteps were identical to those for miRNA. All samples were subjected to three technicalreplicates. Relative expression levels were calculated using the relative quantificationmethod (2−∆∆Ct) [54]. The primers used for qRT-PCR are listed in Supplementary Table S6.Additionally, Pearson correlation coefficients between the relative expression of mRNAsand the TPM of their corresponding miRNA were calculated using the R package corrplot0.92.5. ConclusionsThis study reveals that metabolite accumulation, predominantly oleic acid, in C.esculentus tubers primarily occurs from 65 DAS to 105 DAS. Small RNA sequencing acrossfive developmental stages identified 183 miRNAs, comprising 31 known miRNAs from11 families and 152 novel miRNAs. Significant differences in miRNA expression wereobserved during key periods of metabolite accumulation. The target genes of these miRNAsare significantly enriched in glycerate 3-phosphate synthesis and cell growth processes,and 13 related regulatory networks were identified. For instance, ces-miR156b-3p-GDPHand ces-miRn117-SPL are associated with total oil content and oleic acid synthesis, whileces-miRn17-AGO-ces-miR166-HD-ZIP and ces-miR15-CDC48 regulate tuber growth anddevelopment. miRNAs primarily participate in biological processes by suppressing theexpression of their target genes. These findings provide insights into the key miRNAsrelated to metabolite accumulation and their regulatory networks with target genes incrops where tubers are the main economic organ, enhancing our understanding of thedevelopmental and oil production mechanisms in tiger nut tubers.Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants13233305/s1, Figure S1: Content of various fatty acids oftiger nut tubers at five development stages; Figure S2: KEGG enrichment analysis of differentiallyexpressed miRNAs, Table S1: Fatty acid composition of tiger nut tubers at different developmentalstages; Table S2: Summary of small RNA sequencing quality; Table S3: Basic information of 183 iden-tified miRNAs; Table S4: Expression levels of miRNAs at different developmental stages; Table S5:Target genes of miRNAs and their annotation information; Table S6: Primers for qRT-PCR analysis ofmiRNAs and their target genes; Table S7: Correlation between miRNAs and their target genes.Plants 2024, 13, 3305 13 of 15Author Contributions: Conceptualization, L.W. and S.Y.; validation, Y.G. and S.C.; formal analysis,L.W.; investigation, L.W., Y.G., X.L., W.C. and J.Z. (Jian Zhang); resources, Y.G., X.L., X.N. andW.C.; data curation, Y.G. and L.C.; writing—original draft preparation, L.W.; writing—review andediting, L.W., L.C. and S.Y.; visualization, L.C. and S.C.; supervision, S.Y. and J.Z. (Jun Zhang); projectadministration, L.W. and X.N.; funding acquisition, J.Z. (Jun Zhang). All authors have read andagreed to the published version of the manuscript.Funding: This research was funded by innovation of high-yield and high-quality tiger nut germplasmand breeding of new cultivars from Jilin Provincial Development and Reform Commission (2023C035-5).Data Availability Statement: All small RNA sequencing data have been deposited at the NGDCwebsite (https://ngdc.cncb.ac.cn/), and the bioproject number is PRJCA031104.Conflicts of Interest: The authors declare no conflicts of interest.References1. Sánchez-Zapata, E.; Fernández-López, J.; Angel Pérez-Alvarez, J. Tiger Nut (Cyperus esculentus) Commercialization: HealthAspects, Composition, Properties, and Food Applications. Compr. Rev. Food. Sci. Food Saf. 2012, 11, 366–377. [CrossRef]2. Linssen, J.P.H.; Cozijnsen, J.L.; Pilnik, W. Chufa (Cyperus esculentus): A new source of dietary fibre. J. Sci. Food Agric. 1989, 49,291–296. [CrossRef]3. Nina, G.C.; Goncharov, A.; Batishcheva, N.; Vlasov, S.; Okuskhanova, E.; CÃsarovÃ, M.; Hleba, L.; Ogori, A.F.; Shariati, M.A.Proximate, mineral and functional properties of tiger nut flour extracted from different tiger nuts cultivars. J. Microbiol. Biotechnol.Food Sci. 2019, 9, 653–656. [CrossRef]4. Zhang, S.; Li, P.; Wei, Z.; Cheng, Y.; Liu, J.; Yang, Y.; Wang, Y.; Mu, Z. Cyperus (Cyperus esculentus L.): A Review of Its Compositions,Medical Efficacy, Antibacterial Activity and Allelopathic Potentials. Plants 2022, 11, 1127. [CrossRef]5. Makareviciene, V.; Gumbyte, M.; Yunik, A.; Kalenska, S.; Kalenskii, V.; Rachmetov, D.; Sendzikiene, E. Opportunities for the useof chufa sedge in biodiesel production. Ind. Crop. Prod. 2013, 50, 633–637. [CrossRef]6. Wang, L.; Jing, M.; Ahmad, N.; Wang, Y.; Wang, Y.; Li, J.; Li, X.; Liu, W.; Wang, N.; Wang, F.; et al. Tracing Key MolecularRegulators of Lipid Biosynthesis in Tuber Development of Cyperus esculentus Using Transcriptomics and Lipidomics Profiling.Genes 2021, 12, 1492. [CrossRef]7. Li, T.; Sun, Y.; Chen, Y.; Gao, Y.; Gao, H.; Liu, B.; Xue, J.; Li, R.; Jia, X. Characterisation of two novel genes encoding ∆9 fatty aciddesaturases (CeSADs) for oleic acid accumulation in the oil-rich tuber of Cyperus esculentus. Plant Sci. 2022, 319, 111243. [CrossRef]8. Ji, H.; Liu, D.; Yang, Z. High oil accumulation in tuber of yellow nutsedge compared to purple nutsedge is associated withmore abundant expression of genes involved in fatty acid synthesis and triacylglycerol storage. Biotechnol. Biofuels 2021, 14, 54.[CrossRef]9. Turesson, H.; Marttila, S.; Gustavsson, K.; Hofvander, P.; Olsson, M.E.; Bülow, L.; Stymne, S.; Carlsson, A.S. Characterization ofoil and starch accumulation in tubers of Cyperus esculentus var. sativus (Cyperaceae): A novel model system to study oil reserves innonseed tissues. Am. J. Bot. 2010, 97, 1884–1893. [CrossRef]10. Gao, Y.; Sun, Y.; Gao, H.; Chen, Y.; Wang, X.; Xue, J.; Jia, X.; Li, R. Ectopic overexpression of a type-II DGAT (CeDGAT2-2) derivedfrom oil-rich tuber of Cyperus esculentus enhances accumulation of oil and oleic acid in tobacco leaves. Biotechnol. Biofuels 2021, 14,76. [CrossRef]11. Liu, D.; Ji, H.; Yang, Z. Functional Characterization of Three Novel Genes Encoding Diacylglycerol Acyltransferase (DGAT) fromOil-Rich Tubers of Cyperus esculentus. Plant Cell Physiol. 2020, 61, 118–129. [CrossRef] [PubMed]12. Martin, A.; Adam, H.; Díaz-Mendoza, M.; Zurczak, M.; González-Schain, N.D.; Suárez-López, P. Graft-transmissible induction ofpotato tuberization by the microRNA miR172. Development 2009, 136, 2873–2881. [CrossRef] [PubMed]13. Kitazumi, A.; Kawahara, Y.; Onda, T.S.; De Koeyer, D.; de Los Reyes, B.G. Implications of miR166 and miR159 induction to thebasal response mechanisms of an andigena potato (Solanum tuberosum subsp. andigena) to salinity stress, predicted from networkmodels in Arabidopsis. Genome 2015, 58, 13–24. [CrossRef] [PubMed]14. Bhogale, S.; Mahajan, A.S.; Natarajan, B.; Rajabhoj, M.; Thulasiram, H.V.; Banerjee, A.K. MicroRNA156: A Potential Graft-Transmissible MicroRNA That Modulates Plant Architecture and Tuberization in Solanum tuberosum ssp. andigena. Plant Physiol.2014, 164, 1011–1027. [CrossRef] [PubMed]15. Tang, X.; Bian, S.; Tang, M.; Lu, Q.; Li, S.; Liu, X.; Tian, G.; Nguyen, V.; Tsang, E.W.T.; Wang, A.; et al. MicroRNA–MediatedRepression of the Seed Maturation Program During Vegetative Development in Arabidopsis. PLoS Genet. 2012, 8, e1003091.[CrossRef]16. Wang, J.; Jian, H.; Wang, T.; Wei, L.; Li, J.; Li, C.; Liu, L. Identification of microRNAs Actively Involved in Fatty Acid Biosynthesisin Developing Brassica napus Seeds Using High-Throughput Sequencing. Front. Plant Sci. 2016, 7, 1570. [CrossRef]17. Yin, D.; Li, S.; Shu, Q.; Gu, Z.; Wu, Q.; Feng, C.; Xu, W.; Wang, L. Identification of microRNAs and long non-coding RNAsinvolved in fatty acid biosynthesis in tree peony seeds. Gene 2018, 666, 72–82. [CrossRef]18. Zhang, T.; Li, Z.; Song, X.; Han, L.; Wang, L.; Zhang, J.; Long, Y.; Pei, X. Identification and Characterization of microRNAs in theDeveloping Seed of Linseed Flax (Linum usitatissimum L.). Int. J. Mol. Sci. 2020, 21, 2708. [CrossRef]Plants 2024, 13, 3305 14 of 1519. Zhao, X.; Yang, G.; Liu, X.; Yu, Z.; Peng, S. Integrated Analysis of Seed microRNA and mRNA Transcriptome Reveals ImportantFunctional Genes and microRNA-Targets in the Process of Walnut (Juglans regia) Seed Oil Accumulation. Int. J. Mol. Sci. 2020, 21,9093. [CrossRef]20. Li, J.; Ding, J.; Yu, X.; Li, H.; Ruan, C. Identification and expression analysis of critical microRNA-transcription factor regulatorymodules related to seed development and oil accumulation in developing Hippophae rhamnoides seeds. Ind. Crop. Prod. 2019, 137,33–42. [CrossRef]21. Zhang, L.; Zhang, C.; Yang, B.; Chen, S.; Yang, Z.; Kang, L.; Xu, J.; Sun, T.; Kong, X.; Qu, C.; et al. Comprehensive high-throughputsequencing, evolutionary and functional analyses reveal the conservation and diversification of miR166s in regulating pleiotropictraits between rapeseed and Arabidopsis. Ind. Crop. Prod. 2024, 218, 118817. [CrossRef]22. Ozseyhan, M.E.; Li, P.; Na, G.; Li, Z.; Wang, C.; Lu, C. Improved fatty acid profiles in seeds of Camelina sativa by artificialmicroRNA mediated FATB gene suppression. Biochem. Biophys. Res. Commun. 2018, 503, 621–624. [CrossRef]23. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: Asoftware environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [CrossRef]24. Ezeh, O.; Gordon, M.H.; Niranjan, K. Tiger nut oil (Cyperus esculentus L.): A review of its composition and physico-chemicalproperties. Eur. J. Lipid Sci. Technol. 2014, 116, 783–794. [CrossRef]25. Belide, S.; Petrie, J.R.; Shrestha, P.; Singh, S.P. Modification of Seed Oil Composition in Arabidopsis by Artificial microRNA-Mediated Gene Silencing. Front. Plant Sci. 2012, 3, 168. [CrossRef] [PubMed]26. Zhang, Y.; Zhang, Y.; Thakur, K.; Zhang, F.; Hu, F.; Zhang, J.; Wei, P.; Wei, Z. Integration of miRNAs, Degradome, andTranscriptome Omics Uncovers a Complex Regulatory Network and Provides Insights Into Lipid and Fatty Acid SynthesisDuring Sesame Seed Development. Front. Plant Sci. 2021, 12, 709197. [CrossRef]27. Zhao, Y.; Cao, P.; Cui, Y.; Liu, D.; Li, J.; Zhao, Y.; Yang, S.; Zhang, B.; Zhou, R.; Sun, M.; et al. Enhanced production of seed oil withimproved fatty acid composition by overexpressing NAD+-dependent glycerol-3-phosphate dehydrogenase in soybean. J. Integr.Plant Biol. 2021, 63, 1036–1053. [CrossRef]28. Huang, D.; Koh, C.; Feurtado, J.A.; Tsang, E.W.; Cutler, A.J. MicroRNAs and their putative targets in Brassica napus seedmaturation. Bmc Genom. 2013, 14, 140. [CrossRef] [PubMed]29. Zhao, Y.; Wang, M.; Fu, S.; Yang, W.; Qi, C.; Wang, X. Small RNA Profiling in Two Brassica napus Cultivars Identifies MicroRNAswith Oil Production- and Development-Correlated Expression and New Small RNA Classes. Plant Physiol. 2012, 158, 813–823.[CrossRef]30. Ohto, M.; Fischer, R.L.; Goldberg, R.B.; Nakamura, K.; Harada, J.J. Control of seed mass by APETALA2. Proc. Natl. Acad. Sci. USA2005, 102, 3123–3128. [CrossRef]31. Jofuku, K.D.; Omidyar, P.K.; Gee, Z.; Okamuro, J.K. Control of seed mass and seed yield by the floral homeotic gene APETALA2.Proc. Natl. Acad. Sci. USA 2005, 102, 3117–3122. [CrossRef] [PubMed]32. Zeeman, S.C.; Umemoto, T.; Lue, W.; Au-Yeung, P.; Martin, C.; Smith, A.M.; Chen, J. A Mutant of Arabidopsis Lacking aChloroplastic Isoamylase Accumulates Both Starch and Phytoglycogen. Plant Cell 1998, 10, 1699–1711. [CrossRef] [PubMed]33. Bustos, R.; Fahy, B.; Hylton, C.M.; Seale, R.; Nebane, N.M.; Edwards, A.; Martin, C.; Smith, A.M. Starch granule initiation iscontrolled by a heteromultimeric isoamylase in potato tubers. Proc. Natl. Acad. Sci. USA 2004, 101, 2215–2220. [CrossRef][PubMed]34. Zhang, M.; Wang, C.; Lin, Q.; Liu, A.; Wang, T.; Feng, X.; Liu, J.; Han, H.; Ma, Y.; Bonea, D.; et al. A tetratricopeptide repeatdomain-containing protein SSR1 located in mitochondria is involved in root development and auxin polar transport in Arabidopsis.Plant J. 2015, 83, 582–599. [CrossRef] [PubMed]35. Park, S.; Rancour, D.M.; Bednarek, S.Y. In Planta Analysis of the Cell Cycle-Dependent Localization of AtCDC48A and Its CriticalRoles in Cell Division, Expansion, and Differentiation. Plant Physiol. 2008, 148, 246–258. [CrossRef]36. Shi, L.; Zhang, X.; Shi, Y.; Xu, X.; He, Y.; Shao, G.; Huang, Q.; Wu, J. OsCDC48/48E complex is required for plant survival in rice(Oryza sativa L.). Plant Mol. Biol. 2019, 100, 163–179. [CrossRef]37. De Mendiburu, F.; Agricolae, S.R. Agricolae—Ten years of an open source statistical tool for experiments in breeding, agricultureand biology. PeerJ 2015, 3, e1404. [CrossRef]38. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120.[CrossRef]39. Marcel, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. Embnet. J. 2011, 17, 10–12. [CrossRef]40. Kalvari, I.; Nawrocki, E.P.; Ontiveros-Palacios, N.; Argasinska, J.; Lamkiewicz, K.; Marz, M.; Griffiths-Jones, S.; Toffano-Nioche,C.; Gautheret, D.; Weinberg, Z.; et al. Rfam 14: Expanded coverage of metagenomic, viral and microRNA families. Nucleic AcidsRes. 2020, 49, D192–D200. [CrossRef]41. Bao, W.; Kojima, K.K.; Kohany, O. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob. DNA 2015, 6,11. [CrossRef] [PubMed]42. Kozomara, A.; Birgaoanu, M.; Griffiths-Jones, S. miRBase: From microRNA sequences to function. Nucleic Acids Res. 2018, 47,D155–D162. [CrossRef]43. Kuang, Z.; Wang, Y.; Li, L.; Yang, X. miRDeep-P2: Accurate and fast analysis of the microRNA transcriptome in plants.Bioinformatics 2019, 35, 2521–2522. [CrossRef]Plants 2024, 13, 3305 15 of 1544. Garg, V.; Dudchenko, O.; Wang, J.; Khan, A.W.; Gupta, S.; Kaur, P.; Han, K.; Saxena, R.K.; Kale, S.M.; Pham, M.; et al. Chromosome-length genome assemblies of six legume species provide insights into genome organization, evolution, and agronomic traits forcrop improvement. J. Adv. Res. 2022, 42, 315–329. [CrossRef] [PubMed]45. Dai, X.; Zhuang, Z.; Zhao, P.X. psRNATarget: A plant small RNA target analysis server (2017 release). Nucleic Acids Res. 2018, 46,W49–W54. [CrossRef] [PubMed]46. Cantalapiedra, C.P.; Hernández-Plaza, A.; Letunic, I.; Bork, P.; Huerta-Cepas, J. eggNOG-mapper v2: Functional Annotation,Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol. Biol. Evol. 2021, 38, 5825–5829. [CrossRef]47. Feng, T.; Dechang, Y.; Yuqi, M.; Jinpu, J.; Ge, G. PlantRegMap: Charting functional regulatory maps in plants. Nucleic Acids Res.2019, 48, D1104–D1113. [CrossRef]48. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to thehuman genome. Genome Biol. 2009, 10, R25. [CrossRef]49. Danecek, P.; Bonfield, J.K.; Liddle, J.; Marshall, J.; Ohan, V.; Pollard, M.O.; Whitwham, A.; Keane, T.; Mccarthy, S.A.; Davies, R.M.Twelve years of SAMtools and BCFtools. Gigascience 2021, 10, giab008. [CrossRef]50. Love, M.; Anders, S.; Huber, W. Differential analysis of count data–the DESeq2 package. Genome Biol. 2014, 15, 10–1186.51. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. GenomeBiol. 2014, 15, 550. [CrossRef] [PubMed]52. Wu, T.; Hu, E.; Xu, S.; Chen, M.; Guo, P.; Dai, Z.; Feng, T.; Zhou, L.; Tang, W.; Zhan, L.; et al. clusterProfiler 4.0: A universalenrichment tool for interpreting omics data. Innovation 2021, 2, 100141. [CrossRef] [PubMed]53. Bai, X.; Chen, T.; Wu, Y.; Tang, M.; Xu, Z. Selection and Validation of Reference Genes for qRT-PCR Analysis in the Oil-Rich TuberCrop Tiger Nut (Cyperus esculentus) Based on Transcriptome Data. Int. J. Mol. Sci. 2021, 22, 2569. [CrossRef] [PubMed]54. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−∆∆CTMethod. Methods 2001, 25, 402–408. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Researcher
Cyperus esculentus
miRNAs
GC-MS
oil synthesis
tuber development
Identification of miRNAs Involved in Lipid Metabolism and Tuber Development in Cyperus esculentus L.
Text