To profile the LncRNA transcripts in response to salt stress, we performed the time-series dynamic analysis of RNA-sequencing on quinoa roots exposed to high salinity conditions (Supplementary Table 1). In total, 464 million raw reads were generated, and 430 million clean reads were obtained after cleaning. Approximately 85.30% (367 million) of the clean reads were mapped to the quinoa reference genome and assembled into transcripts. Cleaned data (∼30 Gb) was first mapped to the reference genome, and then alignments were assembled into transcripts in each sample (Supplementary Table 2). Transcripts from all samples were merged to form a unified set of transcripts, which consisted of 188,663 genes and 234,387 transcripts. There were 146,440 transcripts (62.5% of the predicted transcripts) with only one exon. There were 1.24 transcripts per gene and 3.15 exons in a transcript.

Coding potential ability was analyzed by CPC2, from which there were 153,751 non-coding LncRNAs identified (Supplementary Table 3). A genome location study showed these LncRNAs were intensely distributed in 18 chromosomes of C. quinoa, and the distribution intensity in the two ends of each chromosome was higher than in other parts (Figure 1A). According to the relationship between a transcript and the closest reference transcript, LncRNAs were grouped into different classes represented by characters and symbols by GffCompare software (Figure 1B). “u” was the most abundant type (72.8%), showing that most LncRNAs came from the intergenic region. The second most abundant class was the “x” type (12.1%), followed by “i” (6.7%) and “p”(3.3%). These results indicated that LncRNAs seldom overlapped with reference transcripts. As shown in Figures 1C,D, the length of LncRNAs was generally much shorter than the coding RNAs, and LncRNA had less exon contained than the coding ones. The difference in expression pattern between coding RNAs and LncRNAs was measured statistically using a two-tailed Mann-Whitney U-Test (Figure 1E). The expression pattern of LncRNAs was significantly different from that of the coding genes. The overall expression level of LncRNA was significantly lower than that of coding RNAs at 0, 0.5, and 2 h of -salinity treatment (p < 2.2e-16). However, at 24 h of treatment, the overall expression level of LncRNAs increased sharply and was more significant than that of coding RNAs. The expression level of coding RNAs was relatively stable within 24 h of treatment. For further validating the reliability of LncRNAs, 50 randomly selected LncRNAs were aligned against Pacbio full-length cDNA datasets from Jarvis et al. (2017). A 47 of 50 LncRNAs were successfully aligned into full-length transcripts. Among them, 46 full-length cDNAs corresponding to LncRNAs were also no coding potential by CPC2 program prediction. These cross-validation results revealed that a set of high-confidence LncRNAs was obtained in the study (Supplementary Table 8).

Differentially expressed transcripts, including coding RNAs and LncRNAs, were then identified by DESeq2 (Supplementary Table 4). We identified 4,460 DE-LncRNAs, of which 214, 1,731, and 3,102 were identified at 0.5, 2, and 24 h of treatment (Figure 2A). Only 54 LncRNAs were differentially expressed at all the three-time points, occupying 1.2% of total DE-LncRNAs (4,460) (Figure 2B). Totally 6,791 DE-coding RNAs were also identified, 104 (1.5%) differentially expressed at three-time points (Figure 2B). At 0.5 h of treatment 75% (161) DE-LncRNAs and 77% (237) DE-coding RNAs were upregulated, while at 2 h of treatment 48% (831) DE-LncRNAs and 56% (2,197) DE-coding RNAs were upregulated. At 24 h of treatment, as high as 81% (2,513) DE-LncRNAs were upregulated, whereas only 39% (1,462) DE-coding RNAs were upregulated (Figure 2C and Supplementary Table 2).

Previous studies have suggested lncRNAs may play cis regulation role against neighboring genes. To investigate this possibility, we further measured the expression correlation between salt-responsive LncRNAs and their closest neighboring gene in either the 5′ or 3′ direction, yielding a dataset of 1,740 LncRNA-Coding Genes pairs (differentially expressed LncRNAs and their neighboring genes). The correlation analysis showed lncRNAs were strongly and highly significantly correlated with the expression of their closest neighboring gene (r = 0.346, p-value < 2.2e-16) (Figure 2D). Suggesting that lncRNAs may either be involved in cis-acting regulation or are subject to some of the same cis-acting regulatory features as their neighboring genes.

