Full-Length Transcriptome and the Identiﬁcation of lncRNAs Involved in Salicylic Acid-Induced Flowering in Duckweed ( Lemna gibba )

: Long noncoding RNAs (lncRNAs) are crucial components in regulating the ﬂowering of plants. However, the regulatory mechanism of lncRNAs underlying salicylic acid (SA)-induced ﬂowering remains unknown in duckweed (e.g., Lemna gibba L.), an aquatic model species with signiﬁcant potential applications in agriculture and industry. In this work, L. gibba plants were collected at four crucial time points during SA-induced ﬂowering and subjected to PacBio full-length sequencing and strand-speciﬁc RNA sequencing. A total of 474 lncRNAs were identiﬁed, of which 31 were differentially expressed and involved in SA-induced ﬂowering. A trans -regulatory analysis found that these lncRNAs displayed temporal-speciﬁc expression trends and mainly participated in stress metabolism, photosynthesis, jasmonate metabolism, and transport under SA treatment. Five lncRNAs were determined to act as targets of miRNAs that played critical roles in regulating ﬂowering. In addition, ﬁfteen lncRNAs showed co-expression with ﬂowering-related genes, and lncRNA03 and lncRNA25 were identiﬁed as key players involved in ﬂowering via lncRNA-miRNA-mRNA interactions. Finally, twelve lncRNAs related to trans -regulation, miRNA targets, or co-expression with ﬂowering-related genes were veriﬁed by qRT-PCR. These ﬁndings deepen our understanding of lncRNAs in SA-induced ﬂowering in duckweed and provide valuable resources for in-depth functional analysis in the future.


Introduction
Long noncoding RNAs (lncRNAs) are abundant in plants and commonly defined as RNA transcripts longer than 200 nucleotides in length but lacking protein-coding ability.Emerging studies have uncovered that lncRNAs are crucial regulators participating in the stress response [1,2], vernalization [3], and flowering [4].The lncRNA VAS (TaVRN1 alternative splicing) was vernalization-induced, and overexpression of VAS promoted TaVRN1 expression to accelerate flowering in winter wheat [3].The lncRNA COOLAIR was cold-induced and acted as a key repressor of flowering during vernalization by regulating the expression of FLC (FLOWERING LOCUS C) [4].
LncRNAs carry out their biological functions through various mechanisms.They can function in a cis-acting or trans-acting manner to regulate gene expression by interacting with RNAs or DNAs through sequence complementarity, epigenetic modifications, or the regulation of promoter activities [5,6].In cis-acting regulation, the flowering-involved lncRNA MAS interacted with WDR5a (one core component of the COMPASS-like complexes) and activated the expression of MAF4 (MADS AFFECTING FLOWERING4) [7]; lncRNA33732 was positioned down-stream of RBOH (respiratory burst oxidase) and enhanced the expression level of RBOH, leading to increased H 2 O 2 levels in the defense reaction of tomato [8].In trans-acting regulation, the lncRNA TAS3 (trans-acting siRNA3) repressed the transcription of NRT2.4 (nitrate transporter 2.4), influencing nitrogen transport and root development in low-nitrogen environments [9]; the lncRNA FLAIL (floweringassociated intergenic lncRNA) directly bound to target genes that controlled flowering via RNA-DNA interactions [10].LncRNAs can also serve as targets for cleavage by miR-NAs [11,12].Additionally, lncRNAs have recently been found to act as target mimics, affecting the interactions of miRNAs and mRNAs by combining with miRNAs [13].For instance, the lncRNA IPS1 (INDUCED BY PHOSPHATE STARVATION 1) increased the transcript accumulation of miR399-mediated PHO2 (PHOSPHATE 2) by functioning as a target mimic [14].
Massive sequencing technology has led to the discovery of numerous lncRNAs associated with flowering in plants.In Arabidopsis, 2132 lncRNAs were identified by RNA-seq analysis, of which the lncRNA AtLnc428 functioned in temperature-mediated flowering [15].In tomato, 248 novel lncRNAs were predicted by transcriptomic analysis, and 65 were significantly differentially expressed during flowering [16].In chickpea, 1414 lncRNAs were identified by RNA-seq profiling, and 37 could bind to flowering-related genes to maintain inflorescence development [17].In Angelica sinensis, 2327 lncRNAs were identified during vernalization, and 36 were characterized as participants in flowering [18].In Chinese cabbage, 3382 lncRNAs were identified in response to plumule vernalization, and a few of them were involved in promoting flowering [19].In Prunus sibirica, 3423 lncRNAs were identified by whole-transcriptome RNA sequencing, and 1351 were associated with flowering time [20].These findings indicated significant roles of lncRNAs in the flowering of plants.However, no comprehensive analyses have been reported regarding flowering-related lncRNAs in duckweed.
Duckweed is a family of small floating aquatic plants with rapid growth, high starch content, and efficient nutrient uptake.Thus, it has shown great potential applications in animal feed, bioethanol production, and sewage treatment [21][22][23].Duckweed propagates mainly asexually and rarely flowers in natural conditions, which dramatically limits its germplasm collection and heterosis usage [24].Salicylic acid (SA) is critical in duckweed species' flower induction [25][26][27].SA is also a key regulator in the flowering of Arabidopsis, Pharbitis nil, and Malus domestica [28][29][30].Although hundreds of genes and tens of functional pathways (e.g., stress, carbohydrate metabolism, and hormone signaling) have been identified, the roles of lncRNAs remain largely elusive in SA-induced flowering.
In this work, L. gibba plants were collected at four crucial time points during SAinduced flowering and subjected to PacBio full-length sequencing and strand-specific RNA sequencing (ssRNA-seq).Afterward, the flowering-related lncRNAs were systematically identified, their expression patterns were investigated, and their putative functions were analyzed.These findings will expand our knowledge of lncRNAs in SA-induced flowering in duckweed and provide valuable resources for upcoming functional investigation of these lncRNAs.

