Genome-Wide Analysis of Long Non-coding RNAs Involved in Nodule Senescence in Medicago truncatula

Plant long non-coding RNAs (lncRNAs) are widely accepted to play crucial roles during diverse biological processes. In recent years, thousands of lncRNAs related to the establishment of symbiosis, root nodule organogenesis and nodule development have been identified in legumes. However, lncRNAs involved in nodule senescence have not been reported. In this study, senescence-related lncRNAs were investigated in Medicago truncatula nodules by high-throughput strand-specific RNA-seq. A total of 4576 lncRNAs and 126 differentially expressed lncRNAs (DElncRNAs) were identified. We found that more than 60% lncRNAs were associated with transposable elements, especially TIR/Mutator and Helitron DNA transposons families. In addition, 49 DElncRNAs were predicted to be the targets of micro RNAs. Functional analysis showed that the largest sub-set of differently expressed target genes of DElncRNAs were associated with the membrane component. Of these, nearly half genes were related to material transport, suggesting that an important function of DElncRNAs during nodule senescence is the regulation of substance transport across membranes. Our findings will be helpful for understanding the functions of lncRNAs in nodule senescence and provide candidate lncRNAs for further research.

Root nodules are special organs formed by legume-rhizobium symbiosis. Emerging evidence suggests that lncRNAs function as crucial regulators of symbiotic nitrogen fixation (SNF) in nodules. A well-known lncRNA associated with SNF is ENOD40 in M. truncatula (Campalans et al., 2004) which can act as a dual RNA (Bardou et al., 2011) in nodule organogenesis. Another lncRNA is TAS3 RNA in M. truncatula and the miR390/TAS3 pathway plays negative roles in nodulation and nodule development (Hobecker et al., 2017). Recently, thousands of lncRNAs in M. truncatula have been identified to be involved in SNF and possibly regulate mRNA expression in cis way (Pecrix et al., 2018). SNF by nodules lasts for a peroid, peaks at some time in the nodule life-span and declines with the senescence of nodules. Mature indeterminate nodules (such as nodules on M. trunctula roots) are divided into four developmental regions namely the apical meristematic, the infection, the nitrogen fixation and the senescence zones ( Van de Velde et al., 2006). Nodule senescence is a developmental process which is initiated in the senescence zone and advances gradually to the meristematic zone. Although a large number of lncRNAs involved in SNF have been identified, little is known about the lncRNAs related to nodule senescence. Interestingly, several recent reports have suggested that lncRNAs play key roles in leaf senescence (Chao et al., 2019;Huang et al., 2021) and nodule senescence has a relatively high similarity with leaf senescence at the molecular level ( Van de Velde et al., 2006), indicating that lncRNAs are also likely to be important regulators during nodule senescence. However, previous work focused on the identification and functions of multiple genes involved in nodule senescence, the research on ageing-related lncRNAs in root nodules has been lacking. In this study, we conducted high-throughput strand-specific RNA-seq of nodules at 21-and 35-days post inoculation (dpi) with Sinorhizobium meliloti 1021 to investigate and characterize lncRNAs associated with nodule senescence. Our findings will provide new insights into the underlying functions of lncRNAs during nodule senescence.

