The first released available genome of the common ice plant (Mesembryanthemum crystallinum L.) extended the research region on salt tolerance, C₃-CAM photosynthetic conversion, and halophilism

Plant materials and growth conditions

Seeds of the common ice plant (Mesembryanthemum crystallinum) were personally provided by Dr. John C. Cushman from the University of Nevada and stored under coolness and darkness until use. Originally, wild-type seeds were collected from the plants identified by Dr. Klaus Winter, an expert on the common ice plant, on a coastal cliff at the Mediterranean Sea shore close to Caesarea in Israel (around N32° 29′ 43.4″, E34° 53′ 22.8″) in 1978 (Winter et al. 1978). Three voucher specimens of M. crystallinum have been deposited in the Herbarium at the Royal Botanic Gardens Kew (55793.000, K000296094, and K000267571). In this study, our biological materials were recognized as the same plants as those specimens. Experiments, including collecting samples for this study, were conducted in compliance with relevant institutional, national, and international guidelines and laws. The seeds were aseptically sown on a medium for germination containing 4.6 g L^-1 MS salt (mixed salts for Murashige-Skoog medium), 30 g L^-1 sucrose, 1 mL^-1 B5 vitamin (Gamborg et al. 1968), 1 g L^-1 nicotinic acid, 1 g L^-1 pyridoxine hydrochloride, 10 g L^-1 thiamine hydrochlorides, and 100 g L^-1 myo-inositol), 0.80% (w/w) agarose, and pH 5.7. The raising of seedlings was performed according to the methods published by Agarie et al. (2009). The two-week-old seedlings grown in a growth chamber under 12 h of light and 12 h of darkness at 25 °C were transferred to plastic pots filled with the growth medium soils composed of 50% peat moss, 30% cocopeat, and 20% perlite, tailored for the ice plants (Japan Agricultural Cooperatives Ito-Shima, Fukuoka, Japan). The plants were irrigated with a nutrient solution of 1.5 g L^-1 OAT House No. 1, containing primary nutrients including 10% Nitrogen (1.5% as ammoniacal nitrogen and 8.2% as nitrate nitrogen), 8.0% water-soluble phosphoric acid, 27% water-soluble potassium, 4.0% water-soluble magnesium, 0.10% water-soluble manganese, 0.10% water-soluble boron, 0.18% iron, 2.0 × 10^-3% copper, 6.0 × 10^-3% zinc, and 2.0 × 10^-3% molybdenum, in addition to 1.0 g L^-1 No. 2, comprising 11% nitrogen and 23% lime (OAT Agrio Co., Ltd., Tokyo, Japan) in a greenhouse at Kyushu University for five weeks. The plants were treated with the solution including 51 mM NaCl for two weeks. Approximately 0.60 g of tissue from each leaf was collected, quickly frozen in liquid nitrogen, and stored at −80 °C.

DNA extraction, library construction, and sequencing

Total genomic DNA was extracted from the leaf tissue and purified using MagExtractor™-Plant Genome Nucleic Acid Purification Kits (Toyobo Co., Ltd., Shiga, Japan), according to the manufacturer’s instructions. The DNA samples were fragmented by sonication and used to construct short insert paired-end libraries construction using NEBNext^® Ultra™DNA Library Prep Kits for Illumina (New England Biolabs Ltd., Ipswich, MA, USA). Briefly, in the end-repair step, fragmented DNA was phosphorylated at the 5′ end and adenylated at the 3′ end. During the ligation step, full-length circulated adaptor sequences were ligated to the fragments. After adaptor cleavage, purification and size selection were performed. The indexed PCR products were taken to obtain the final sequencing libraries. The mean insert size for paired-end libraries was 300 bp. The paired-end (2 × 150 bp) sequencing was conducted on an Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA).

Clean read preparation and genome size estimation