Plant Materials and SA Treatment
Lemna gibba clone 7741 was used as plant material in this work.It was derived from the Biological Collection Centre of the Institute of Tropical Bioscience and Biotechnology, Chinese Academy of Tropical Agricultural Sciences, Hainan, China.As described before [27], three-frond colonies were cultivated in 250 mL flasks containing 100 mL of MH medium (pH = 5.9) supplemented with 20 µM SA and grown under 16 h light vs. 8 h dark fluctuations at 40 µmol•m −2 •s −1 photosynthetic active radiation and 25 ± 1 • C. Colonies cultivated on a medium without SA supplementation were used as the control.
Based on our previous studies [26,27], day 0 (D0, the control), day 12 (D12, one day before the day of first flowering), day 13 (D13, the day of first flowering), and day 20 (D20, the day with highest flowering rate) represented four critical time points during SA-induced flowering.Therefore, L. gibba colonies were collected with three replicates at these time points.The samples were frozen immediately in liquid nitrogen and kept at −80 • C until used for full-length transcriptome sequencing and ssRNA-seq.

Full-Length Transcriptome Analysis
Full-length transcriptome sequencing was carried out by the Annoroad Gene Technology Corporation (Beijing, China).Equal amounts of total RNA from D0, D12, D13, and D20 were mixed and used for PacBio library construction using the SMRTbell Template Prep Kit (Pacific Biosciences, Menlo Park, CA, USA).PacBio sequencing was performed on a PacBio sequel platform.
Full-length transcripts were obtained with the IsoSeq3 software (https://github.com/PacificBiosciences/IsoSeq3) (accessed on 6 September 2021).Circular consensus sequence (CCS) reads were derived from the sub-reads with the removal of primers and unwanted combinations.Full-length non-concatemer (FLNC) transcripts were generated by removing polyA tails and artificial concatemers.Finally, these FLNC transcripts were clustered and polished to produce consensus sequences.The generated high-quality, full-length transcripts were then used for ssRNA-seq analysis and lncRNA identification.

Library Preparation and ssRNA-Seq Analysis
The ssRNA-seq libraries were constructed and sequenced by the Annoroad Gene Technology Corporation (Beijing, China).The quality and concentration of total RNA were assessed using the Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) and the Agilent 2100 Bioanalyzer system (Agilent, Palo Alto, CA, USA).The ssRNA-seq library was prepared using the TruSeq TM RNA sample prep Kit (Illumina, San Diego, CA, USA) and subjected to rRNA removal using the Ribo-Zero Magnetic kit.These libraries were sequenced on the Illumina Hiseq-4000 platform to produce 150 bp paired-end reads.Each sample was analyzed with three biological replicates.
Clean reads were mapped to full-length transcripts for gene expression quantification.Reads Per Kilobase per Million mapped reads (RPKM) were calculated to estimate the expression level [27].Differentially expressed (DE) genes were identified by DESeq2 [31] setting |log2FC (fold-change)| > 1 and a false discovery rate (FDR) < 0.05.

LncRNA Identification
LncRNAs were identified via a vigorous pipeline, as previously described [32,33].Firstly, transcripts with a length less than 200 bp, an ORF length greater than 300, and read coverage below two were eliminated.Secondly, transcripts with potential protein-coding ability were discarded according to an assessment using the Coding Potential Calculator (CPC) [34] and Coding/Noncoding Index (CNCI) [35].Thirdly, the transcripts containing Pfam domains were also eliminated according to hidden Markov models.DE lncRNAs were identified using DESeq2 [31] with parameters of FDR < 0.05 and |log2FC| > 1.