Plant Materials
M. truncatula A17 seeds were surface-sterilized in 75% ethanol for 5 min and 2% sodium hypochlorite solution for 15 min, before washing 5-6 times with sterile water. The seeds were placed on 1.5% (w/v) agar plates in 4 • C for 1 day. After germination in a greenhouse (20 • C/25 • C and 16 h/8 h light/dark) for 1-2 days, the seedlings were planted in sterilized sand and watered with Fahraeus nitrogen-free nutrient solution (Fahraeus, 1957). S. meliloti 1021 was inoculated after the cotyledons were expanded. Nodules were collected from the taproots of 40 plants at each dpi.
Paraffin sections of nodules were carried out for microscopic observation. Nodules at different dpi were cut longitudinally and fixed with FAA more than 24 h. After dehydrated with gradient ethanol and cleared with dimethylbenzene, the nodules were embedded in paraffin and made into sections. Then the slides were stained with toluidine blue and observed with the Olympus light microscope.
Library Preparation for Long Non-coding RNA-Seq Total RNA was obtained using plant RNA extraction kit RN40 (Aidlab, Beijing, China). Nanodrop2000 (Thermo Fisher, Waltham, MA, United States) was employed to verify RNA concentration and purity. Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States) was used to verify RNA integrity. Ribo-off rRNA Depletion Kit N409-2 (Vazyme, Nanjing, China) was used to remove rRNA. The VAHTS Total RNA-seq Library Prep Kit for Illumina (Vazyme, Nanjing, China) was used for library construction. The libraries were sequenced on an Illumina NovaSeq 6000 platform (PE150).

Identification and Analysis of Long
Non-coding RNA Clean data were produced by removing reads containing adapter and low-quality reads from raw data. HISAT2 v2.0.4 (Kim et al., 2019) was used for sequence alignment. The transcriptome was assembled using StringTie v1.3.1 (Pertea et al., 2015) and Scripture based on the reads mapped to the reference M. truncatula genome MedtrA17_4.0. The assembled transcripts were compared using the Cuffcompare v2.1.1 program (Trapnell et al., 2010). LncRNAs were screened for using the following criteria: (1) Transcripts less than 200 nt were removed. (2) LncRNA transcripts were evaluated for their potential proteincoding with CPC (Kong et al., 2007), CNCI (Sun et al., 2013), Pfam and CPAT  platforms, and the intersection of the four methods were retained.

Identification of Differentially Expressed
Long Non-coding RNAs, Target Gene Prediction and Transposable Element Analysis of Long Non-coding RNAs Differential expression analysis was performed using the DESeq R package v1.10.1 (Anders and Huber, 2010). LncRNAs or mRNAs with p < 0.05 and Fold Change ≥ 1.5 were considered to be differentially expressed. For target gene prediction, Perl scripts were used to search adjacent genes within ± 100 kb of lncRNAs as the cis-target genes, while trans-target gene prediction was based on the complementary sequence using LncTar (Beattie, 2014) prediction program. The target lncRNAs of microRNAs were predicted using TargetFinder (v1.0; Fahlgren and Carrington, 2010). Extensive de-novo TE Annotator (EDTA; Ou et al., 2019) 1 was used for Transposable Element (TE) annotation. The lncRNA overlapping with TE-site but not completely inside a TE was confirmed as TE-associated lncRNA (Wang D. et al., 2017).

Gene Annotation and Functional Analysis of Differentially Expressed Long Non-coding RNAs Target Genes
Gene function was annotated based on the databases of Nr (NCBI nonredundant protein sequences 2 ), Pfam (Protein family 3 ), KOG/COG (Clusters of Orthologous Groups of proteins 4 ), Swiss-Prot (a manually annotated and reviewed protein sequence database 5 ), KEGG (Kyoto Encyclopedia of Genes and Genomes 6 ) and GO (Gene Ontology 7 ). GO enrichment analysis was implemented by the TopGO R packages, and KOBAS (Mao et al., 2005) software was used for KEGG pathway analysis.

Quantitative Real-Time PCR
Total RNA was extracted by TRizol regent and random reverse primer was used for reverse transcription of lncRNA and mRNA. The qPCR was performed on qTOWER 3.0 real-time PCR System (Analytik, Jena, Germany). Primers are listed in Supplementary  Table 1. The relative expression levels of genes were calculated by 2 − CT method.

Statistical Analysis
Statistical analysis was performed using SPSS 19.0 software (IBM, Chicago, IL, United States). Two groups of data were analyzed using the unpaired two-tailed t-test. Significance analysis of length and expression level between TE-and non-TE-lncRNAs was performed by Wilcoxon test.