The mean insert size was calculated using REAPR (v1.0.18) (Hunt et al. 2013), and raw paired-end sequences were filtered based on the frequency of 21-mer sequences using the program Musket (v1.1) (Liu et al. 2013). The key parameter values were as follows: musket -omulti output -inorder pair1.fastq pair2.fastq. Sequence reads that appeared rarely or abnormally frequently were removed to obtain clean read data. In the corrected reads, unique and duplicate read numbers in the corrected reads were measured using fastqc (v0.11.9) (Simon 2010). The clean data were used for an estimate of genome size as follows. K-mers were counted and exported to histogram files using jellyfish (v2.3) (Marçais and Kingsford 2011) [key parameter: jellyfish histo reads.jf]. GenomeScope2.0 (Ranallo-Benavidez et al. 2020) corresponding key parameters were applied to calculate the genome sizes using k-mers lengths of 21 and 25.

De novo genome assembly and quality evaluation

The reads were assembled using ALGA (v1.0.3; Swat et al. 2021) with the default parameter --error-rate = 0.02. long DNA fragments 1 to 10 kb in length were combined, and gaps between them were filled with unknown bases (Ns) using Redundans (v0.14a; Pryszcz and Gabaldón 2016), a software program for scaffolding, with default parameter values. The genome coverage of reads was estimated using Mosdepth program (Pedersen and Quinlan 2018). The completeness of the assembled genome was evaluated based on the content of orthologs in higher plants, using the benchmarking universal single-copy orthologs (BUSCO) program (v5.0; Manni et al. 2021). The lineage dataset was embryophyta_odb10 (creation date: 2020-09-10, number of BUSCOs: 1614). We also searched for core genes in the genome sequences of nine other plant species: Kewa caespitosa, Pharnaceum exiguum, Macarthuria australis, Solanum chaucha, Populus trichocarpa, Arabidopsis thaliana, and Oryza sativa using BUSCO. The first three species belong to the same order, Caryophyllales, to which the ice plants belong. Genome information was obtained from the NCBI (see Note 1 “Address to genome information”, Sato et al. 2023b). The number of bases, sequences, sequences in several base number ranges, and maximum base length of the final draft genome sequences was calculated using gVolante (v2.0.0) (Nishimura et al. 2017). BLASTN (v2.2.31+; McGinnis and Madden 2004) was used to investigate the number of cDNA sequences identified by transcriptome (Lim et al. 2019), and registered DNA sequences (retrieved from NCBI, last accessed February 2022) were aligned to the final assembled genome sequence.

Phylogenetic tree creation among multiple plant species using 18S ribosomal DNA sequences

The 18S ribosomal genes were extracted using barrnap (v0.9; Seemann 2018) from the obtained genome sequences of the ice plant. As comparative objectives, 25 kinds of 18S ribosomal genes from general crops (Japanese radish [Raphanus sativus], Soybean [Glycine max], Japanese trefoil [Lotus japonicus], Barrelclover [Medicago truncatula], Adzuki bean [Vigna angularis], Banana [Musa acuminata], Barley [Hordeum vulgare], Sorghum [Sorghum bicolor], Bread wheat [Triticum aestivum], Maize [Zea mays], Apple [Malus domestica], Peach [Prunus persica], Coffee tree (Arabica var.) [Coffea arabica], Coffee tree (Robusta var.) [C. canephora], Clementine [Citrus clementina], Orange [C. sinensis], Poplar, Tobacco [Nicotiana tabacum], Tomato [Solanum lycopersicum], Eggplant [S. melongena], Potato [S. tuberosum] and Grape [Vitis vinifera]) were selected using the SILVA database 138.1 (Release. 2020-08; Pruesse et al. 2007). After joining all ribosomal DNA sequences into one file, a molecular phylogenetic tree was created using implemented in NGPhylogeny.fr (Lemoine et al. 2019) (Released in 2019). SH-aLRT (Shimodaira-Hasegawa-approximate likelihood ratio test) (Shimodaira and Hasegawa 1999) was used to determine the molecular phylogenetic tree.

Detection of repetitive regions

Repetitive sequences were detected, and custom repeat libraries involving transposable elements and long terminal repeat-retro transposons were generated using RepeatModeler2 (v2.0.2; Flynn et al. 2020) and TEclass (v2.1.3; Abrusán et al. 2009). Known repeat sequences were detected and classified in the assembled genome sequence with reference to the Repbase library (Bao et al. 2015) and the custom repeat libraries, using RepeatMasker (v4.1.2-p1; Smit et al. 2013-2015). The capital letters in the genome sequences were replaced with small characters as soft masking.