Prediction of LncRNA Targets
To determine lncRNA targets in trans-acting regulation, DE genes and DE lncRNAs were merged into an expression matrix to perform a co-expression network analysis by WGCNA [36].Genes and lncRNAs displaying similar expression patterns were clustered in the same group and might be involved in trans-regulation.
To explore the biological function of lncRNAs, L. gibba full-length transcripts were annotated and assigned to different functional categories by using the MapMan tool [37].The significance of functional enrichment was calculated by Fisher's exact test, as previously reported [38].

Identification of LncRNAs Functioning as MiRNA Targets
The DE lncRNAs discovered in this work and the miRNAs derived from miRBase (Release 22) were uploaded to PsRobot [39] and TAPIR [40] for the identification of lncR-NAs acting as miRNA targets.The interactive networks of lncRNAs and miRNAs were visualized by Cytoscape [41].

Identification of LncRNAs Involved in Flowering
The sequences of DE genes were used to perform a BLASTP search against the Arabidopsis FLOR-ID database [42] to determine flowering-related genes in L. gibba.The DE lncRNAs and flowering-related genes were clustered by Pearson's correlation of their expression profiles.Those lncRNAs displaying high similarity (r > 0.8) to the expression of flowering-related genes were considered to be involved in flowering.The interactive networks of lncRNAs and genes were also visualized by Cytoscape [41].

qRT-PCR Analysis
Total RNA was extracted from L. gibba plants using the RNA-prep Pure Plant Kit (TIANGEN Biotech, China).The first strand of cDNA was reverse-transcribed using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China).To confirm the expression levels derived from ssRNA-seq, a total of twelve lncRNAs were investigated by qRT-PCR using SYBR-green (TaKaRa, Dalian, China) on the Stratagene Mx3005P machine (Stratagene, San Diego, CA, USA).The qRT-PCR procedures were as follows: 95 • C for 30 s, followed by 40 cycles of 95 • C for 10 s, 60 • C for 20 s, and 72 • C for 20 s.The actin gene was applied as an internal control (Supplemental Table S1).qRT-PCR was carried out with three biological replicates, and the relative expression levels were calculated by the 2 −∆∆Ct method [27].

Full-Length Transcriptome and LncRNA Identification
L. gibba colonies were induced to first flower by SA treatment on day 13 (D13), and the flowering rate reached the maximum value on D20 (Figure 1A).Therefore, bulked RNAs were isolated from samples taken at four crucial time points (D0, D12, D13, and D20) and then sequenced by the PacBio Iso-Seq platform to obtain wide transcriptome coverage.In total, 737,364 CCS reads were generated with a mean length of 1470 bp.By removing primers and polyA tails, these CCS reads were processed into 701,677 full-length nonchimeric reads, which were further clustered and polished to produce 42,281 high-quality, full-length transcripts (mean length = 1405 bp) for lncRNA identification.
In total, 987, 734, and 1123 lncRNAs were identified using the CPC, CNCI, and Pfam tools, respectively (Figure 1B).Of these, 474 were detected with all three tools.About onehalf and 30% of these lncRNAs were 201-600 bp and 601-1000 bp in length, with means of 400 bp and 780 bp (Figure 1C), respectively.The percentage of expressed lncRNAs (RPKM > 1) gradually increased from D0 to D20 (Figure 1D).Subsequently, the 474 lncRNAs were used to analyze their transcriptional responses under SA treatment.

ssRNA-Seq and DE LncRNA Identification
The samples collected on D0, D12, D13, and D20 under SA treatment were also used to perform ssRNA-seq analysis.By removing low-quality reads and adapter sequences, approximately 850 million clean reads were obtained, and 84% of them were mapped to full-length transcripts for expression quantitation.

Functional Analysis of LncRNAs in Trans-Acting Regulation
To reveal the putative roles of lncRNAs in trans-acting regulation, 1897 DE genes and 31 DE lncRNAs were applied to perform co-expression network analysis.In total, nine coexpression groups/modules (M1-M9) were detected based on the expression patterns (Figure 2A).

ssRNA-Seq and DE LncRNA Identification
The samples collected on D0, D12, D13, and D20 under SA treatment were also used to perform ssRNA-seq analysis.By removing low-quality reads and adapter sequences, approximately 850 million clean reads were obtained, and 84% of them were mapped to full-length transcripts for expression quantitation.