To infer the potential biological functions of the LncRNAs, a weighted gene co-expression network consisting of both LncRNAs and coding RNAs based on expression profiles was constructed by the WGCNA program (Langfelder and Horvath, 2008). The soft-thresholding power was predicted and defined as 7, as it was the lowest power for which the scale-free topology fit index reached 0.90 (Figures 3A,B). There were finally 36 modules (Figure 3C) generated, and they were named from M1 to M36. The relationship between modules and salinity treatment was calculated, of which seven modules were highly relevant to salinity treatment (Figure 3D). M30 (r = 0.93, p = 1.57e-05) and M12 (r = 0.99, p = 9.12e-10) were upregulated significantly at 0.5 and 2 h, respectively; M17 (r = 0.90, p = 5.48e-05) and M32 (r = −0.99, p = 4.71e-09) modules upregulated and downregulated at 24 h, respectively; M13 (r = 0.90, p = 5.53e-05) upregulated at 0.5 and 2 h both; M11 (r = 0.90, p = 6.35e-05) upregulated at both 2 and 24 h, while M33 (r = 0.90, p = 6.35e-05) downregulated at the same time. No module was significantly regulated at all the three-time points. The percentage of LncRNAs in salinity-responsive modules ranged from 20 to 40%. The genes with the highest connectivity in each module were selected as hub genes. Totally, 210 hub genes were identified from these seven salinity-responsive modules listed above, which constituted a subnetwork (Figure 3E and Supplementary Table 6). The hub genes of M30 included both TF genes and LncRNAs, and so did M13, M12, and M11 modules. This implied that LncRNAs within these modules might interact with transcript factors and their role in salinity response. However, TF genes were not included in the hub genes of M17 and M32 modules. Instead, more than half of the hub genes were LncRNAs in them. In M17 module as high as 23 LncRNAs were at the hub position (Supplementary Table 6).

Gene ontology analysis of coding RNAs was performed to predict the function of the LncRNAs within the same module (Supplementary Table 5). The most representative GO terms of high-salinity responsive modules are shown in Figure 4A. A large part of the GO terms was enriched in all the seven high-salinity responsive modules. They were related to various critical biological processes, including biological regulation (GO:0065007), regulation of gene expression (GO:0010468), response to stress (GO:0006950), regulation of metabolic process (GO:0019222), transcription (GO:0006350) and developmental process (GO:0032502). Transport (GO:0006810) and metabolic process (GO:0008152) were enriched in six salinity-responsive modules. A small part of the GO terms was enriched in only two or three modules such as cell cycle (GO:0007049). Within these modules, there were several genes directly responding to salt stress. For instance, in M11, a transcript encoding a salt tolerance zinc finger (C2H2 type) which is highly homologous to the Arabidopsis gene (AT1G27730.1); in M12, there is a transcript encoding calcineurin B-like protein 1, which is highly homologous to Arabidopsis gene (AT4G17615.1), it might function as a positive regulator of salt and drought responses and as a negative regulator of cold response, and mediates the activation of AKT1 by CIPK proteins (CIPK6, CIPK16, and CIPK23) in response to low potassium conditions; In M13, one transcript homologous to Arabidopsis gene AT5G12010 was included, which might respond to salt stress and function in ABA-activated signaling pathway.

To validate the reliability of salt-responsive LncRNAs, 15 salt-responsive LncRNAs were randomly selected and then subjected the 0 and 2 h salt-treated samples to quantitative real-time PCR (qRT-PCR) to compare expression changes between replicated control and salt-treated. As a result, 11 of 15 salt-responsive LncRNAs were successfully detected and showed a high degree of consistency (r = 0.877, p = 0.000392) between RNA-seq and qRT-PCR (Figure 4B and Supplementary Table 7).