Search for genomic sequences coding transfer RNA (tRNA) and micro-RNA (miRNA)

The tRNA genes were identified in the draft common ice plant genome using tRNAscan-SE2.0 (v2.0.9) (Chan et al. 2021). The tRNA data of other nine plant species—Arabidopsis, rice, tomato, poplar, horseradish, potato, grape, soybean, and coffee tree (robusta species)—were obtained from the PlantRNA database (Cognat et al. 2013). The percentages of arbitrary tRNAs against the total tRNAs in the genome were calculated and compared to the ice plants’ values with those of the other species. Smirnov-Grubbs’ outlier tests were performed to select tRNAs more significantly involved. The test statistic T was calculated using the following equation:

The miRNA loci in the genome sequence were identified using the cmscan command in Infernal (v1.1.4; Nawrocki and Eddy 2013) using Rfam.

Gene prediction

The BRAKER2 pipeline (v2.1.5; Brůna et al. 2021) was used for the prediction of genes in the common ice plant genome. Amino acid sequences were translated from the transcriptome profile reported by Lim et al. (2019) and used as additional reference data for the prediction of genes. BRAKER2 was used with the default parameters (–softmasking). The total sequences, total bases, total amino acids, and N50 were computed based on the resulting fasta-format files containing information about the genes, coding sequences, and amino acids using seqkit (v2.0.0; Shen et al. 2016) [key parameter: seqkit stats]. Protein BLAST searches (E-value < 1e-5) were conducted using DIAMOND (v2.0.13.151; Buchfink et al. 2021) against the NCBI-non-redundant protein sequences (retrieved from NCBI in March 2022), Uniprot-swissprot (retrieved in March 18), Ensemble TAIR10 (retrieved in March 2022), and NCBI poplar amino acid sequence databases (retrieved from NCBI in March 2022).

Protein domain searches

The protein domains in the genome were identified using the Pfam (v33.1) database (Mistry et al. 2021) with E-value < 1e-3, using HMMER (v3.1b2; Potter et al. 2018). The protein databases of rice, maize, and poplar from the NCBI (last accessed February 2022) were used in the domain for a detailed classification of the PKinase family, the iTAK (v18.12) web tool (Zheng et al. 2016; last accessed February 2022) was utilized. The ratio of families with a high ratio of genes to total genes in the ice plant was compared with that of the same families in the other plants. For statistical analysis, we used Smirnov-Grubbs’ outlier tests. The following equation was used to obtain the test statistic T:

Finally, BLASTP was used to compare proteins generated from the ice plant genome and those from Arabidopsis, rice, maize, and poplar and renamed TAIR10 ID. These IDs were subjected to gene ontology (GO) enrichment analysis using DAVID (updated in 2022; accessed on March 24; Sherman et al. 2022) based on a modified Fisher exact probability test with E-value < 0.05.

Genome sequencing and de novo genome assembly

Short insert reads data (300 bp; Figure S1-(A), Sato et al. 2023b) with an estimated coverage of 50.92× and a ratio of unique to duplicate reads of about 1.63:1 was obtained by removing erroneous reads of raw paired-end data from the Illumina platform (BioProject: PRJDB13817; BioSample: SAMD00508673) (Table S1). The M. crystallinum genome size was estimated to be 366 to 369 Mb, with very low heterozygosity (about 0.010%) following an analysis of the frequency of 21 and 25-mers, using GenomeScope2.0 (Figure S1-(B) and (C), Sato et al. 2023b). The M. crystallinum final draft assembly included 286 Mb in 49,782 scaffolds with a scaffold N50 of 10,562 bp (Table S2). The BUSCO tool revealed 1,509 (93.49%) of 1,614 embryophyte library core genes, with 1,223 (75.77%) of these being ‘Complete’ matches in the genome. The completeness and contiguity of the genome were greater than the shotgun assembled M. australis and S. chaucha genomes (Figure S2: Sato et al. 2023b). Around 24,081 (99.5%) of the 24,204 transcripts from the transcriptome assembly of M. crystallinum leaves (Lim et al. 2019), and all 135 DNA sequences registered in the NCBI, were aligned to the assembled genome (Supplementary Dataset S1: Sato et al. 2023c).