Functional Analysis of LncRNAs in Trans-Acting Regulation
To reveal the putative roles of lncRNAs in trans-acting regulation, 1897 DE genes and 31 DE lncRNAs were applied to perform co-expression network analysis.In total, nine co-expression groups/modules (M1-M9) were detected based on the expression patterns (Figure 2A).Group M1 contained five lncRNAs.The expression of lncRNAs and genes in this group was significantly depressed from D0 to D20 under SA treatment.Category enrichment revealed that those lncRNAs and genes were enriched during stress (Figure 2B).RD22 (related to desiccation and abscisic acid), PRB1 (related to ethylene and SA responses), and MLO1 and MLO5 (related to cell death and the defense response) were included in M1 (Supplemental Table S2).Many heat shock proteins, including HSP70, CR88, CLPB3, HSP91, HSP90.1, and HSP101h, were also included in this group.
There were three lncRNAs in group M3.The expression of lncRNAs and genes gradually increased on D12 and D13 but suddenly declined on D20.These lncRNAs and genes were enriched in the Calvin cycle, light reaction, and nitrate transport.FBA2, GAPA2, and several rubisco small subunits involved in the Calvin cycle were included in M3.ATPC1 (involved in ATP synthase), PGR5 (involved in electron flow), and a few subunits related to photosystem I (LHCA1, PSAG, PSAE2) and photosystem II (LHB1B1 and PPL2) were also included (Supplemental Table S2).
There were nine lncRNAs in group M6, in which the lncRNAs and genes showed the highest expression levels on D20.The enriched categories of this group included Group M1 contained five lncRNAs.The expression of lncRNAs and genes in this group was significantly depressed from D0 to D20 under SA treatment.Category enrichment revealed that those lncRNAs and genes were enriched during stress (Figure 2B).RD22 (related to desiccation and abscisic acid), PRB1 (related to ethylene and SA responses), and MLO1 and MLO5 (related to cell death and the defense response) were included in M1 (Supplemental Table S2).Many heat shock proteins, including HSP70, CR88, CLPB3, HSP91, HSP90.1, and HSP101h, were also included in this group.
There were three lncRNAs in group M3.The expression of lncRNAs and genes gradually increased on D12 and D13 but suddenly declined on D20.These lncRNAs and genes were enriched in the Calvin cycle, light reaction, and nitrate transport.FBA2, GAPA2, and several rubisco small subunits involved in the Calvin cycle were included in M3.ATPC1 (involved in ATP synthase), PGR5 (involved in electron flow), and a few subunits related to photosystem I (LHCA1, PSAG, PSAE2) and photosystem II (LHB1B1 and PPL2) were also included (Supplemental Table S2).
There were nine lncRNAs in group M6, in which the lncRNAs and genes showed the highest expression levels on D20.The enriched categories of this group included jasmonate metabolism, signaling, phosphate transport, and potassium transport (Figure 2B).LOX1 and seven LOX5 homologs related to jasmonate biosynthesis were included in M6.Many genes related to phosphate transport (PHT1;2a, PHT1;2b, PHT1;4a, PHT1;4b, and PHT1;4c) and potassium transport (HAK5a to HAK5f, and GORK) were also included (Supplemental Table S2).
There were nine lncRNAs in group M7.These lncRNAs and genes showed high expression levels on D13 and D20, and they were enriched in tetrapyrrole synthesis.A few genes (including HEMA1, GATB, CRD1, UPM1, HEME1, and HEME2) related to tetrapyrrole synthesis were included in this group (Supplemental Table S2).
In contrast, no lncRNAs were found in M2, M8, or M9, although the genes from these groups were enriched in the light reaction and wax metabolism, respectively.
Together, these results suggest that the DE lncRNAs mainly participate in stress metabolism, photosynthesis (including the Calvin cycle, light reaction, and photorespiration), tetrapyrrole synthesis, jasmonate metabolism, and transport via trans-acting regulation under SA treatment.

Functional Analysis of LncRNAs as MiRNA Targets
LncRNAs can influence the regulatory efficiency of miRNAs by acting as their competitive targets.Thus, it is necessary to investigate the roles of DE lncRNAs in functioning as miRNA targets.
Five lncRNAs were identified as targets of ten miRNAs (Supplemental Table S3).MiR169 is crucial in promoting early flowering in response to various abiotic stresses [43].In this work, lncRNA01 was targeted by six members of the miR169 family (Figure 3A), indicating that lncRNA01 might be involved in SA-induced flowering of L. gibba via the interaction with miR169.
In contrast, no lncRNAs were found in M2, M8, or M9, although the genes from these groups were enriched in the light reaction and wax metabolism, respectively.
Together, these results suggest that the DE lncRNAs mainly participate in stress metabolism, photosynthesis (including the Calvin cycle, light reaction, and photorespiration), tetrapyrrole synthesis, jasmonate metabolism, and transport via transacting regulation under SA treatment.