Chenopodium quinoa cultivar QQ056 were acquired from the USDA-ARS National Plant Germplasm System (NPGS) with permission for scientific research. Detailed information on the variety could be found at npgsweb.ars-grin.gov (accession: PI 584524). The experiments were performed in a phytotron at 28°C under 16 h/8 h (day/night) photoperiod. Seeds of quinoa were surface sterilized, germinated for 7 days, and then moved into a half-strength Hoagland solution. Seedlings of 28-day-old were transferred to half-strength MS with 300 mM NaCl for salt treatment. Roots of quinoa were gathered at 0, 0.5, 2, and 24 h of treatment and stored at –80°C for further investigation. Three biological repeats were used at each time point.

Total RNA was isolated, purified and concentrated with an RNAprep Pure Plant Kit (Tiangen, China). A Thermo 2000 Bioanalyzer evaluated the concentration and quality of RNA with an RNA NanoDrop (Thermo Scientific, United States). cDNA libraries were conducted using TruSeq mRNA Sample Prep Kit (Illumina, United States). RNA sequencing (RNA-seq) was performed on the Illumina NovaSeq 6000 platform (Illumina, United States) at Annoroad Gene Technology (Beijing) Co., Ltd., China.

Raw data were trimmed with low-quality bases and short reads (< 50 bp) using Fastp v0.23.2 (Chen et al., 2018). Cleaned reads were then mapped into the quinoa reference genome (Jarvis et al., 2017) by splice-aware alignment method using STAR v2.7.10a with two-pass mode (Dobin et al., 2013). The mapped reads of each sample were separately assembled by the reference annotation-based transcript (RABT) assembly algorithm and then combined with known transcript annotation (Jarvis et al., 2017) into an updated GTF file using StringTie v2.2.1 (Pertea et al., 2016). Finally, the abundance of transcripts was quantified and normalized with HTseq (Putri et al., 2022) and DESeq2 (Love et al., 2014). The transcripts with fragments per kilobase of transcript per million fragments mapped (FPKM) < 1 in more than three samples were excepted for downstream analysis.

Computational prediction of quinoa LncRNAs was followed as described by Lv et al. (2019) with some custom modification. Sequences of transcripts were firstly retrieved using Gffread v0.12.2 (Pertea and Pertea, 2020). Then transcripts were evaluated in their coding potential using CPC v2.0 (Kang et al., 2017). Default parameters were used. Finally, non-coding transcripts longer than 200 bp were considered as LncRNAs. The updated GTF file was also compared with the reference GTF file using gffCompare (Pertea et al., 2016) to generate class codes representing the position information between the updated transcript and the closest reference transcript.

For differential expression analysis, we compared every time point (0.5, 2, and 24 h of treatment) with CK (0 h) using DESeq2 (Trapnell et al., 2010) based on the negative binomial distribution. Differentially expressed should be fulfilled the following criteria: (I) fold change > = 1; (II) adjusted p value < 0.05.

A weighted co-expression network was further constructed for linking coding RNAs and LncRNAs using the WGCNA program (Zhang and Horvath, 2005; Langfelder and Horvath, 2008). The correlation between module eigengenes and salt treatment was then calculated. Modules that showed a significant correlation (| r| > 0.9, p < 0.001) with a specific time of treatment were recognized as salt-responsive modules (Xue et al., 2013). The hub genes of each module were worked out based on the Topological overlap matrix (TOM).

Gene Ontology (GO) enrichment analysis was performed with agriGO v2.0 toolkit (Tian et al., 2017). Fisher’s exact test was applied for the enrichment analysis, and the false discovery rate (FDR) was assessed using the Yekutieli method. GO terms with an FDR less than 0.05 were considered significant.

To validate the reliability and accuracy of the LncRNA analysis, ten salt-responsive LncRNAs were chosen and measured by quantitative real-time PCR (qRT-PCR). The first-strand cDNA of 0 and 2 h samples were synthesized using the PrimeScript TM RT-PCR Kit (Takara, Japan). We conducted the qRT-PCR on an ABI StepOnePlus Real-Time PCR System (Applied Biosystems, CA, United States) with FastStart Universal SYBR Green Master (Roche, Germany) according to the manufacturer’s protocol. The primers were designed using Primer Premier v.5.0 software (Premier Biosoft International, CA, United States). CqEF1a (AUR62020767) was used as an internal standard to normalize the relative expression level and determine expression values based on the 2^–ΔΔCt method. The primers for qRT-PCR were presented in Supplementary Table 7.