Nodules at 35 Days Post Inoculation Displayed Senescence
Nodules at 21 and 35 dpi were collected to determine their developmental stage. At 21 dpi, nodules were small and pink, while at 35 dpi, a small proximal section gained a green color, indicating that aging has occurred ( Figure 1A). Paraffin sections stained with toluidine blue were performed to observe the developmental zones. The cells in nitrogen fixation zone of 21 dpi nodules remained healthy with a large number of bacteroids. While at 35 dpi, a small distinct senescence zone was present at the proximal region of the fixation zone ( Figure 1B). In this region, the number of infected cells reduced and the loss of cellular content was observed, indicating the degradation of bacteroids ( Figure 1C). Moreover, some infected cells were abnormal with a very large vacuole. According to the above observations, we determined that 35 dpi nodules have entered the aging stage.
We compared the expression level, transcript and ORF length, exon number, and the isoform number of lncRNAs with mRNAs. The results reflected the different characterizations between lncRNAs and mRNAs. The average expression level of mRNAs was 1.8 times that of lncRNAs ( Figure 2D). Most lncRNAs have a transcript length of less than 1000 nt (71.1%; Figure 2E) and an Frontiers in Plant Science | www.frontiersin.org ORF length ≤ 100aa ( Figure 2G). In contrast, the average length of mRNAs was 2467nt ( Figure 2F) and the ORF of most mRNAs was longer than 100aa ( Figure 2H). The majority of lncRNAs contained less than three exons (Figure 2I), while about 76.5% mRNAs have more than three exons ( Figure 2J). For isoform number, the presence of one or two isoforms is the most common case for both lncRNAs and mRNAs ( Figure 2K).

Identification and Functional Analysis of Long Non-coding RNAs Related to Nodule Senescence
A total of 126 DElncRNAs including 67 up-regulated and 59 down-regulated lncRNAs were identified (N35 vs N21; Figures 3A,B). The distribution of DElnRNAs in chromosomes displayed the same preference with total lncRNAs. Moreover, although chromosome 6 had a small number of DElncRNAs, most of them were up-regulated ( Figure 2C). Many lncRNAs function by regulating gene expression, so the prediction of their target genes can provide insight into their biological roles. A total of 1911 putative cis-regulated and 28 trans-regulated target genes of DElncRNAs were predicted. GO terms analysis of these target genes showed significant differences between the mature and senescent nodules. Notably, for the top 20 terms of cellular component, integral component of membrane was the most significantly enriched term by both cis ( Figure 3C) and trans ( Figure 3D) target genes. KEGG pathway analysis revealed that cis target genes were enriched in RNA polymerase, purine and pyrimidine metabolism, as well as flavonoid and amino acids biosynthesis pathways (Supplementary Figure 1). The trans target genes were enriched in MAPK signaling, plant hormone signal transduction, plant-pathogen interaction and isoflavonoid biosynthesis pathways (Supplementary Figure 1).

Identification of Long Non-coding RNAs Targeting Memebrane-Related Genes and Transcription Factors
Among the target genes of DElncRNAs, 48 genes were identified to be differentially expressed (DEmRNAs) between N35 and N21. TopGO analysis of the DEmRNAs demonstrated that 13 of the 48 DEmRNAs were membrane associated ( Figure 4A). Furthermore, six of the 13 DEmRNAs encoded membrane proteins related to transport function which are two casparian strip membrane proteins, the SNARE protein SYP132, an EamA domain protein, a nitrate transporter NRT1(PTR) and a transmembrane protein (Supplementary Table 4).
Since transcription factors (TFs) play important roles during nodule senescence, we investigated the TF genes targeted by DElncRNAs. Altogether, 41 TF genes belonging to 13 families were identified, of which, MYB (11, 26.8%) and MADS (8, 19.5%) constituted the largest two families containing the high number of target genes ( Figure 4B). We found 7 of the 41 TF genes were differentially expressed between N21 and N35 including three