Functional Analysis of LncRNAs as MiRNA Targets
LncRNAs can influence the regulatory efficiency of miRNAs by acting as their competitive targets.Thus, it is necessary to investigate the roles of DE lncRNAs in functioning as miRNA targets.
Five lncRNAs were identified as targets of ten miRNAs (Supplemental Table S3).MiR169 is crucial in promoting early flowering in response to various abiotic stresses [43].In this work, lncRNA01 was targeted by six members of the miR169 family (Figure 3A), indicating that lncRNA01 might be involved in SA-induced flowering of L. gibba via the interaction with miR169.MiR156 influences plant flowering by controlling the transition from the vegetative to the reproductive phase [44], while miR396 is involved in the plant response to vernalization and flower development [45].Interestingly, lncRNA03 was targeted by both miR396 and miR156i (Figure 3B), strongly suggesting the involvement of lncRNA03 in SA-induced flowering through miRNA-lncRNA interactions.MiR156 influences plant flowering by controlling the transition from the vegetative to the reproductive phase [44], while miR396 is involved in the plant response to vernalization and flower development [45].Interestingly, lncRNA03 was targeted by both miR396 and miR156i (Figure 3B), strongly suggesting the involvement of lncRNA03 in SA-induced flowering through miRNA-lncRNA interactions.
Similarly, lncRNA07 and lncRNA09 were targeted by miR156h from the miR156 family (Figure 3C,D).In addition, these two lncRNAs and lncRNA25 were targeted by miR5021, which is pollen-specific and involved in abiotic stress [46,47], indicating that these lncRNAs might function in stress-induced flowering.
Collectively, these findings revealed possible roles of lncRNAs in SA-induced flowering of L. gibba by acting as miRNA targets.

Identification of LncRNAs Involved in Flowering
In total, 15 flowering-related DE genes were found by BLASTP against the Arabidopsis FLOR-ID database.Then, these genes and DE lncRNAs were used to perform a coexpression analysis to identify lncRNAs involved in flowering.
LncRNA01, lncRNA06, lncRNA16, lncRNA17, and lncRNA25 were co-expressed with eight flowering-related genes (Figure 4A and Supplemental Table S4).Their expression levels were significantly induced by SA treatment on D12 but were then depressed on D13 and D20 (Figure 4B).LHY, lncRNA06, and CCA1a were the top three hub genes with the most connections to others.In addition, lncRNA06, lncRNA01, lncRNA17, and lncRNA25 showed connections to as many as four flowering-related genes, indicating their crucial roles in SA-induced flowering.
Agronomy 2023, 13, x FOR PEER REVIEW 8 of 14 Similarly, lncRNA07 and lncRNA09 were targeted by miR156h from the miR156 family (Figure 3C,D).In addition, these two lncRNAs and lncRNA25 were targeted by miR5021, which is pollen-specific and involved in abiotic stress [46,47], indicating that these lncRNAs might function in stress-induced flowering.
Collectively, these findings revealed possible roles of lncRNAs in SA-induced flowering of L. gibba by acting as miRNA targets.

Identification of LncRNAs Involved in Flowering
In total, 15 flowering-related DE genes were found by BLASTP against the Arabidopsis FLOR-ID database.Then, these genes and DE lncRNAs were used to perform a coexpression analysis to identify lncRNAs involved in flowering.
LncRNA01, lncRNA06, lncRNA16, lncRNA17, and lncRNA25 were co-expressed with eight flowering-related genes (Figure 4A and Supplemental Table S4).Their expression levels were significantly induced by SA treatment on D12 but were then depressed on D13 and D20 (Figure 4B).LHY, lncRNA06, and CCA1a were the top three hub genes with the most connections to others.In addition, lncRNA06, lncRNA01, lncRNA17, and lncRNA25 showed connections to as many as four flowering-related genes, indicating their crucial roles in SA-induced flowering.LncRNA22 and lncRNA26 were co-expressed with three flowering-related genes, namely, CCA1b, STO, and GLK1 (Figure 4C).The expression of these genes and lncRNAs was significantly induced by SA treatment on D13 but was depressed on D20 (Figure 4D).
Together, these findings revealed a dozen flowering-involved lncRNAs, which tended to function at different stages during SA-induced flowering of L. gibba.

Expression Verification of LncRNAs
Twelve lncRNAs related to trans-regulation, miRNA targets, or co-expression with flowering-related genes were examined by qRT-PCR.The correlation coefficients ranged from 0.78 to 0.99 between the two methods of qRT-PCR and ssRNA-seq (Figure 5 and Supplemental Table S1), supporting reliable expression levels of lncRNAs investigated by ssRNA-seq.LncRNA22 and lncRNA26 were co-expressed with three flowering-related genes, namely, CCA1b, STO, and GLK1 (Figure 4C).The expression of these genes and lncRNAs was significantly induced by SA treatment on D13 but was depressed on D20 (Figure 4D).
Together, these findings revealed a dozen flowering-involved lncRNAs, which tended to function at different stages during SA-induced flowering of L. gibba.

