To dissect the molecular mechanisms underlying Ulva gametogenesis, the global gene expression levels of U. mutabilis were measured by RNA-Seq as a function of time. We sampled five-time points: 0 h, 6 h, 24 h, 48 h and 72 h after induction of gametogenesis (Fig. 1) under standard induction conditions (CW) and two-time points (6h and 24h after induction of gametogenesis) that omitted the washing step (CNW), with three biological replicates for each time point and condition. RNA-sequencing resulted in 215 x 106 reads, >70% mapped to the U. mutabilis genome. The total number of mapped reads for each sample and the raw reads counts for each gene were listed in Table S2 and Table S3, respectively. Transcript profiles were highly reproducible among the three biological replicates at each time point based on the Pearson R2 test and PCA analysis (Table S4, Fig. S1).
Gene expression pattern during gametogenesis
In total, 8296 distinct genes (62.2% of annotated genes) were expressed during gametogenesis, with relatively low variation in the total number of expressed genes between the time points: ranging between 7146 (0 h) and 7949 expressed genes (72 h). We identified 6056 differentially expressed genes (DEGs) between any two of the five time points (0 h, 6 h, 24 h, 48 h, 72 h) during the gametogenesis process, representing 45% of genes of the annotated genome. Each time point could be characterized by a set of genes describing the various steps of gametogenesis (Fig. 1). Hierarchical clustering grouped genes in 5 clusters for further analyses based on the stabilization of the SSE and the average silhouette width values (Fig. S2).
The identified clusters were characterized by nearly unique sets of enriched GO terms and reflected the identified phases of the gametogenesis process well (Fig. 2). Cluster 1 grouped genes expressed in the vegetative phase with decreasing expression levels during gametogenesis. Cluster 2 was the only cluster with high expression during the determination phase but low expression during the differentiation and swarming phase. Cluster 3 and 4 grouped genes with high expression during the differentiation phase. The main difference between both clusters consisted of genes of cluster 3 being downregulated during the swarming phase. In contrast, genes in clusters 4 were characterized by a high expression level during the swarming phase also. Cluster 5 contains the least number of genes, upregulated only in the swarming phase.
DEGs related to photosynthesis and carboxylic acid biosynthesis were upregulated during the determination phase. Next, genes involved in DNA replication and microtubule-based movements have a high expression during the early differentiation process of gametogenesis as well as during the gamete formation. However, the different enriched GO terms between the two clusters apply primarily to DNA replication which indicates that DNA replication is completed in the early phase of the differentiation process, before gametes formation starts. And the final swarming phase was mainly relevant to the pigment synthesis and cellular amino acid biosynthesis which was the final step for the gametes formation.
Initiation of gametogenesis versus response to fragmentation
By adding an additional control group whereby thalli were fragmented but the sporulation inhibitors were not removed (CNW, chopping no washing), we aimed to disentangle the effect of fragmentation on gene expression from the induction of gametogenesis. By analyzing DEGs between CW and CNW group at 6 h and 24 h, we identified 901 and 1137 significantly differential expressed genes (FDR<0.05, log2FC>1 or <-1) at 6 h and 24 h, respectively. DEGs included 493 up-regulated and 408 down-regulated genes at 6 h and 652 up-regulated and 485 down-regulated genes at 24 h. Nearly half of the DEGs at 6 h have the same expression profile at 24 h compared to the CNW group (Fig. 3). Of the 493 up-regulated and 408 down-regulated genes at 6 h, 206 and 219 genes showed the same expression profiles at 24 h, indicating that DEGs with a function in the determination phase tended to have continuous expression profiles during early gametogenesis. The enriched GO terms for the up-regulated DEGs at 6 h are mainly related to the photosynthesis, ATP biosynthetic process and oxidation-reduction processes (Fig. 4), which are in accordance with cluster 2 enrichment. In contrast, GO terms significantly enriched at 24 h are mainly related to the microtubule-based movement, DNA replication and cilium organization, suggesting that the DNA replication initiation for the gametogenesis starts from 24 h, consistent with the expression pattern of cluster 3 (Fig. 2). The over-represented GO terms for the down-regulated genes were related to the response to oxidative stress and oxidation-reduction process both at 6 h and 24 h after induction of gametogenesis.
To identify the possible key initiator for Ulva gametogenesis, we analyzed the differentially expressed TFs at 6 h which signifies the early responder to the washing treatment exclusively by comparing the chopping and washing (CW) group with the chopping without washing (CNW) group. We identified 12 up-regulated and 6 down-regulated transcription factors (TFs) at 6 h (Table 1). Several homologs of these TFs, as represented by the top blast hits in the Ulva-PLAZA database, are involved in gamete formation in Arabidopsis or Chlamydomonas.
Table 1. Differential expressed TFs of induced (CW) versus non-induced (CNW) thalli at 6 h. Log2FC>1: upregulated, LogFC<-1: downregulated. The homologs of the TFs of A. thaliana and C. reinhardtii are listed in the "AtGID" and "CrID" columns respectively. Gene function annotation are summarized from TAIR or Phytozome databases.
GeneID
|
AtGID
|
CrID
|
Family/Domain
|
Log2FC
|
Gene function annotation
|
UMSL012_0085
|
AT4G18770
|
CR03G00530
|
MYB
|
2.55
|
MYB98: expressed in the synergid cells, affect the female gametophytea
|
UMSL129_0009
|
N/A
|
N/A
|
MYB
|
1.85
|
N/A
|
UMSL039_0026
|
AT2G43770
|
CR12G01870
|
DEAD
|
2.21
|
FAP 52: flagellar associated proteinb
|
UMSL054_0061
|
N/A
|
N/A
|
DEAD
|
2.81
|
N/A
|
UMSL003_0532
|
AT5G58080
|
CR02G02460
|
GARP_G2-like
|
1.08
|
ARR18: member of response regulatora
|
UMSL103_0003
|
AT5G43990
|
N/A
|
SET
|
1.53
|
SUVR2: regulation of eukaryotic gene expression and chromatin structurea
|
UMSL008_0028
|
AT5G57390
|
CR01G02250
|
AP2/EREBP
|
1.49
|
AIL5: involved in germination and seedling growtha
|
UMSL032_0028
|
AT2G01830
|
CR06G00550
|
DEAD
|
1.08
|
AHK4: cytokinin-binding receptor that transduces cytokinin signalsa
|
UMSL057_0048
|
AT5G53040
|
CR03G09530
|
RWP-RK
|
1.62
|
GR: promotes zygote elongation and basal cell fatesa
|
UMSL024_0043
|
AT3G28730
|
CR12G13130
|
HMG
|
4.77
|
ATHMG: binds to the promoter of repressor of flowering: FLCa
|
UMSL033_0026
|
AT2G25170
|
CR03G03890
|
SWI/SNF_SNF2
|
1.51
|
CHD3: involved in post-germination repression of embryonic developmenta
|
UMSL003_0659
|
AT3G27730
|
CR10G06300
|
DEAD
|
1.52
|
MER3: DNA helicase required for interference-sensitive meiotic crossovera
|
UMSL151_0016
|
AT3G16770
|
CR16G05640
|
AP2
|
-1.32
|
ATEBP: ethylene response factora
|
UMSL085_0029
|
AT2G17520
|
CR08G03590
|
zf-CCCH
|
-1.27
|
ATIRE1-2: involved in the regulation of the ER stress responsive genesa
|
UMSL051_0112
|
AT5G41370
|
CR06G12570
|
DEAD
|
-1.14
|
ATXPB1: encodes a DNA repair proteina
|
UMSL003_0378
|
AT1G19270
|
CR10G01560
|
DEAD
|
-1.11
|
DA1: controls the initiation of axillary meristemsa
|
UMSL053_0033
|
AT3G58680
|
CR02G09000
|
MBF1
|
-1.03
|
FAP280: flagellar associated proteinb
|
UMSL011_0269
|
AT1G01040
|
N/A
|
DEAD
|
-1.03
|
ABNORMAL SUSPENSOR 1: encodes a Dicer homologa
|
a Access from TAIR database (https://www.arabidopsis.org/)
b Access from Phytozome database (https://phytozome.jgi.doe.gov/pz/portal.html)
We acquired the homologs of the Arabidopsis and Chlamydomonas of the upregulated TFs at 6 h by blast. With reported function of the homologs, we tried to find the conserved TFs function in plants gametogenesis. Particularly, an extensive studied TF containing RWP-RK domain aroused our interests for further characterization. The RWP-RK gene in Chlamydomonas, which is minus dominance (MID) is responsible for switching on the minus-programme and switching off the plus-programme in gamete differentiation [2]. Studies in Arabidopsis indicate that the RWP-RK TFs control cell differentiation during female gametophyte development [24]. We investigated the expression pattern of the RWP-RK TF (UMSL057_0048) in more detail by means of qRT-PCR analysis, sampling tissue fragments undergoing gametogenesis every 2 hours. At each timepoints the expression was compared against a control treatment not undergoing gametogenesis. UMSL057_0048 is up-regulated in the early determination phase and throughout the determination and differentiation phase with maximum expression levels at 6 h (Fig. 5a). More specifically, UMSL057_0048 shows a quick response to the chopping and washing treatment with slight upregulation at 2 h after induction and significantly upregulation at 4 h.
In an additional experiment we contrasted the expression at 6 h of UMSL057_0048 to fragmented thalli from which the SI were not removed (CNW) and unfragmented thalli that were washed to remove the SI in the medium and hence were partially induced (NCW). qRT-PCR results indicate UMSL057_0048 is not upregulated at 6 h in treatments where gametogenesis is not induced (CNW and NCNW). In contrast in both treatments where gametogenesis is induced (CW and NCW), UMSL057_0048 is significantly upregulated (Fig. 5b). The Ulva genome contains one additional gene with an RWP-RK domain, UMSL048_0014. The latter, however, is not expressed (CPM value ~0) during gametogenesis.