Identification of Long Non-coding RNAs Targeting Well-Studied Genes Involved in Nodule Senescence
Since some ageing-related genes were reported to play important roles in nodule senescence (de Zelicourt et al., 2012;Berrabah et al., 2014;Pierre et al., 2014;Dhanushkodi et al., 2018;Deng et al., 2019;Trujillo et al., 2019), the lncRNAs probably targeting them (within ± 100 kb of the associated genes) were investigated. We firstly examined the expression of these genes in the data of RNA-seq. As the early molecular markers of nodule senescence, two cysteine proteinase genes, MtCP6 and MtVPE, were significantly upregulated in N35. A NAC family TF gene MtNAC969 which is a negative regulator of nodule senescence also showed upregulated expression. While the expression of the rest genes has no significant difference between N21 and N35. The result was consistent with previous reports. A total of 34 lncRNAs were  predicted targeting to 13 senescence-associated genes ( Table 1). Of these, 9 DElncRNAs including 7 down-regulated and 2 upregulated lncRNAs were identified. As seen from Table 1, most genes could be targeted by multiple lncRNAs. For instance, MtNAC969 was probably targeted by six lncRNAs, two of them showed differential expression. MtCP6 and MtVPE were targeted by one DElncRNAs, respectively. The result implied that lncRNAs played regulatory roles in nodule senescence by targeting some key senescence-related genes.

Prediction of Differentially Expressed Long Non-coding RNAs Targeted by MicroRNAs
Plant lncRNAs can play regulatory roles by acting as the target mimicry of miRNA, so it is necessary to predict the senescence-associate lncRNAs targeting by miRNAs. In total, 49 DElncRNAs were predicted to be targeted by 93 miRNAs (Supplementary Table 5). We found that some miRNAs could target more than one DElncRNAs. As a well-known regulator of multiple physiological processes, miR156 probably target three DElncRNAs. In addition to miR156, some other miRNAs such as miR172 and miR168 were also found to interact with DElncRNAs. Conversely, some DElncRNAs could also bind multiple miRNAs. For example, MSTRG.16162.3 was predicted to bind with eight miRNAs which belonged to four miRNA families. High-throughput sequencing of miRNAs showed that 36 of the above 93 miRNAs were differentially expressed (Supplementary Table 5) including 19 known miRNAs such as miR156, miR172, miR1509 and miR2629, and 17 novel miRNAs.
The percentage of TE-lncRNAs in polyA+ RNA was closed to that in polyA-RNA (63.1% vs 61.6%; Figure 5A). Compared with non-TE-lncRNAs, the average length of TE-lncRNAs was significantly longer and their expression level was relatively lower (Figures 5B,D). Furthermore, for DElncRNAs, the difference in length between TE-and non-TE-lncRNAs was greater ( Figure 5C). We investigated the family of TEs in lncRNAs according to their sequence homology with known TEs. In terms of quantity, the proportion of TE-lncRNAs associated with DNA transposons (2004, 70.3%) was much larger in M. truncatula than that in many plant species (Wang D. et al., 2017;Golicz et al., 2018b). The family which contributed the most to TE-lncRNAs was classified as Helitron, followed by TIR/Mutator, LTR/unknown, LTR/Gypsy and LTR/Copia. Similarly, for the TEs in DElncRNAs, the top three families were also Helitron, Mutator and LTR/unknown, but a small number of LTR/Copia and LTR/Gypsy elements were identified ( Figure 5E).
In addition, we investigated the contribution of different TE families to lncRNA length ( Table 3). We found that the lncRNAs related to LTR/Gypsy accounted for the most, which was different with the result calculated by quantity. Interestingly, the family that contributed most to DElncRNAs was TIR/Mutator rather than LTR/Gypsy. Furthermore, for both total lncRNAs and DElncRNAs, the contribution of DNA transposons to lncRNAs was greater than their contribution to the genome. Especially, DElncRNAs originated from TIR/Mutator accounted for 27.03% which was nearly three times the proportion of TIR/Mutator elements in all the TEs in the genome. However, the proportion of TE-lncRNAs derived from LTR/Gypsy or LTR/Copia was lower than their proportion in the genome.