Expression Verification of LncRNAs
Twelve lncRNAs related to trans-regulation, miRNA targets, or co-expression with flowering-related genes were examined by qRT-PCR.The correlation coefficients ranged from 0.78 to 0.99 between the two methods of qRT-PCR and ssRNA-seq (Figure 5 and Supplemental Table S1), supporting reliable expression levels of lncRNAs investigated by ssRNA-seq.LncRNAs have been revealed to play important regulatory roles in the flowering of many species, including Arabidopsis [4], wheat [3], Lactuca sativa [48], Chinese cabbage [19], and chickpea [17].However, the roles of lncRNAs in flowering have not been reported so far in duckweed, which dramatically limits its germplasm collection and heterosis usage [26,27].In this study, a total of 474 lncRNAs were systematically identified in L. gibba during SA-induced flowering.The lncRNA number was approximately 1.5-fold greater than that identified in L. sativa [48] but was less than that reported in A. sinensis, P. sibirica, and chickpea [17,18,20].To some extent, these findings demonstrate the impact of the species, sequenced depth, and used criteria on the identification of lncRNAs.
LncRNAs can function in a temporal-specific or tissue-specific manner [32,33].In this work, 31 lncRNAs were differentially expressed during SA-induced flowering.Of these, five, five, and four were, respectively, identified at flower induction (D12/D0), flower initiation (D13/D12), and flowering stages (D20/D13) [27], but none of them were found to be common among these stages (Figure 1F).Moreover, the trans-regulatory co-expression analysis further supported the regulation of lncRNA expression via a temporal-specific pattern (Figure 2A), consistent with those previously described [32,33].Collectively, these findings strongly suggest that lncRNAs are critical players involved in the flowering of L. gibba.This study is the first report on SA-responsive lncRNAs in duckweed and provides a new angle to explore the mechanisms of SA-induced flowering.

Functional Analysis of LncRNAs during SA-Induced Flowering
LncRNAs can exert their functions in a trans-acting manner by influencing the transcript level of genes positioned at distant genomic locations [9,10].In this study, coexpression networks were conducted to survey the functions of lncRNAs according to the functional enrichment of the co-expressed genes.
The expression of five lncRNAs and 310 co-expressed genes in M1 significantly decreased from D0 to D20, in accord with the inhibition of L. gibba growth under SA treatment [27].Moreover, the 310 co-expressed genes were enriched in stress metabolism (Figure 2B).Given that L. gibba flowering induced by SA was associated with stress [27], our results highly suggest that these lncRNAs might regulate flowering via stress-related pathways.The results are partially supported by a recent study that reported that lncR-NAs interacted with stress-related target genes in response to SA treatment [49].The lncRNAs also showed co-expression with genes relevant to photosynthesis, tetrapyrrole synthesis, jasmonate metabolism, signaling, and transport (Figure 2B), indicating similar roles of lncRNAs in L. gibba during SA-induced flowering.LncRNAs have been evidenced to play comparable roles in the flowering of other species: in A. sinensis, Liu et al. [18] found the co-expression of lncRNAs with genes related to carbohydrate metabolism, biosignaling, energy, and transport during vernalization; in P. sibirica, Xu et al. [20] revealed that lncRNAs participated in secondary metabolism, hormone signaling, and starch and sucrose metabolism during flowering; in L. sativa and chickpea, lncRNAs were reported to participate in flowering based on their co-expression with protein-coding genes related to photoperiodism, transport, and hormone signaling [17,48].Most importantly, fifteen lncRNAs were co-expressed with flowering-related genes (including LHY, CCA1a, PI, and SVP) (Figure 4), providing us with excellent hints to further characterize the function of lncRNAs by examining their interactions [4,10].
LncRNAs can also exert their roles in cis-acting regulation to influence the transcript levels of neighboring genes [7,8].However, no such analyses were performed in this work due to the unavailability of a reference genome for L. gibba.

