The genome and preliminary single-nuclei transcriptome of Lemna minuta reveals mechanisms of invasiveness

Abstract

The ability to trace every cell in some model organisms has led to the fundamental understanding of development and cellular function. However, in plants the complexity of cell number, organ size, and developmental time makes this a challenge even in the diminutive model plant Arabidopsis (Arabidopsis thaliana). Duckweed, basal nongrass aquatic monocots, provide an opportunity to follow every cell of an entire plant due to their small size, reduced body plan, and fast clonal growth habit. Here we present a chromosome-resolved genome for the highly invasive Lesser Duckweed (Lemna minuta) and generate a preliminary cell atlas leveraging low cell coverage single-nuclei sequencing. We resolved the 360 megabase genome into 21 chromosomes, revealing a core nonredundant gene set with only the ancient tau whole-genome duplication shared with all monocots, and paralog expansion as a result of tandem duplications related to phytoremediation. Leveraging SMARTseq2 single-nuclei sequencing, which provided higher gene coverage yet lower cell count, we profiled 269 nuclei covering 36.9% (8,457) of the L. minuta transcriptome. Since molecular validation was not possible in this nonmodel plant, we leveraged gene orthology with model organism single-cell expression datasets, gene ontology, and cell trajectory analysis to define putative cell types. We found that the tissue that we computationally defined as mesophyll expressed high levels of elemental transport genes consistent with this tissue playing a role in L. minuta wastewater detoxification. The L. minuta genome and preliminary cell map provide a paradigm to decipher developmental genes and pathways for an entire plant.

Figure 1

Lemna minuta growth highlighting its anatomical analogs to other plants. A, L. minuta (lm5633) growing as a dark green mat of fronds in a sewage slough at Cotton Creek Park, Encinitas, CA, USA. The inset highlights the small size (∼1 mm) of several fronds on a fingertip. B, Cartoon of the L. minuta generational contribution that leads to the dense frond mat: mother (M; green), daughter (D; purple), GD (blue), great grand-daughter (GGD; orange), and great, great grand-daughter (GGGD) fronds (pink). The gray line represents how they are connected through exponential growth. C, A single M frond view with the attached D fronds that have GD, GGD, and GGGD fronds nested in the two meristematic pockets. D, A “nested doll” view of one M frond and the maturing generations of D, GD, GGD, and GGGD in the pocket. E, One interpretation of the Lemna frond is that the M frond is an axillary stem that has two poaches or bracts where the D frond is attached by a stipe or internode. The D, GD, GGD, and GGGD progression then is similar to a branching structure of a generic plant, and the root-like structure at each subsequent axillary node is equivalent to an adventitious root. The arrows indicate that multiple internodes will emerge from the axillary stem over time (Landolt, 1986).

Figure 2

The highly heterozygous L. minuta (lm5633) genome resolved into 21 chromosomes. A, K-mer (k = 31) frequency plot for lm5633 reveals two peaks consistent with a high level of heterozygosity (2%). B, Assembly graph of a 40 Mb region visualized with Bandage shows both the heterozygous branches as well as repeat tangles, “hairballs.” C, High-throughput chromatin conformation capture (HiC) contact map resolving the polished lm5633 assembly into 21 chromosomes (darker red more contacts, lighter red to white less contacts). Chromosomes were sorted by size and renumbered before gene prediction.

Figure 3

Gene family expansion in lm5633 is driven by TD. A, Lm5633 chromosomes aligned to sp9509 chromosomes with syntenic blocks (gray lines) anchoring positions between the two genomes. Chromosomes are the correct ratio between one another but are not to scale between the two species. A minus sign after the number means the chromosome has been flipped for visualization purposes. B, The single MYB transcription factor LATE ELONGATED HYPOCOTYL (blue line) and PSEUDO-RESPONSE REGULATOR 7 (green line) are in tight linkage and TE fragments (orange) resulting in the region expansion in lm5633. Gray lines connect other syntenic genes (blue, forward; green reverse). C, Boric acid channel TD (green line) also shows the expansion of the lm5633 genome due to TE fragments (orange). D, multi-dimensional scaling (MDS)-based visualization of Semantic similarity between significant GO terms associated with TDs in lm5633. GO terms that are more semantically similar (shared words) will be closer together in the scatter plot. Size of the circle is the log frequency and the color (red high, and yellow low) is the log FDR.

Figure 4

Whole genome evolution shows consistent gene family contractions in lm5633 ancestry. A, Gene family contractions (red) and expansions (blue) along the phylogenetic tree leading to minimized Duckweed genomes are shown at each node and total protein family contractions and expansions for each species (right). Lm5633 shows the greatest similarity and gene family conservation with wa8730. B, Self versus self and C, lm5633 versus other Ks plots to elucidate the WGD history of lm5633. Dating is based on the mean Ks peak for each comparison using all paralogs/orthologs. Plant species used: A. trichopoda (Amborella), A. thaliana, Vitis vinifera (grape), Ananas comosus (pineapple), O. sativa (rice), Oropetium thomaeum (oropetium), Zostera marina (zostera), C. esculenta (taro), S. intermedia (si7747), S. polyrhiza (sp9509), W. australiana (wa8730), and L. minuta (lm5633).

Figure 5

snRNA-seq of clonally propagating whole lm5633 plants. A, The most abundant transcripts reveal that highly expressed genes are generally related to photosynthesis across all cell types while some highly expressed genes lack robust annotation. The box plot describes the mean fraction of reads assigned to the 269 cells, center line, and the first and third quartiles at the left and right box limits, respectively. Additional points outside whiskers represent outlier cells. B, UMAP embedding of 13 k-means clusters in this snRNA-seq dataset. Trajectory analysis (black line) reveals a complex network where multiple branches result in putative terminally fated cell types. C, Log fold change in marker genes expressed in 269 nuclei forming 13 specific clusters shows robust separation of cell types for functional analysis. Expression level (red, high; blue, low; larger circle, more cells; smaller circles, fewer cells).

Figure 6

Cell-type definitions supported by GO categories. Upset plot showing GO term uniqueness associated with each cell cluster based on annotated marker genes. The bars on top count the number of GO terms unique to a specific set of clusters defined by the dots on the bottom. For example, there are 340 GO terms associated with marker genes for cell cluster 6, root epidermis, that are not associated with any marker genes for any other cell cluster. Overall, unique GO terms are associated with individual clusters (i.e. cell types) suggesting each cell type’s marker genes have a functionally unique GO annotation profile. The strongest co-occurrence of GO terms appears between the meristem–replication and meristem-axillary meristem clusters.

Figure 7

Trajectory analysis suggests 13 cell types within the L. minuta plant. A, UMAP tree embedding of snRNA-seq clusters describing the cell types from clusters described in Figure  5. B, Cartoon of L. minuta based on (Banaszek and Musiał, 2011) coloring specific cell types per the UMAP. C, Live image of lm5633 showing the connected mother and daughter fronds, with one daughter frond visible in the pocket.

Figure 8

Metal transport and accumulation-associated genes in L. minuta are expressed in putative mesophyll cells. snRNA-seq expression profiling of metal tolerance-related genes suggest a cell-type-specific expression pattern with greater gene expression localized to green tissue types in lm5633. Each peak represents transcript expression of the respective gene within each predicted cell type. The scale (right; y-axis) shows the range of gene expression (TPM), and all of the cells are on the x-axis. For example, Lm5633.a03.5.g03640 (top) encodes BOR4 localized to mesophyll and replicating cells described further in the text.