Phylogenetic tree based on 18S ribosomal DNA sequences from 26 kinds of plants

We performed a phylogenic analysis using 18S ribosomal DNA (rDNA) among related species. The seven types of 18S rDNA were chosen from the ice plant genome sequence. Also, the other 25 plant species’ 18S rDNA sequences (see Phylogenetic tree creation among multiple plant species using 18S ribosomal DNA sequences in Methods) were retrieved from the ribosomal RNA database in SILVA and were aggregated with the ice plant 18S rDNAs. Based on these sequences, a molecular phylogenetic tree of 18S rDNA was constructed using PhyML+SMS/One Click in NGPhylogeny.fr (Figure S3, Sato et al. 2023b). The results showed that five species of 18S rDNAs were relatively closely related to poplar’s 18S rDNA.

Search and classification of repetitive regions

In the 286.0 Mb M. crystallinum genome, 2,423 distinct repetitive sequences families, accounting for 110.9 Mb (38.8%) of the genome, were identified using custom repeat libraries and Repbase (Bao et al. 2015). This ratio was smaller than Shen et al.’s (2022) reported value (48.04%). In decreasing order of frequency, the annotated repetitive elements were unclassified 78.0 Mb (27.27%), retroelements 21.9 Mb (7.64%), long interspersed nuclear elements (LINE) 12.5 Mb (4.37%), long terminal repeats (LTR) 9.35 Mb (3.27%), and simple repeats 7.26 Mb (2.54%). Some retroelements were classified into subfamilies, including L1/CIN4 12.4 Mb (4.34%) and RTE/Bov-B 0.85 Mb (0.03%) in the LINE, and Ty1/Copia 4.90 Mb (1.71%) and Gypsy/DIRS1 4.35 Mb (1.52%) in the LTR (Table 1).

Detection of tRNA and miRNA coding genes from the genome

A total of 971 tRNAs, excluding pseudogenes, were detected in the assembled genome, and were sorted into several groups based on codon designation. The codon with the most abundant tRNA was isoleucine and the least was tryptophan (Figure S4, Sato et al. 2023b). The number of tRNAs was as follows: Arabidopsis 585, rice 505, poplar 505, tomato 723, horseradish 500, potato 736, grape 391, and soybean 700. Interspecific comparisons using the Smirnov-Grubbs outlier test and focusing on these eight species indicated that the abundance of isoleucine was significantly highest and that of tryptophan was significantly lower (P < 0.05; Figure 1 and Table S3, Sato et al. 2023b).

Figure 1. Comparison of the percentage of tRNAs in 9 plant species including ice plant.

tRNAs significant differently abundant from the other 8 species by Smirnov-Grabs outlier test, are shown in black, and the other tRNAs are shown in gray. Bars indicate ice plant, Arabidopsis, rice, tomato, poplar, horseradish, potato, grape, and soybean from the left of each series. Asterisks indicate statistical significance: * P < 0.05, n = 9.

In addition, miRNAs loci were identified from the genome with reference to the Rfam database, to obtain miRNA profiling independent of their expression levels. MiRNAs are 21 to 24 nt molecules that regulate post-transcriptional mRNA modification, playing important roles in plant growth and tolerance to environmental stress. 100 miRNA loci were identified and categorized into 25 families. The RNA family with the largest number of loci was MIR169 (25), followed by mir-399 (16), MIR159 (8), and mir-166 (7). mRNAs targeted by miRNA families were predicted (Table S4). For instance, MIR169 family miRNAs were presumed to bind to mRNAs encoding nuclear factor gamma subunit A (NF-YA) (Chiang et al. 2016). Overall, 13 types of 25 miRNA families were likely to target mRNAs encoding transcription factors: MYB33, MYB65, HD-ZIP, WRKY, AP2-like, NAC, ARFs, IAR3, ARF16, OsSPL14, SPL, GRF2, and HLH. The rest of the targeted mRNAs are anticipated to have functions in processes such as miRNA maturation, mRNA cleavage, or metal binding.

