Chromosome-level genomes of multicellular algal sisters to land plants illuminate signaling network evolution

Abstract

The filamentous and unicellular algae of the class Zygnematophyceae are the closest algal relatives of land plants. Inferring the properties of the last common ancestor shared by these algae and land plants allows us to identify decisive traits that enabled the conquest of land by plants. We sequenced four genomes of filamentous Zygnematophyceae (three strains of Zygnema circumcarinatum and one strain of Z. cylindricum) and generated chromosome-scale assemblies for all strains of the emerging model system Z. circumcarinatum. Comparative genomic analyses reveal expanded genes for signaling cascades, environmental response, and intracellular trafficking that we associate with multicellularity. Gene family analyses suggest that Zygnematophyceae share all the major enzymes with land plants for cell wall polysaccharide synthesis, degradation, and modifications; most of the enzymes for cell wall innovations, especially for polysaccharide backbone synthesis, were gained more than 700 million years ago. In Zygnematophyceae, these enzyme families expanded, forming co-expressed modules. Transcriptomic profiling of over 19 growth conditions combined with co-expression network analyses uncover cohorts of genes that unite environmental signaling with multicellular developmental programs. Our data shed light on a molecular chassis that balances environmental response and growth modulation across more than 600 million years of streptophyte evolution.

Figure 1:

(A) Three cells of a vegetative filament of SAG 698-1b (top) compared to one cell of a vegetative filament of SAG 698-1a (bottom, both samples 1 month old). Scale bar: 20 μm; C chloroplast; N nucleus; P pyrenoid. One-cell filament contains two chloroplasts and one nucleus. (B) Chromosome counting on light micrographs of SAG 698-1b fixed and stained with acetocarmine at prophase (0.5 months old); count was also performed in metaphase and telophase. Green dots represent the 20 chromosomes which were counted after rendering a stack of ~ 100 images (scale bar: 10 μm); see Figure S1 for the original images. (C) Confocal laser scanning image of one SAG 698-1b cell (0.5 month). Scale bar: 20 μm. (D) Chromosome-level assembly of SAG 698-1b genome. Concentric rings show chromosome (Chr) numbers, gene density (blue), repeat density (yellow), RNA-seq mapping density (log10(FPKM) (dark green), and GC% density (red). Red and green links show respectively intra- and inter-chromosomal syntenic blocks. (E) Comparison of genome properties for 13 algal and 3 land plant species. The time-calibrated species tree was built from 493 low-copy genes (all nodes supported by >97% non-parametric bootstrap; numbers at nodes are estimated divergence times (mean ± standard deviation) (see Table S1F for details). Data for bar plot can be found in Table S1I,J.

Figure 2:

Comparative genomics of 13 algal and 3 land plant genomes. (A) Gene family expansion and contraction patterns estimated by CAFE using Orthofinder-identified orthogroups and the time-calibrated phylogeny of Figure 1E. Key nodes are indicated on the tree and circles denote significant expansions and contractions (circle size reflects the number of expanded/contracted orthogroups/OGs). (B) Pfam domain enrichment for genes on the node leading to Zygnematophyceae and Embryophyta (Z&E). (C) Functional (GO) enrichment for the Z&E node. (D) Orthogroups overlap among Chlorophyta, Embryophyta, Zygnematophyceae, and other streptophyte algae. (E) Enriched GO terms in the 493 orthogroups exclusive to Zygnematophyceae and Embryophyta. (F) Pfam domain overlap among Chlorophyta, Embryophyta, Zygnematophyceae, and other streptophyte algae. (G) Exclusive Pfam domains found only in Zygnematophyceae and Embryophyta. One Pfam family WI12 was studied with phylogenetic analysis, suggesting a possible horizontal gene transfer from bacteria and expression response to stresses. (H) Pfam domain combination overlap among Chlorophyta, Embryophyta, Zygnematophyceae, and other streptophyte algae. (I) Exclusive Pfam domain combinations in Zygnematophyceae and Embryophyta.

Figure 3:

Cell wall innovations revealed by protein family analyses. (A) Heatmap of homolog presence in 76 enzyme subfamilies (rows) across 17 plant and algal genomes (except for Coleochaete scutata that we used its transcriptome). Enzyme subfamilies are grouped by polysaccharide and colors indicate their biochemical roles; phylogenetic patterns compatible with gene gain that might have involved horizontal gene transfer (HGT) are with asterisks. (B) Counts of subfamilies and gene percentages (with respect to the total annotated genes) across the 17 species. Shown in the plot is the gene percentage x 100. (C) Co-expression network of SAG 698-1b containing 25 genes (most belonging to the 76 analyzed subfamilies) involved in cell wall polysaccharide syntheses. (D) Phylogeny of GT2 across the 17 species. Major plant CesA/Csl subfamilies are labeled by the SAG 698-1b homolog and newly defined subfamilies are in red. Ten bacterial beta-glucan synthase (BgsA), and fungal mixed-linkage glucan (MLG) synthase (Tft1) homologs are included to show their relationships with plant CesA/Csl subfamilies. (E) Gene expression of 11 SAG 698-1b GT2 genes across 19 experimental conditions (3 replicates each); highly expressed genes are in red. (F) Phylogeny of GT2 with ZcCesA1 (Zci_04468) homologs retrieved by BLASTP against NCBI’ NR (E-value < 1e-10); colors follow D and >5,000 bacterial homologs from >8 phyla are collapsed (blue triangle; the three major phyla are indicated). (G) Phylogeny of GT2 with ZcCslP1 (Zci_0910) homologs retrieved by BLASTP against NCBI’s NR (E-value < 1e-40); colors follow D and 279 bacterial CesA homologs (blue triangle) and 363 fungal MLG synthase homologs (turquoise triangle) are collapsed (see Data S1–12 for details). (H) Phylogeny of GT2 with ZcCslN (Zci_08939) homologs retrieved by BLASTP against NCBI’s NR (E-value < 1e-10); colors follow D and bacterial homologs are collapsed (blue triangle); Pfam domain organization are shown on the right (see Data S1–13 for details).

Figure 4.

Gene co-expression modules in Zygnema circumcarinatum SAG 698-1b. (A) Heatmap of per-module co-occurrence frequencies among genes associated with plant-microbe (p-m) interaction, calcium signaling, phytohormone, stress, transporters, cell division, and diverse phytohormones (ethylene, cytokinin, abscisic acid/ABA, auxin/AUX, jasmonic acid/JA, salicylic acid/SA); based on 150 out of 406 total modules showing co-ocurrence of at least two functional categories. (B) Modules reflecting connectivity among signaling pathways, (C) biotic interaction, (D) ethylene, stress, and growth, (E) cytokinin, strees & growth, (F) cell division, (G) GRAS-domain containing genes. In gene networks, nodes are genes (size proportional to number of neighbors) and edges reflect co-expression (width proportional to Pearson’s correlation coefficient and colors are those of interconnected genes); numbers below indicate network number; gradient colors of the edges highlight the two dominant gene categories indicated in the KEY. The full gene co-expression results can be accessed in our online portal (https://zygnema.sbs.ntu.edu.sg/).

Figure 5.. Phylogenetic distribution of (A) proteins involved in phytohormone biosynthesis, signaling, and phenylpropanoid biosynthesis. and (B) transcription factors.

CK, cytokinin; ETH, ethylene; ABA, abscisic acid; AUX, auxin; SL, strigolactone; JA, jasmonic acid; GB, gibberellic acid; SA, salicylic acid; BR, brassinosteroids; PPP, phenylpropanoid; TF transcription factors; TR transcriptional regulators; PT putative transcription-associated proteins. For phytohormone-related proteins, homolog numbers were inferred from maximum likelihood gene family trees estimated from significant BLASTP hits (E-value<1e-6) using Arabidopsis queries. Note that the high number of homologs found in Penium margaritaceum are likely due to the large genome of 3.6 Gb and >50K annotated proteins. Transcription factors were identified by TAPscan.

Figure 6:

Microexon prediction in 16 plant and algae genomes. (A) Heatmap of 45 conserved microexon-tags predicted by MEPmodeler with default parameters. Microexon rate is the rate of true microexons among all predicted results in the cluster, e.g., green cell indicates that 100% microexons with both two flanking introns are present, red indicates all microexon sequences are parts of large exons and none of them could be considered as microexons, and the others are between 0 and 1. A gray cell indicates missing data (a microexon-tag could not be found). Numbers on the right column indicate the predicted clusters containing at least one true microexon (see Yu et al., 2022 for more detail). (B) RNA-seq evidence of the 1 bp microexon in Cluster 2. (C) RNA-seq evidence of 1 bp microexon in Cluster 2 two adjacent microexons (5 and 12 bp) in Cluster 7. In B and C, the RNA-seq of condition p881sControl2 was used; RNA-seq read depth and gene annotation are shown; the number in each intron indicates the junction reads and the arrows point to microexons. (D) Exon-intron structures of microexon-tag Cluster 7 in 14 plant genomes. The structure was predicted by relaxing the stringency in M. viride genome and by doing TBLASTX search in S. muscicola genome (all the three copies are intronless in this microexon-tag), respectively. The others are predicted with default parameters.