Quantitative Real-Time Validation
To verify the data of RNA-Seq, eight DElncRNAs were selected randomly for qRT-PCR detection (Figures 6A-H). The expression trends of the DElncRNAs were consistent with the  results of RNA-Seq (Figure 6I), which indicated the reliability of expression analysis. Furthermore, we selected three interesting lncRNAs to check the co-expression tendency of lncRNAs and their putative target genes by qRT-PCR. In 15, 21, 28, 35, and 45 dpi nodules, the expression tendency of lncRNA MSTRG.14267.1 and gene-LOC25492610 (encoding a lysinespecific demethylase) was highly consistent (Figure 6J). While the expression profiles of lncRNA MSTRG28751.1 and gene-LOC25502666 (encoding a bHLH transcription factor) presented an opposite trend ( Figure 6K). Additionally, the expression trends of two genes (senescence-associated gene, newGene_6237 and transmembrane protein gene, newGene_6245) were both consistent with that of MSTRG.31647.1 (Figure 6L).

Long Non-coding RNAs Are Involved in Regulating Nodule Senescence
Nodule senescence leads to the decrease of nitrogen fixation efficiency and affect crop yield. Thus, one effective measure to ensure crop yield is to prolong the period of nitrogen fixation by delaying the onset of nodule senescence. Investigating the regulatory mechanism underlying nodule senescence can provide potential targets for such. However, little research has focused on the lncRNAs related to nodule senescence. In this study, we identified 126 putative nodule senescence-related lncRNAs by strand-specific RNA-seq. Our findings can provide insight into the functions of lncRNAs in nodule senescence and provided candidate targets for nodule senescence regulation.
Transposable Element-Associated Long Non-coding RNAs Play Important Roles in Nodule Senescence TEs are widely distributed in plant genome. Previous work has demonstrated that a number of plant lncRNAs are derived from TEs. Lots of TE-lncRNAs have been identified in Arabidopsis, rice (Wang D. et al., 2017), maize (Lv et al., 2019), cotton (Zhao et al., 2018), tomato (Wang et al., 2016) and soybean (Golicz et al., 2018b). Similarly, our results revealed that more than 60% lncRNAs contained TE sequences, and the proportion was higher in DElncRNAs. In many plant species, TE-lncRNAs mainly originate from retrotransposon especially LTR family (Wang D. et al., 2017). But surprisingly, we found in M. truncatula nodule, Helitron and TIR/Mutator families contributed the most to lncRNAs in quantity and length, respectively, which was distinct with the results in maize (Lv et al., 2019), rice (Wang D. et al., 2017), cotton (Zhao et al., 2018) and soybean (Golicz et al., 2018b), but have some similarities with the finding in Arabidopsis (Wang D. et al., 2017). The result suggested that the contribution of different TE families varied according to plant species and growth conditions. TEs act as the functional domain of lncRNAs (Johnson and Guigo, 2014) and current researches showed that TE-lncRNAs often work as regulators both in plant response to abiotic stress (Wang D. et al., 2017;Lv et al., 2019) and plant development including fruit ripening (Wang et al., 2016) and the control of seedling height (Zhao et al., 2018). However, whether they are involved in nodule senescence remained unknown. An interesting finding in this work is that the contribution of TIR/Mutator to DElncRNAs in length was significantly greater than that of other families, while Helitron contributed the most in quantity. Mutator elements were firstly found in maize and their homologues are distributed in many other plant species, which are called Mutator-like element (MULE). MULEs are able to selectively capture host gene fragments in Arabidopsis (Yu et al., 2000), maize (Talbert and Chandler, 1988) and rice (Juretic et al., 2005). Genes associated with MULEs play important roles in plant growth and development. For instance, in Arabidopsis MULEderived genes acted as the transcriptional regulator of the genes involved in light response (Hudson et al., 2003). The mutation of genes related to MULEs caused delays in plant growth and flowering and reproductive defects (Joly-Lopez et al., 2012). Helitron elements have been reported to change the function and the expression level of genes Liu et al., 2020). However, the contribution of MULE and Helitron to lncRNAs and plant aging is unknown. Our results suggested TE-lncRNAs derived from MULE and Helitron are involved in nodule senescence.