Gene prediction and annotation

Genes (34,223), coding sequences (35,702), and amino acid regions (35,702) were predicted from the soft-masked draft M. crystallinum scaffolds ab initio using a homology-based pipeline in BRAKER2 using transcriptome data (Table S4). The representative value on bases showed that coding sequence regions cover at least 10.6% (30.4 Mb) of the total genome sequence. In comparison to several databases on 25 plant species’ genes registered in PGDBj (Asamizu et al. 2014; last accessed in March 2022), the ice plants’ genes were as abundant as those of Sorghum bicolor and Arabidopsis lyrate. Additionally, summarized data indicated that the M. crystallinum gene number was 16 times larger than those of S. bicolor and A. lyrate, equivalent to about 27.6% of the number of genes of Triticum aestivum (bread wheat) and 3.31-fold greater than that of Pyropia yezoensis (bangia) (Figure S5, Sato et al. 2023b). Each translated protein sequence was used in a BLASTP search with the DIAMOND program (Buchfink et al. 2021) against four kinds of protein sequence databases. In order of the proportion of homologous amino acid sequences identified, they were NCBI-non-redundant (82.59%), poplar (70.65%), TAIR10 (65.39%), and Swiss-prot (56.05%; Table 2) (Supplementary Dataset S2: Sato et al. 2023d). To simplify gene ID conversion to GO terms, the results, including TAIR ID, were used in the functional estimation.

Functional estimation and comparison of genomes

A Pfam domain search based on the Pfam (Mistry et al. 2021) database identified 3,703 domains in 23,521 (97.1%) genes. The most frequently occurring domain was the protein kinase domain (PKinase), at 2.18%, followed by a domain of unknown function (DUF) 4238 (1.85%), reverse transcriptase (RVT)_1 (1.45%), PPR domain-containing protein (PPR)_2 (1.42%), and protein tyrosine and serine/threonine kinase (PK_Tyr_Ser-Thr) (1.18%) (Supplementary Dataset S3: Sato et al. 2023e). The PKinase family was further classified into 94 kinase families using iTAK (Zheng et al. 2016). The top 30 kinase families with the largest number of ice plant genes are shown in descending order in Figure 2. Compared to the other four plant species, the proportion of 12 families was significantly higher (P < 0.05), and eight families—DUF4238, RVT_1, RVT_2, RVT_3, Retrotrans_gag_2, zf_RVT, Retrotrans_gag_3, Retrotrans_gag—contained retroelement domains that could be attributed to a retrotransposable element (Figure 3). The annotated genes were assigned to GO classifications based on TAIR ID in three groups —biological process (BP), cellular component (CC), and molecular function (MF)—and were categorized into 403 GO terms using the gene functional classification tool in the DAVID web service. The proportion of genes assigned to 94 GO terms did not differ significantly among five plant species (P > 0.05; Figure S6 to S8, Sato et al. 2023b), indicating that they are essential to plant survival. These findings confirmed that the ice plant genome constructed in this study contained conserved genes to some extent. 18 GO terms were identified only from ice plants, although the number of genes was small (Table S6), involving virus resistance, pollen tube development, and fat biosynthesis (BP); cytoplasmic vesicle and “soluble NSF attachment protein receptor” (SNARE; CC); and O-acyltransferase for transferring fatty acids (MF).

Figure 2. Top 30 Pkinase subfamilies classified in descending order of the number of genes included in them.

The family with the highest number of genes is shown in black, and the other families are shown in white.

Figure 3. Top 30 gene families obtained from amino acid sequences detected in ice plant four other plant species.

The top row for each family shows ice plant, Arabidopsis, rice, maize, and poplar. The independence of the proportion of genes belonging to a family in the ice plant is displayed using the Smirnov-Grubbs rejection test. Asterisks (*) indicate statistical significance: P < 0.05, n = 5. Independence is shown in red if the proportion is independently high in ice plant, in blue if it is low, and in gray if there is no difference.