Interactive Networks of lncRNA-mRNA-miRNA in SA-Induced Flowering
In addition, lncRNAs can function as miRNA targets, influencing miRNA-mediated gene regulation [12].Numerous plant miRNAs have been well studied; therefore, their functions can be utilized to investigate the potential roles of lncRNAs.In this study, five DE lncRNAs were found to be targets of ten miRNAs (Figure 3).LncRNA25 was targeted by miR5021, which is pollen-specific and involved in abiotic stress [46,47].Four floweringrelated genes (PI, ATJ3, CCA1a, and SVP), which were co-expressed with lncRNA25, were also targeted by miR5021 (Figure 6A and Supplemental Table S5).These results suggest that lncRNA25 participates in flowering by mediating miR5021 to regulate the expression of these flowering-related genes.Recently, it has become evident that flowering is also influenced by stress, which is recognized as the third category of the flowering response, in addition to vernalization and photoperiod flowering [50].Under stress conditions, genes related to photosynthesis and energy metabolism were dramatically influenced [27].Consistently, a few genes related to stress (HSP70.2,HSP90.1, and NCED4), photosynthesis (RCA1, RCA2, and DRT112), and the TCA cycle (CA1 and MTE2) were also identified as targets of miR5021 (Figure 6A and Supplemental Table S5).Meanwhile, these genes were coexpressed with lncRNA25, indicating that the flowering mediated by lncRNA25 interacting with miR5021 was stress-involved.This result is consistent with the involvement of stress in the SA-induced flowering of L. gibba [27].
gene regulation [12].Numerous plant miRNAs have been well studied; therefore, their functions can be utilized to investigate the potential roles of lncRNAs.In this study, five DE lncRNAs were found to be targets of ten miRNAs (Figure 3).LncRNA25 was targeted by miR5021, which is pollen-specific and involved in abiotic stress [46,47].Four floweringrelated genes (PI, ATJ3, CCA1a, and SVP), which were co-expressed with lncRNA25, were also targeted by miR5021 (Figure 6A and Supplemental Table S5).These results suggest that lncRNA25 participates in flowering by mediating miR5021 to regulate the expression of these flowering-related genes.Recently, it has become evident that flowering is also influenced by stress, which is recognized as the third category of the flowering response, in addition to vernalization and photoperiod flowering [50].Under stress conditions, genes related to photosynthesis and energy metabolism were dramatically influenced [27].Consistently, a few genes related to stress (HSP70.2,HSP90.1, and NCED4), photosynthesis (RCA1, RCA2, and DRT112), and the TCA cycle (CA1 and MTE2) were also identified as targets of miR5021 (Figure 6A and Supplemental Table S5).Meanwhile, these genes were co-expressed with lncRNA25, indicating that the flowering mediated by lncRNA25 interacting with miR5021 was stress-involved.This result is consistent with the involvement of stress in the SA-induced flowering of L. gibba [27].Similarly, lncRNA03 was targeted by miR156i, which was reported to play a major role in flowering by controlling the transition from the vegetative to the reproductive phase [44].MiR156 also targeted SPL members and activated the expression of floweringrelated genes (e.g., LFY and AP1) to promote flowering [51,52].Five flowering-related genes (AP1, FT1, GLK1, GRP7, and TEM1a) co-expressing with lncRNA03 were included in this co-expression network.Of these, GLK1 was targeted by miR156i (Figure 6B).Several genes related to stress (ABI1, UVR8, and PAL1) and starch metabolism (BMY3, APE2, and ADG1) were co-expressed with lncRNA03 and identified as targets of miR156i as well (Figure 6B and Supplemental Table S5).These results suggest that lncRNA03 participates in flowering by mediating miR156i to regulate the expression of flowering- Similarly, lncRNA03 was targeted by miR156i, which was reported to play a major role in flowering by controlling the transition from the vegetative to the reproductive phase [44].MiR156 also targeted SPL members and activated the expression of flowering-related genes (e.g., LFY and AP1) to promote flowering [51,52].Five flowering-related genes (AP1, FT1, GLK1, GRP7, and TEM1a) co-expressing with lncRNA03 were included in this co-expression network.Of these, GLK1 was targeted by miR156i (Figure 6B).Several genes related to stress (ABI1, UVR8, and PAL1) and starch metabolism (BMY3, APE2, and ADG1) were co-expressed with lncRNA03 and identified as targets of miR156i as well (Figure 6B and Supplemental Table S5).These results suggest that lncRNA03 participates in flowering by mediating miR156i to regulate the expression of flowering-related genes.Moreover, this flowering process was accompanied by the expression regulation of genes relevant to stress and carbohydrate metabolism, in accord with previous studies [50,53].
Together, these findings unveiled complex networks of lncRNAs, mRNAs, and miR-NAs in SA-induced flowering of L. gibba, leading to the proposal of a schematic model (Figure 6C).However, the functions of lncRNAs have not been verified by experiments, and further studies are required to explore their roles in depth.