Long Non-coding RNAs Regulate Nodule Senescence by Interacting With miRNAs
In plants, miRNAs are regarded as the control center of diverse biological processes including plant aging (Werner et al., 2021). During plant flowering, the aging pathway was regulated by miR156 and its target, SQUAMOSA PROMOTER BINDING-LIKE (SPL) transcription factors (Gou et al., 2019). SPL genes can increase the expression of miR172, and miR156/SPLs/miR172 constitute the regulatory network of aging pathway (Wei et al., 2017;Werner et al., 2021). Additionally, miR168 was involved in seed senescence in barley (Puchta et al., 2021). In legumes, miR156 and miR172 were essential regulators of nodulation in soybean (Yan et al., 2013;Yun et al., 2022), common bean (Nova-Franco et al., 2015) and alfalfa (Aung et al., 2017). However, the involvement of miRNAs in nodule senescence has not been reported. In our study, miR156, miR172 and miR168 displayed differential expression in N21 and N35, suggesting their roles in nodule senescence. LncRNAs can bind to miRNAs as competing endogenous targets. Here four DElncRNAs were predicted to be targeted by miR156 or miR172, which indicated that DElncRNAs could function in nodule aging by interacting with miRNAs.

Long Non-coding RNAs Regulate Nodule Senescence by Affecting the Material Transport Across Membrane
TopGO analysis of DEmRNAs targeted by DElncRNAs showed that 13 DEmRNAs were associated with membrane component and six of them encoded proteins related to material transport. Of these proteins, SYP132, a Qa-SNARE, was reported to mediate the fusion between vesicles and the target cell membrane. A SNARE in tobacco was responsible for the regulation of cell membrane ion channels (Leyman et al., 1999). SYP132 in A.thaliana mediated the endocytosis of H + ATPase (Xia et al., 2019). Previous work suggested that MtSYP132 localized to the symbiosome membrane and the membrane around the infection threads, indicating its roles in nodulation and nodule development. The symbiosome membrane provides a medium for communication between the bacteroids and host cells and transmembrane ion transport across the symbiotic membrane is crucial for the function and survival of bacteroids (Catalano et al., 2007). Because MtSYP132 was up-regulated during senescence, it was likely to regulate the transport of some special aging-related molecules across the symbiosome membrane. Besides SYP132, two casparian strip membrane proteins and a NRT1(PTR) protein were also up-regulated. Casparian strip membrane proteins mediated the deposition of casparian strip which regulated the transport of water and inorganic salts, and defects in its development led to increased solute leakage Calvo-Polanco et al., 2021). NRT1(PTR) proteins are known as nitrate transporter (Miller et al., 2007), which can also transport other substrates (Waterworth and Bray, 2006;Krouk et al., 2010). A NRT1(PTR) protein in M. truncatula was reported to be essential for lateral root growth and nodule development (Yendrek et al., 2010). In summary, we speculated that lncRNAs played a role in nodule senescence by affecting the material transport across membrane.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the SRA: https: //www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA810777, accession number PRJNA810777.

AUTHOR CONTRIBUTIONS
YL, LZ, and LY conceived the project and design the protocol. XQ, JY, and LY performed the experiments. LY, TH, TW, ZL, and YL performed the data analysis. YL and LZ wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the National Natural Science Foundation of China (31760062 and 32060062) and Guangxi Natural Science Foundation (2018GXNSFAA138118).

ACKNOWLEDGMENTS
We are very grateful to Professor Youguo Li who provided Sinorhizobium meliloti 1021 and plant seeds.

SUPPLEMENTARY MATERIAL
The