Genome sequence and salinity adaptation of the desert shrub Nitraria sibirica (Nitrariaceae, Sapindales)

Abstract

The genetic bases of halophytes for salinity tolerance are crucial for genetically breeding salt-tolerant crops. All natural Nitrariaceae species that exclusively occur in arid environments are highly tolerant to salt stress, but the underlying genomic bases to this adaptation remain unknown. Here we present a high-quality, chromosome-level genome sequence of Nitraria sibirica, with an assembled size of 456.66 Mb and 23,365 annotated genes. Phylogenomic analyses confirmed N. sibirica as the sister to all other sampled representatives from other families in Sapindales, and no lineage-specific whole-genome duplication was found except the gamma triplication event. Still, we found that the genes involved in K⁺ retention, energy supply, and Fe absorption expanded greatly in N. sibirica. Deep transcriptome analyses showed that leaf photosynthesis and cuticular wax formation in roots were enhanced under salt treatments. Furthermore, many transcription factors involved in salt tolerance changed their expressions significantly and displayed tissue- and concentration-dependent signalling in response to salt stress. Additionally, we found vacuolar Na⁺ compartmentalization is an ongoing process under salt treatment, while Na⁺ exclusion tends to function at high salt concentrations. These genomic and transcriptomic changes conferred salt tolerance in N. sibirica and pave the way for the future breeding of salt-tolerant crops.

1. Introduction

Soil salinization is a global problem that threatens crop growth and production and hampers modern agriculture’s sustainable development.^1–3 High salt concentrations induce ionic imbalance and hyperosmotic stress that may lead to poor plant growth or death.^4,5 However, many plants occurring in salty or arid environments have developed a series of genetic changes to prevent salt toxicity and maintain osmotic homeostasis.^6–13 For example, in the desert poplar, the up-regulated expression of ion transport-related genes and the expansion of genes involved in maintaining intracellular homeostasis make this species salt-tolerant.⁷ In grasses, salt tolerance is enhanced by some duplicated genes resulting from whole-genome duplication (WGD).¹² However, the genomic changes of many other salt-tolerant lineages remain largely unexplored.

The family Nitrariaceae in the broad sense (including Tetradiclidaceae and Peganaceae), in the order Sapindales, comprise around four genera and 20 species that mainly occur in arid and salty environments and show high salt tolerance.^14,15 However, other families of this order, including Rutaceae, Anacardiaceae, and Aceraceae are distributed in subtropical and tropical environments, and most species are sensitive to salt stress. Therefore, Nitrariaceae comprise a specific lineage to examine genomic changes for salt tolerance compared with other Sapindales families. Furthermore, the family has attracted research interests in many other aspects, such as secondary chemistry, biogeography, and floral morphology.^16–19 For example, Nitrariaceae plants are rich and diverse in metabolites, of which numerous chemical constituents isolated from the genera Nitraria and Peganum are shown to have anticancer and antioxidant effects^18,19; In addition, although a small family but it has multiple intercontinental discontinuous distribution patterns, making it an ideal plant group to investigate biogeographical pattern in drylands.¹⁷ However, up to now, no genome sequence is available for Nitrariaceae, although de novo genomes have been reported for other species across several families in Sapindales.^9,20,21 In this study, we aimed to assemble the genome sequence of the type species of the family, Nitraria sibirica Pall. This shrub, occurring in the desert and salt-alkaline regions like other species, can tolerate extremely saline and drought conditions (Fig. 1A).^22–24 Both physiological and biochemical tests have also confirmed such traits.^24,25 For example, mesophyll vacuoles of leaves act as the main Na⁺ sink to store excess sodium ions,²⁶ while the roots have superior K⁺ retention ability.^26,27

Here, we report a high-quality chromosome-scale genome of N. sibirica from a combination of several sequencing platforms and strategies, including Illumina short reads sequencing, Oxford Nanopore Technologies (ONT) long reads sequencing, and chromosome conformation capture (Hi-C) technology. RNA-seq-based transcriptomics was used to optimize gene structural prediction and annotation and explore differentially expressed genes (DEGs), co-responsive expression patterns, and their regulatory networks under different salt concentrations across various tissues. As the first genome sequence of the family Nitrariaceae, we aimed to use it to address the following questions. (1) Did special genomic changes, for example, polyploidization that usually increases the salinity adaptation,^12,28,29 occur in this lineage? (2) If not, did gene families expand their copies as found in the other salt-tolerant species?⁷ Especially what roles do transcription factors (TFs) play in response to salt stress? Moreover, (3) do TFs and their associated regulatory networks display tissue- or concentration-dependent signalling in response to salt stress?

2. Materials and methods

2.1. Plant materials

Plant samples for this study were collected from salinized land (E: 103°34ʹ17.9472″, N: 39°0ʹ30.5172″) near Minqin Desert Botanical Garden, Gansu Province, China. The collected tissues of roots, stems, leaves, and flowers of N. sibirica were immediately stored in liquid nitrogen and subsequently used for genome or transcriptome sequencing. Plant seeds collected from the same population were germinated and planted in the greenhouse for subsequent salt stress treatment and transcriptome sequencing.