Conclusions
In this work, a total of 474 lncRNAs were systematically identified in L. gibba, of which 31 were differentially expressed during SA-induced flowering.The potential roles of lncRNAs were extensively analyzed for their roles in trans-acting regulation, miRNA targets, and co-expression with flowering-related genes.Notably, lncRNA03 and lncRNA25 were identified as key players involved in flowering via lncRNA-miRNA-mRNA interactions.The findings provide valuable resources for the in-depth functional analysis of lncRNAs participating in the SA-induced flowering of duckweed.In the future, we will explore the regulatory mechanisms of lncRNAs in duckweed flowering.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/agronomy13102631/s1:Table S1.Primer sequences and the association between qRT-PCR and ssRNA-seq data.Table S2.Summary of DE lncRNAs and their coexpressed genes in response to SA-induced flowering.Table S3.Information on lncRNAs functioning as miRNA targets.Table S4.Summary of lncRNAs and the co-expressed flowering-related genes in Figure 4. Table S5.Information on genes in the lncRNA-mRNA-miRNA networks in Figure 6.
Author Contributions: J.Z., Z.D. and L.F. conceptualized the study; L.F., D.T. and X.S. performed the experiments and analyzed the data; L.F., J.Z. and Z.D. were responsible for writing-original draft preparation and reviewing and editing.All authors have read and agreed to the published version of the manuscript.
Funding: This project was funded by the Natural Science Foundation of Hainan Province (322MS123), the High Level Talents Project of Hainan Basic and Applied Research Program (321RC639), and the Natural Science Foundation of China (31800301).

Figure 2 .
Figure 2. Functional prediction of lncRNAs by co-expression analysis.(A) The expression profiles of DE lncRNAs and mRNAs.A total of nine co-expressed groups (M1-M9) were detected.(B) Category enrichment of the groups presented in (A).

Figure 2 .
Figure 2. Functional prediction of lncRNAs by co-expression analysis.(A) The expression profiles of DE lncRNAs and mRNAs.A total of nine co-expressed groups (M1-M9) were detected.(B) Category enrichment of the groups presented in (A).

Figure 3 .
Figure 3. Analyses of lncRNAs functioning as miRNA targets.(A-C) Interactive networks of miRNAs and lncRNAs.The miRNAs are indicated by blue triangles, while the lncRNAs are indicated by pink circles.The node sizes represent connections to others.(D) Examples of sequence alignment between lncRNAs and miRNAs in (A-C).

Figure 3 .
Figure 3. Analyses of lncRNAs functioning as miRNA targets.(A-C) Interactive networks of miRNAs and lncRNAs.The miRNAs are indicated by blue triangles, while the lncRNAs are indicated by pink circles.The node sizes represent connections to others.(D) Examples of sequence alignment between lncRNAs and miRNAs in (A-C).

Figure 4 .
Figure 4. Investigation of lncRNAs co-expressed with flowering-related genes.(A,C,E) Interactive networks of lncRNAs and the co-expressed flowering-related genes.LncRNAs and genes are represented by pink circles and green diamonds, respectively.The node sizes represent connections to others.(B,D,F) The expression profiles of lncRNAs and flowering-related genes in (A,C,E).

Figure 4 .
Figure 4. Investigation of lncRNAs co-expressed with flowering-related genes.(A,C,E) Interactive networks of lncRNAs and the co-expressed flowering-related genes.LncRNAs and genes are represented by pink circles and green diamonds, respectively.The node sizes represent connections to others.(B,D,F) The expression profiles of lncRNAs and flowering-related genes in (A,C,E).

Figure 5 .
Figure 5. Expression verification of twelve lncRNAs.Expression levels were normalized among samples.The data are shown as mean ± standard deviation from three biological replications.The blue and red lines represent the results derived from qRT-PCR and ssRNA-seq, respectively.

Figure 5 .
Figure 5. Expression verification of twelve lncRNAs.Expression levels were normalized among samples.The data are shown as mean ± standard deviation from three biological replications.The blue and red lines represent the results derived from qRT-PCR and ssRNA-seq, respectively.

4 .
Discussion 4.1.LncRNA Is a Critical Player Involved in the Flowering of L. gibba

Figure 6 .
Figure 6.Network analysis of lncRNAs, mRNAs, and miRNAs.(A,B) Representative interaction networks of lncRNAs, mRNAs, and miRNAs relevant to SA-induced flowering of L. gibba.LncRNAs, mRNAs, and miRNAs are indicated by circles, diamonds, and triangles, respectively.mRNAs assigned to different functional categories are marked by distinct colors.The node sizes represent connections to others.(C) A schematic diagram of lncRNA-mRNA-miRNA networks for SA-induced flowering of L. gibba.

Figure 6 .
Figure 6.Network analysis of lncRNAs, mRNAs, and miRNAs.(A,B) Representative interaction networks of lncRNAs, mRNAs, and miRNAs relevant to SA-induced flowering of L. gibba.LncRNAs, mRNAs, and miRNAs are indicated by circles, diamonds, and triangles, respectively.mRNAs assigned to different functional categories are marked by distinct colors.The node sizes represent connections to others.(C) A schematic diagram of lncRNA-mRNA-miRNA networks for SA-induced flowering of L. gibba.