2.2. Sequencing and genome assembly

Total DNA was extracted from the fresh young leaves using the sodium dodecyl sulphate method and purified with the QIAGEN^® Genomic kit (Hilden, Germany) according to the standard operating procedure. The DNA library for Illumina short reads with 300 bp insert sizes were constructed using the Illumina DNA Prep Kit and sequenced by Illumina HiSeq 4000 platform in 150-bp paired-end mode. The genomic library of long-reads was constructed and sequenced utilizing PromethION, an Oxford Nanopore Technology sequencer, to obtain raw sequencing data with 20 Kb insert sizes. Then a common standard was used to filter raw reads (removal of sequencing adapter and trimming of low-quality bases with an ‘average quality score < 7’). The GC content distribution range of filtered long reads is shown in Supplementary Fig. 1. The Hi-C sequencing was carried out as follows. Interoperable DNA fragments were obtained by cross-linkage of 2% formaldehyde with sampled DNA, and Dpn II restriction enzymes were used to digest chromatin, then the Illumina Nova-seq platform was used for sequencing in 150-bp paired-end mode.³⁰

The genome size of N. sibirica was estimated with a k-mer depth-frequency distribution (k-mer = 17) using KmerFreq_AR as implemented in SOAPdenovo2³¹ from Illumina 150 bp paired-end read data. For de novo genome assembly, the NextDenovo v2.3.0 software (https://github.com/Nextomics/NextDenovo), an overlap layout-consensus (OLC) assembled method specially developed for ONT reads, was used to generate the primary assembly for N. sibirica as follows. First, the program NextCorrect was used to self-correct the original subreads, and then the consensus sequences (CS reads) were obtained. Second, CS reads were compared with the NextGraph module to capture each correlation between them. The preliminary genome was assembled based on the correlation relationships of CS reads. Finally, to improve the accuracy of the assembly, the contigs were refined with Racon v1.3.1³² using ONT long reads and Nextpolish v1.2.4³³ using Illumina short reads, respectively. After, we assessed the potential contamination by comparing the assembled contigs with the National Center for Biotechnology Information (NCBI)—non-redundant (NR) nucleotide database (https://www.ncbi.nlm.nih.gov/nucleotide/) using BLAST with an E-value cutoff of 1e – 5. As shown in Supplementary Table 1, no external contaminations (endophytes and bacteria, etc.) were found and four contigs derived from mitochondrion and chloroplast was removed from subsequent chromosome-level assembly. For the Hi-C-assisted assembly, 370 million paired-end reads were generated. The raw Hi-C data were then quality controlled using Hi-C-Pro v2.8.1,³⁴ and valid interaction paired reads were used for subsequent scaffolding. Subsequently, the assembled scaffolds were divided into bins of 100 Kb size, and the number of contacted Hi-C read pairs between each pair of bins was calculated. Based on the contact information, the scaffolds were further clustered, sorted, and oriented to chromosomes by LACHESIS.³⁵ Finally, some errors in position and orientation were manually adjusted based on chromatin interaction frequency by using Juicebox (https://github.com/aidenlab/Juicebox). The completeness of genome assembly was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO) v4.0.5³⁶ and Core Eukaryotic Gene Mapping Approach (CEGMA) v2.³⁷ We evaluated the accuracy of the assembly by mapping all the Illumina paired-end reads to the assembled genome using Burrows–Wheeler Aligner (BWA) v0.7.12.³⁸ Mapping rate and genome coverage were assessed using SAMtools v0.1.1855.³⁹ Besides, base accuracy of the assembly was calculated with BCFtools v1.8.0.⁴⁰

2.3. Genome annotation

Tandem repeats were annotated using GMATA v2.2⁴¹ and Tandem Repeats Finder (TRF) v4.07b.⁴² Transposable elements (TEs) were identified using a combination of ab inito and homology-based methods. Briefly, an ab inito repeat library was first predicted using MITE-hunter⁴³ and RepeatModeler v1.0.11⁴⁴ with default parameters. For further identification of the repeats throughout the genome, RepeatMasker was used to search for known and novel TEs by mapping sequences against the de novo repeat library and Repbase.⁴⁵ Finally, overlapping transposable elements of the same repeat class were collated and combined. In addition, we used tRNAscan-SE v2.0,⁴⁶ Infernal v1.1.2,⁴⁷ and RNAmmer v1.2⁴⁸ to predict non-coding RNAs (including miRNAs, tRNAs, rRNAs, and snRNAs).

Three independent approaches were used for protein-coding gene prediction in a repeat-masked genome, including ab initio prediction, homology search, and RNA-seq-based. In detail, GeMoMa v1.6.1⁴⁹ was used to align the homologous peptides from four closely related species with high-quality published genomes (Acer yangbiense, Citrus maxima, Athaliana thaliana, and Vitis vinifera) to the repeat-masked assembly of N. sibirica to obtain the gene structure information. Meanwhile, Augustus v3.3.1⁵⁰ was run with default parameters for ab initio gene prediction with the training set. For RNA-seq-based gene prediction, filtered mRNA-seq reads were aligned to the N. sibirica reference genome using STAR v2.7.3a⁵¹ with default parameters. The transcripts were then assembled using Stringtie v1.3.4d,⁵² and open reading frames (ORFs) were predicted using PASA v2.3.3.⁵³ Finally, comprehensive gene sets obtained by different strategies were integrated by EVidenceModeler (EVM) v1.1.1.⁵³ For gene function annotation, BLASTP (E-value < 1 – e05) was used in two public databases, including SwissProt (https://web.expasy.org/docs/swiss-prot_guideline.html) and the NCBI-NR protein database (https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/) to obtain functional annotations of predicted protein-coding genes. Protein structural domains were predicted by InterProScan v5.23.⁵⁴ Gene Ontology (GO) entries for each gene were searched using Blast2GO v3.0,⁵⁵ and pathway information was assigned using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/).

2.4. Phylogenomic analyses

To explore the phylogenetic position of Nitrariaceae, we collected the protein-coding gene sequences of 11 other sequenced plant species in eudicots: including four Sapindales species (A. yangbiense, Citrus clementina, Poncirus trifoliata, and Pistacia vera) and seven additional eudicot species (A. thaliana, Corchorus capsularis, Caragana korshinskii, Cucumis sativus, Populus trichocarpa, Tripterygium wilfordii, and V. vinifera). The detailed genome information and sources are listed in Supplementary Table 2. Only the longest transcript of each gene was used for subsequent analysis. We used OrthoFinder v2.4.0⁵⁶ to cluster the protein sequences of N. sibirica and the other 11 species to putative gene families (or ortho groups) based on sequence similarity with default parameters. In total, 554 strict single-copy gene families were obtained. MAFFT v7.471⁵⁷ was used to perform multiple sequence alignments of the protein sequences in each sing-copy gene family. Then, the protein sequence alignments were converted into corresponding DNA sequence alignments using PAL2NAL v14.⁵⁸ We used maximum likelihood and coalescent-based methods for phylogenetic analyses. For maximum likelihood tree reconstruction, the DNA sequence alignments of single-copy gene families were concatenated into a super-gene matrix, and then model selection and tree inference were automatically performed with IQTREE v2.1.2.⁵⁹ For coalescent-based phylogenetic analysis, the trees of each gene family were constructed using IQTREE and subsequently summarized using ASTRALL-III v5.7.5.⁶⁰ In addition, the chloroplast genome of N. sibirica was assembled and annotated with GetOrganelle v1.7.5⁶¹ and CPGAVAS2,⁶² respectively. The chloroplast genomes of other representative species were downloaded from the public database (Supplementary Table 3). Phylogenetic analysis of chloroplast genomes was performed using the maximum likelihood method as in the nuclear data. The divergence time estimation of the 12 species was done using the MCMCtree pipeline of the PAML program v4.9.⁶³ We used three calibration points obtained from the TimeTree database (http://timetree.org/), including the crown group of Sapindales (56–99 Mya), the split of C. capsularis–A. thaliana (84–95 Mya), and the crown group of V. vinifera and other species (106–125 Mya).

2.5. Comparative genomic analyses

To estimate the evolutionary dynamics of gene families, we employed the computational analysis of gene family evolution (CAFÉ) v4.2 program⁶⁴ to explore gene family expansion and contraction events across the time-calibrated phylogeny based on a random birth and death model. GO enrichment analysis of the expanded and extracted gene families was performed with TBtools.⁶⁵

We used Ks-based age distributions of paralogous genes and intra- and inter-genomic synteny to analyze the WGD history in N. sibirica. Three other Sapindales species (A. yangbiense, C. clementina, and P. trifoliata) with published chromosomal-level genomes were used for comparisons, and V. vinifera, a lineage without any further WGDs after the gamma event, was used an outgroup. Ks values of paralogous genes were calculated using the Nei–Gojobori method⁶⁶ by the yn00 package of PAML v4.9.⁶¹ Intra- and inter-genomic synteny comparisons were performed using WGDI with default parameters.⁶⁷

2.6. Identification of key gene members relevant to salinity adaptation

We used BLASTP and hidden Markov model (HMM) search⁶⁸ to identify the gene members in ion transport, tricarboxylic acid cycle, nicotianmine biosynthesis, and salt overly sensitive pathways. Protein sequences of gene members of each gene family in A. thaliana were retrieved from the Arabidopsis Information Resource (TAIR) (https://www.arabidopsis.org) or relevant references. The retrieved A. thaliana proteins were then searched against the protein sequences in other four Sapindales species (A. yangbiense, C. clementina, P. trifoliata, and P. vera) and V. vinifera using BLASTP with E-value < 1e – 10 and sequence identity > 30%. Meanwhile, the HMM profiles with protein domains of each family were downloaded from the Pfam website (http://pfam.xfam.org), and then batch searched in each species with HMMER v3.3.2.⁶⁸ All candidates of each gene family were used for maximum likelihood tree construction using IQTREE⁵⁹ as previously described. Any outliers in the phylogeny were manually checked and removed, and the remaining genes were regarded as true members of a gene family.

2.7. Transcriptome sequencing and analyses

Two months N. sibirica seedlings of equal growth conditions were used for salt treatments. It has been suggested that the optimal salt concentration for N. sibirica growth is around 200 mM NaCl, and it could tolerant at least a 400 mM NaCl salt treatment.^22,24,26 Therefore, we set two salt concentrations: moderate with 200 mM NaCl treatment and high with 400 mM NaCl treatment. Three major tissues (roots, stems, and leaves) were sampled to examine the gene expression profiles in response to salt stress. Salt treatments were performed by adding 50 mM NaCl every half day to avoid overstimulation until it reached the two treatment concentrations (T200 and T400) and a further week of stress. Seedlings without salt treatment served as control. A total of three biological replicates were set up. Roots, stems, and leaves of NaCl-treated and control seedlings were collected and sorted in liquid nitrogen for subsequent sequencing. Total RNA extraction, library construction, and sequencing were performed by Novogene Co, Ltd (Beijing, China) on the Illumina Nova-seq platform.

A total of 27 transcriptomes from three tissues with two salt concentrations and controls were generated and yielded 193.38 Gb (Supplementary Table 4). For each sample, the quality of raw reads was evaluated with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and low-quality reads was trimmed with fastp v0.20.0.⁶⁹ After quality control, clean reads were mapped to the N. sibirica genome using HISAT2 v 2.1.0.⁷⁰ StringTie v1.3.4d⁵³ was used to calculate the transcripts per million (TPM) of genes in each sample. Differentially expressed genes between treatments and controls in different tissues were identified using the DESeq2 v3.11.⁷¹ The candidate genes with at least 2-fold differential expression levels in comparisons and a false discovery rate (FDR) value < 0.05 were identified as DEGs. To check whether differentially expressed TFs evolved under positive selection, we analyzed the accumulation of non-synonymous changes by comparing N. sibirica with three salt-sensitive species (A. yangbiense, C. clementina, and P. trifoliata) in Sapindales following the procedures of Ma et al. (2013).⁷ The R package WGCNA v1.69⁷² was used to perform weighted gene co-expression network analysis and the resulting networks were visualized using Cytoscape v3.9.1.⁷³ We analyzed the average expression level (TPM) tendency of genes in each tissue (root, stem, and leaf) under different salt concentrations (control, T200, and T400) on the Omicshare Tools website (https://www.omicshare.com/). The KOBAS website (http://kobas.cbi.pku.edu.cn/kobas3/) was used to perform a KEGG enrichment analysis of the genes in each cluster.

3. Results

3.1 Chromosomal-level genome assembly of N. sibirica

Nitraria sibirica is a diploid halophyte (2n = 2x = 24),⁷⁴ with an estimated genome size of 482.77 Mb based on genome survey analysis (Supplementary Table 5; Supplementary Fig. 2). For de novo whole-genome sequencing, we produced 38 Gb (~79 × coverage) Nanopore long-read sequencing data, 32 Gb (~67 × coverage) short reads by Illumina sequencing, and 55 Gb (~116 × coverage) of Hi-C paired-end reads (Supplementary Tables 6–8). After self-correction of long reads, contig-level primary assembly, polishing with short reads, removing contaminations, and further scaffolding (see Materials and Methods), we obtained a correct genome assembly totalling 456.66 Mb. The assembled genome is slightly smaller than the estimated genome, consisting of 170 contigs with an N50 value of 9.07 Mb and a longest contig of 21.14 Mb (Table 1, Supplementary Table 9). We further used Hi-C contact information to construct a chromosomal-level genome assembly that resulted in 98.44% of the genomic sequences assigned unambiguously to discrete chromosome locations (Fig. 1B and Supplementary Fig. 3), including 12 pseudochromosomes with sizes ranging from 24.3 to 56.2 Mb (Supplementary Table 10).

To assess the accuracy of the final genome assembly, the genomic data from Illumina short reads and transcriptome reads (i.e. those from roots, stems, leaves, flowers, fruits, and various tissues under salt treatments) were mapped to the assembled genome, resulting in mapping rates of above 99.5% and 92.3%, respectively (Supplementary Table 4). Furthermore, 98.14% of the 1,614 BUSCO genes were detected in the assembled genome (Supplementary Table 11). Similarly, 97.98% of the ultra-conserved core eukaryotic genes, based on CEGMA analysis, were completely detected in the assembly (Supplementary Table 12). The genome assembly achieved ‘reference quality’ with an estimated long-terminal repeat (LTR) assembly index (LAI) value of 12.17 (Supplementary Table 13).⁷⁵ Together, these results suggest that the continuity and integrity of the N. sibirica genome assembly are both of high quality.

3.2 Genome annotation of N. sibirica

We annotated 23,365 protein-coding gene models (Table 1) from the N. sibirica genome by integrating de novo, homology, and RNA-seq-based prediction methods (see Materials and Methods). Of these, 19,918 high-confidence genes (85.25% of the total genes) were supported by transcriptome gene expressive evidence (Supplementary Table 14). The average gene length, coding-sequence length, and exons per gene are 4,154.54 bp, 1,274.13 bp, and 5.91, respectively (Table 1), like in other eudicot species (Supplementary Fig. 4). About 87.17% of the total genes were annotated as functional by at least one of the following protein-related databases: GO database (54.71%), KEGG database (32.6%), NCBI-NR protein database (86.48%), the Swiss-Prot database (73.4%), and Clusters of Orthologous Groups for Eukaryotic complete Genomes (KOG) database (49.29%) (Supplementary Table 15). BUSCO assessment showed that our annotation recovered 97.71% of the gene set (Supplementary Table 16).

The N. sibirica genome is composed of 62.85% of repetitive elements (Table 1), which is slightly lower than in A. yangbiense (68.0%)²¹ and P. vera (70.7%),⁹ but much higher than in C. clementina (45%).²⁰ With retrotransposons (Class I) and transposons (Class II) transposable elements (TEs) representing 45.61% and 12.79% of the genome, respectively (Supplementary Table 17). LTR retrotransposons are the most abundant TEs, representing 40.64% of the assembled genome (Supplementary Table 17). Moreover, most of these LTRs burst around one million years ago (Mya) (Supplementary Fig. 5). In addition to protein-coding genes and repetitive elements, we predicted 734 transfer RNAs (tRNA), 264 microRNAs (miRNA), 865 small nuclear RNAs (snRNA) and 3,789 ribosomal RNAs (rRNA; Supplementary Table 18).

3.3. Phylogenomic inference of the evolutionary position of Nitrariaceae

To resolve the phylogenetic position of Nitrariaceae, we extracted 554 strict single-copy nuclear genes from the N. sibirica genome and other representative species with published genomes, including four of Sapindales—C. clementina (Rutaceae), P. trifoliata (Rutaceae), P. vera (Anacardiaceae), and A. yangbiense (Aceraceae), six of other eudicots lineages—C. capsularis (Malvales), A. thaliana (Brassicales), P. trichocarpa (Malpighiales), T. wilfordii (Celastrales), C. korshinskii (Fabales), and C. sativus (Cucurbitales), and one outgroup species—V. vinifera (Vitales; Supplementary Table 2). We employed maximum likelihood and coalescent-based methods to reconstruct the phylogeny of Nitrariaceae and other representative lineages in eudicots (see Materials and Methods). The two methods retrieved the same topology, that is, N. sibirica of Nitrariaceae as a sister to all other representatives of the remaining families in Sapindales (bootstrap support (BS) = 100%; Fig. 2 and Supplementary Fig. 6A). In addition, we also assembled a complete chloroplast genome of N. sibirica and performed plastid phylogenomic analysis for N. sibirica and related species. The results showed that the topology was consistent with the nuclear phylogeny for the placement of N. sibirica, indicating the robustness of Nitrariaceae as a member of Sapindales (BS = 100%; Supplementary Fig. 6B). Based on this well-resolved phylogeny, the divergence time between N. sibirica and other Sapindales species was estimated at around 75 Mya (Supplementary Fig. 7).

3.4. Comparative genomic analyses

The distributions of synonymous substitutions per site (Ks) and intragenomic syntenic analysis suggested that the genomes of Sapindales species have not undergone any additional genome duplication, apart from the common whole-genome triplication (gamma) event shared by all core-eudicots species (Fig. 2A and Supplementary Fig. 8).⁷⁶ Intergenomic syntenic analysis between N. sibirica and three Sapindales species (A. yangbiense, C. clementina, and P. trifoliata) and V. vinifera revealed an apparent 1:1 syntenic depth ratio, which further confirmed that no lineage-specific WGD occurred in N. sibirica and also other species of the order (Fig. 2B).

Since no WGD exists, we examined the expansion and contraction of gene families in N. sibirica. We found that many gene families in the N. sibirica genome have contracted (i.e. 4,481) and while only a few gene families have expanded (i.e. 1,073) (Fig. 2C; Supplementary Table 19). The contracted gene families are mainly involved in essential molecular functions (i.e. binding) and biological processes (i.e. nucleic acid metabolism) (Supplementary Fig. 9). However, GO enrichment indicated that all of the expanded gene families are involved in essential biological processes related to salinity adaptation, including potassium ion transport (GO: 0006813), and potassium ion transmembrane transport (GO:0071805), nicotianamine metabolic processes (GO:0030417), nicotianamine biosynthetic process (GO: 0030418), tricarboxylic acid (TCA) cycle metabolic process (GO:0072350), TCA biosynthetic process (GO:0072350), amine biosynthetic process (GO: 0009308) and carbohydrate metabolic process (GO: 0005975) (corrected P-value ≤ 0.01; Supplementary Table 20 and Supplementary Fig. 10).

The roots of N. sibirica have superior K⁺ retention ability,²⁶ which is carried out through K⁺ channels and transporters^2,77 that comprise a series of gene families, such as K⁺ uptake permease/high-affinity K⁺ transporter/K⁺ transporter (KUP/HAK/KT), Two-pore K⁺ channel (TPK), High-affinity K⁺ transporter (HKT), Shaker-type potassium channels (KV-like), K⁺ efflux antiporters (KEA), and Na⁺/H⁺ exchangers (NXH). Copy numbers of these gene families revealed that in N. sibirica, only the HAK5 subfamily was significantly expanded, which belongs to the KUP/HAK/KT gene family (Fig. 3A; Supplementary Table 21). Phylogenetic analysis of the HUP/HAK/KT gene family showed that the orthologs of HAK5 increased from one or two copies in other species of the Sapindales to six copies in N. sibirica, five of which were generated by tandem duplications (Fig. 3B). Three copies (LG02.979, LG02.984, and LG04.306) displayed up-regulated expression pattern under salt stress (Supplementary Fig. 11).

As an important Fe chelator in plants, nicotianamine (NA) is vital for Fe deficiency toleration in high pH soils, such as saline–alkali soils.⁷⁸ In N. sibirica, the Nicotianmine Synthase (NAS) gene family was identified as significantly expanded to include 16 copies, about fourfold higher than in other species of the Sapindales (Fig. 3A and C). The phylogenetic analysis found that most of these genes were the product of tandem duplications (Fig. 3C). Notably, under salt stress, most of these NAS genes were highly expressed in roots but not in stems and leaves (Supplementary Fig. 12), suggesting that it is likely that enhanced nicotianmine biosynthesis allowed roots to absorb sufficient Fe from the soil.

The TCA cycle is one of the most fundamental processes for cell energy production.⁷⁹ It consists of nine metabolic steps and eight enzymes, of which citrate synthase (CSY), isocitrate dehydrogenase (IDH), and α-ketoglutarate dehydrogenase (α-KGDHC) are three rate-limiting enzymes that regulate the cycling rate (Supplementary Fig. 13). A total of seven CSY genes were identified in N. sibirica, while only two or three genes were in the other Sapindales species (Fig. 3A; Supplementary Table 22). Transcriptomic analysis showed that the expression of three CSY genes was induced by salt treatment (Supplementary Fig. 13).

3.5. Salt-stressed gene expressions based on transcriptomes

We generated RNA-seq data sets across various tissues under different salt concentration treatments. We analyzed RNA from roots, stems, and leaves of N. sibirica seedlings, which were first grown in vermiculite medium (watered 1/2 Hoagland’s solution) for two months. Then seedlings were supplemented for one week with moderate (200 mM, termed as T200) and high (400 mM, termed as T400) concentrations of NaCl, respectively, to identify genes for salinity responses (see Materials and Methods). Seedlings grown continuously in vermiculite medium poured with 1/2 Hoagland’s solution were used as control. With three replicates per sample, 27 transcriptomes were generated to calculate expression profiles. Over 92.3% of quality-filtered reads were properly mapped to the N. sibirica genome with strong transcriptome correlation across tissues (Supplementary Table 4; Supplementary Fig. 14), indicating our data are viable for subsequent analyses.

Expression tendency analyses (see Materials and Methods) of genes in each tissue (root, stem, and leaf) under different salt concentrations (control, T200, and T400) yielded eight tightly grouped clusters, which display co-responsive expression patterns in each condition (Fig. 4A). Consistently up-regulated or down-regulated gene clusters in response to increased salt concentrations treatment are particularly interesting. For example, the expression of each gene in cluster_8 (including 596 genes in leaves, 554 genes in roots, and 107 genes in stems) consistently increased with increasing salt concentration (Fig. 4A). KEGG pathway enrichment analyses revealed that genes of cluster_8 in leaves were enriched in photosynthesis-related pathways (e.g. porphyrin and chlorophyll metabolism, photosynthesis—antenna proteins, and photosynthesis; Fig. 4B), indicating that photosynthesis was enhanced with the increase of salt concentration. In roots, 554 consistently up-regulated genes of cluster_8 were functionally enriched in cutin, suberine and wax biosynthesis pathways (Fig. 4B), which probably are involved in the cuticular wax formation and maintenance of cellular structural stability under high osmotic pressure.

In contrast, genes involved in plant-pathogen interaction displayed down-regulated expression patterns in both Cluster_1 and Cluster_2 as salt concentration increased (Fig. 4), suggesting this type of biological process is weakened under salt stress. From Clusters 5–8 displaying increased expression under either T200 or T400 or both conditions, we observed co-enrichment of phenylpropanoid biosynthesis and phenylalanine metabolism pathways (Fig. 4B). It has been proposed that phenylpropanoids contributed to plant response to biotic and abiotic stresses from various aspects.⁸⁰ More importantly, metabolic pathways and biosynthesis of secondary metabolites are the most enriched pathways in almost all clusters (Fig. 4B), suggesting active metabolic levels are crucial for salt stress tolerance.

The salt overly sensitive (SOS) pathway, including SOS3, SCaBP8, SOS2, and SOS1, is critical for maintaining cellular ion homeostasis through regulating Na⁺ exclusion.^1,81 Therefore, we analyzed the expression changes of related gene members under salt stress in N. sibirica (Supplementary Fig. 15). SOS1 is a key member in the SOS pathway that squeezes excess Na⁺ out of the cell.⁸¹ The expression of SOS1 slightly decreased at a moderate salt concentration (T200) but significantly increased at a high salt concentration (T400) (Supplementary Fig. 15) in both roots and leaves, suggesting the enhanced Na⁺ exclusion under high salt conditions. In addition to Na⁺ exclusion, vacuolar Na⁺ compartmentalization, mainly mediated by NXH1, a vital Na⁺/H⁺ transporter, is another mechanism to reduce excess Na⁺ in the cytoplasm.^1,2,79 We found that the expression of NXH1 increased under salt treatments and maintained at high levels in both T200 and T400 treatments (Supplementary Fig. 15).

3.6. Key transcription factors involved in response to salt stress and tissue- and concentration-dependent signalling

As TFs play central roles in triggering the expression of a series of salt-responsive genes, we analyzed the differentially expressed TFs (Fig. 5A) and their associated regulatory networks (Fig. 5B and C) in the above-identified DEGs. Among the 14 differentially expressed TF families identified, AP2/ERF-AP2 family has the highest number of TF members, all of which were significantly up-regulated, indicating the important role of salt tolerance of this family (Fig. 5A). For example, the expression level of NsiCBF4, a regulator of drought response in Arabidopsis,⁸² was remarkably increased in leaves (8.02-fold at T200) and roots (9.78-fold at T200 and 7.51-fold at T400). Elevated expressions were also found in NsiERF109, a factor involved in retarding programmed cell death under salt stress,⁸³ and NsiERF040, a member involved in generating thickened cell walls.⁸⁴ Several TFs related to abscisic acids (ABA), such as NsiERF030, NsiWRKY33, and NsiWRKY46, were up-regulated under salt stress.

In addition, we found that seven differentially expressed TFs in N. sibirica bears the signatures of positive selection (Fig. 5B), compared with the other Sapindales species distributed in subtropical and tropical environments (see Materials and Methods). These include three functionally related to abiotic stress response (NsiSZF1, NsiSCL13, and NsiWRKY40) and four unrelated or unknown function TFs (NsiNAC036, NsiDof1, NsiHB16, and NsiMYB). By investigating the expression levels of these TFs’ orthologs in a salt-sensitive species of Sapindales, we found most of these TFs did not show differential expression under salt stress (Supplementary Table 23) except NsiSZF1 and NsiWRKY40, functionally involved in response to salt and wounding, respectively. These contrasting expression patterns provide critical evidence for the rapid evolution of TFs confer salinity adaptation of Nitrariaceae.

Moreover, we observed a significant difference in differentially expressed TFs between T200 and T400 treatments in leaves, with 21 significantly up-regulated and three TFs down-regulated under T200 treatment, while only seven TFs down-regulated and no TFs up-regulated under T400 treatment, suggesting concentration-dependent signalling in response to salt stress. In roots, the number of differentially expressed TFs was similar between T200 (13 TFs) and T400 (14 TFs), with some TFs especially expressed in T200 or T400 treatments (Fig. 5A). Weighted gene co-expression network analysis (WGCNA) further revealed that the most differentially expressed TFs were in the black and red modules (Supplementary Fig. 16). The expression of key transcription factors and their associated regulatory networks displayed tissue- and concentration-specific patterns. For example, HHO3, WRKY33, and ERF041 in the black module, bHLH35, PAT1, SCL13, and bZIP55 in the red module were only expressed in the leaves, while NAC036, ERF030 in the black module, and MYB77, ERF3 in the red module were only expressed in the roots (Fig. 5B and C). More interestingly, some TFs showed differences in expression patterns in different tissues in response to different salt treatments; for example, bHLH25, WEKY53, and WRKY46 were only expressed at a high salt concentration (T400) in roots, but at a low salt concentration (T200) in leaves, and ERF109, CBF4, SZF1, and ERF040 were expressed at both low and high salt concentrations in roots, but only at low salt concentrations in leaves (Fig. 5B and C). These results suggest that different tissues of N. sibirica may have different regulatory mechanisms in response to various salt concentrations.

4. Discussion

In the ‘published plant genome’ database (https://www.plabipd.de/plant_genomes_pa.ep, last accessed on 1 February 2023), 141 families of angiosperms have published genome data, which account for only about 1/3 of the total families of angiosperms. In this study, we presented the whole genome sequence of N. sibirica, which is the first chromosomal-level genome in the family Nitrariaceae. BUSCO, CEGMA, and LAI assessments all indicated that the N. sibirica genome reached ‘reference genome’ quality (Supplementary Tables 11–13). This genome provides a valuable resource for filling genomic gaps in flowering plants at the family level. Based on phylogenomic analyses, we confirmed that Nitrariaceae is part of Sapindales^{14,15,85–87} rather than Zygophyllales.⁸⁸ Our results also suggested that N. sibirica is the sister to all other sampled representative species of Sapindales occurring in non-arid environments (Fig. 2), providing a basis for their genomic analyses.

All species representing different lineages of Sapindales did not show WGDs (Fig. 2A and B), suggesting that polyploidization could not account for genomic changes for salinity adaptation in N. sibirica. Interestingly, many gene families decrease their copies in N. sibirica. This is also found in salt-tolerant mangroves.⁸⁹ However, many gene families involved in salt stress expanded in N. sibirica, including genes involved in potassium ion transport, TCA cycle, and nicotianamine biosynthesis, similar to those genomic changes in other salt-tolerant species.^7,12,89 These significantly expanded genes functioned at multiple levels in response to salt stress. For example, the tandem duplications of HAK5 gene (six copies; Figs. 3B and 6) in N. sibirica (only 2–3 copies in other Sapindales species) may partly explain the superior K⁺ retention ability of N. sibirica roots as revealed by previous physiological examination.^26,27 Furthermore, the tremendous expansion and synchronized expression of NAS genes may have enhanced nicotianmine biosynthesis for sufficient Fe uptake to maintain normal growth in saline-alkali soils (Fig. 3 and Supplementary Fig. 12).⁷⁶ Our results highlight the importance of gene expansion in salt stress adaptations in this and other halophytes.

We further found that many genes involved in salt stress changed their expressions. Analysis of a series of transcriptomes revealed that genes in photosynthesis-related pathways consistently increased expression with increased salt concentrations (Fig. 4B), suggesting the enhanced photosynthesis ability of N. sibirica under salt stress (Fig. 6). This phenomenon is contrasted with the observation in crops. For example, salt would lead to toxic effects on photosynthesis-related components and reduced growth rates in barley.⁹⁰ In addition to increasing photosynthetic carbon fixation, the expansion and synchronized expression of CSY genes, a key rate-limiting enzyme in the TCA cycle, may also enhance the energy supply of N. sibirica under salt stress (Figs. 3A and 6). These enhancements are partly due to the enormous energy consumption required to regulate intracellular environmental homeostasis and to maintain normal growth and development. This finding is supported by previous metabolomic evidence.²⁴

Furthermore, we found up-regulated genes related to cuticular wax formation in roots (Fig. 4B), which is favourable for cellular structural stability against high osmotic pressure caused by salt stress. Preventing the accumulation of excess Na⁺ in the cytoplasm through Na⁺ exclusion and vacuolar Na⁺ compartmentalization is critical for maintaining ion homeostasis under salt stress.^1–3 Analyzing the expression patterns of core genes involved in these two processes indicated that enhanced Na⁺ exclusion occurred under high salt stress, in addition to vacuolar Na⁺ compartmentalization, a process likely play roles in both mild and high salt stresses (Fig. 6).

Moreover, many TFs related to salinity adaptations changed their expression significantly under salt stress (Figs. 5A and 6). These TF genes are mainly from the AP2/ERF-AP2, bHLH, bZIP, C2H2, C3H, GRAS, MYB, NAC, and WRKY families, of which AP2/ERF-AP2 is likely of particular importance for salt tolerance in terms of the number and degree of differentially expressed TFs. Because of the regulatory roles of TF, the expression changes of these genes could largely initiate the responses of a series of downstream genes.⁹¹ By analyzing the expression pattern of differentially expressed TFs and their associated regulatory networks, we found tissue- and concentration-dependent signalling in response to different salt concentrations across various tissues (Fig. 5B and C). Importantly, we identified that seven of these TFs experienced rapid evolution by positive selection compared with other salt-sensitive Sapindales species (see Results), potentially contributing to expression changes and salt tolerance in N. sibirica. Such adaptive divergence of sequences within genes participating in salt adaptation was also reported in a desert poplar.⁷

In summary, our results suggest that the integration of multiple genomic and transcriptomic changes might have together contributed to the high salt tolerance of N. sibirica (as summarized in Fig. 6). From the aboveground part, increased photosynthesis and energy supply could provide more organic matter and energy to regulate intracellular environmental homeostasis and maintain normal growth and development. In the belowground part, the increased ability for K⁺ retention and Fe uptake through the expansion of key genes are crucial for maintaining cellular ion homeostasis and chlorophyll synthesis, respectively. Moreover, vacuoles of the succulent leaves acted as the main Na⁺ sink to store excess sodium ions. In addition to vacuolar Na⁺ compartmentalization, under a high level of salt stress, Na⁺ exclusion also play a significant role (Fig. 6). In addition, the expression of vital transcription factors (e.g. ERF, NAC, and bHLH) and their associated regulatory networks also play essential roles for N. sibirica to cope with salt stress. Together, our findings shed new light on the genetic bases of salt tolerance in N. sibirica, enhanced our understanding of salinity adaptation in halophytes, and pave the way for breeding salt-tolerant crops in the future.

Acknowledgements

We received support for computational work from the Big Data Computing Platform for Western Ecological Environment and Regional Development and the Supercomputing Center of Lanzhou University. We would like to thank Dr Xinxing Fu for her helpful discussions of this work and Mr Weidong Tang for his suggestions on plant material collection.

Funding

This work was supported by the National Natural Science Foundation of China (grant No. 32100170) and Fundamental Research Funds for the Central Universities (lzujbky-2021-49).

Conflict of Interest

The authors declare that they have no competing interests.

Data availability

The genomic and transcriptomic sequencing data were deposited in the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/) under the BioProject accession number PRJNA939936. The genome assembly and annotation files are available at the public figshare database under the link https://figshare.com/search?q=10.6084%2Fm9.figshare.22490920.

References

Author notes