Identification of Metabolic Pathways Differentially Regulated in Somatic and Zygotic Embryos of Maritime Pine

Abstract

Embryogenesis is a complex phase of conifer development involving hundreds of genes, and a proper understanding of this process is critical not only to produce embryos with different applied purposes but also for comparative studies with angiosperms. A global view of transcriptome dynamics during pine somatic and zygotic embryogenesis is currently missing. Here, we present a genome-wide transcriptome analysis of somatic and zygotic embryos at three developmental stages to identify conserved biological processes and gene functions during late embryogenesis. Most of the differences became more significant as the developmental process progressed from early to cotyledonary stages, and a higher number of genes were differentially expressed in somatic than in zygotic embryos. Metabolic pathways substantially affected included those involved in amino acid biosynthesis and utilization, and this difference was already observable at early developmental stages. Overall, this effect was found to be independent of the line (genotype) used to produce the somatic embryos. Additionally, transcription factors differentially expressed in somatic versus zygotic embryos were analyzed. Some potential hub regulatory genes were identified that can provide clues as to what transcription factors are controlling the process and to how the observed differences between somatic and zygotic embryogenesis in conifers could be regulated.

Keywords: Pinus pinaster; conifers; gene regulatory networks; pine embryogenic lines; somatic embryogenesis; transcription factor; transcriptome; zygotic embryogenesis.

FIGURE 1

Experimental design. (A) Schematic representation of comparative developing stages of somatic: early-stage translucent (ES1), pre-cotyledonary opaque (ES2), and cotyledonary (ES3) and zygotic embryos: pre-cotyledonary (PC), early cotyledonary (EC) and cotyledonary (C). (B) Scheme of RNA extracted from samples and subsequent NGS-sequencing. (C) Differential expression analysis using bioinformatics tools.

FIGURE 2

Global differential expression analyses. Venn diagrams of differentially expressed genes in zygotic (Z) versus somatic (S) embryos in three consecutive sets of similar developmental stages during late embryogenesis: PC and ES1, EC and ES2, C and ES3, respectively. The number of genes showing a fold-change >2 between consecutive stages of development is shown. (A) Genes overrepresented in zygotic (Z) compared the somatic (S) stages are considered. (B) Genes underrepresented in the zygotic are considered.

FIGURE 3

Gene co-expression network analysis. (A) Mean expression profile of genes in each correlation module representing correlative developing stages (CPM, counts per million mapped reads). The number of genes included in each module is shown by its name. Zygotic developmental stages: PC, EC, and C represented in blue. Somatic developmental stages: ES1, ES2, and ES3 are represented in red. Data are means of three values from the RNA seq reactions. (B) Heat map showing the relationships between the modules and embryo developmental stages. The correlation, positive (light to dark red) or negative (light to dark blue) with either somatic or zygotic embryos data globally represented is indicated.

FIGURE 4

Identification of hub genes and transcription factors. (A) Hierarchical clustering and heat map of hub genes differentially expressed in somatic: ES1, ES2, and ES3 and zygotic: PC, EC, and C embryos. The heat map represents a total number of 1657 hub genes with connectivity >10% (see section “Materials and Methods”). The hierarchical clustering was performed using Ward’s minimum variance method. (B) Co-expression network interactions of TF hub genes differentially expressed during pine embryo development. TFs included are: NAC domain protein (NAC), Myb-related protein (MYB), ethylene-responsive TF (ERF), protein BEARSKIN, high mobility group B protein (HMG), histone H3.3, histone H2A, chromatin licensing and DNA replication factor 1 (CDT1), DNA polymerase (DNApol), regulator of nonsense transcripts 1-like protein (UPF1), bHLH, LINE-1 reverse transcriptase isogeny (LINE-1) and type-B response regulator (B-ARR). Cytoscape version 3.9.1 platform was used to visualize the network (Shannon et al., 2003). Square symbols mean TF, blue edges indicate negative co-expression associations and red edges indicate positive co-expression associations. Colors inside each symbol agree with the co-expression module where the gene is included (see Figure 3).

FIGURE 5

Graphical view of the B-ARR gene co-expression network. Genes included are those encoding: E2F, E2F-like protein; PHR3, protein phosphate starvation response 3; PERK10, proline-rich receptor-like protein kinase; PPI, peptidyl-prolyl cis-trans isomerase; GST, glutathione S-transferase; F-box, F-box/kelch-repeat protein; JAZ, jasmonate ZIM-domain protein; PEPC, phosphoenolpyruvate carboxylase; PP, serine/threonine-protein phosphatase; PD, prephenate dehydrogenase; ADH, alcohol dehydrogenase; HSK, homoserine kinase; FGH, S-formylglutathione hydrolase; MDH, malate dehydrogenase; CCCH, zinc finger CCCH domain; HD-Zip, homeobox-leucine zipper-like protein; RABA1f, ras-related protein; ISA1, isoamylase 1; NIR, ferredoxin–nitrite reductase; HSFA1E, heat stress TF A-1e; GST, glutathione S-transferase; ARF, auxin response factor; CAT, cationic amino acid transporter, usually multiple amino acids move in and out transporter (UMAMIT), starch synthase (SS) and telomeric repeat binding protein (TRF). Cytoscape version 3.9.1 platform was used to visualize the network (Shannon et al., 2003). Square symbols mean TF, triangles mean no-TF, blue edges indicate negative associations and red edges indicate positive associations. Color inside each node means the co-expression module where the gene is included. Gray symbols represent genes that in the B-ARR interaction network are only expressed in the cotyledonary stage (C) of the zygotic embryos.

FIGURE 6

Expression patterns of regulatory hub genes. Quantitative analysis of transcript levels for genes that are differentially expressed in zygotic embryos (Zy) at early cotyledonary (EC) and cotyledonary stages (C) versus four somatic lines, PN519, PP7, PP8, and PP10, at pre-cotyledonary opaque (ES2) and cotyledonary (ES3) stages. Genes included in the analysis were those encoding NAC black, NAC magenta, and MYB brown in reference to the co-expression module where the TF is included. Each value is the mean ± standard deviation of three biological replicates. Asterisks indicate significant differences between zygotic and somatic samples calculated using Student’s t-test (P < 0.05). The expression level for all genes was normalized using Actin-7 as a reference.

FIGURE 7

Differential expression of genes involved in arginine metabolism. Quantitative analysis of transcript levels for genes that are differentially expressed in zygotic embryos (Zy) at early cotyledonary (EC) and cotyledonary stages (C) versus four somatic lines, PN519, PP7, PP8, and PP10, at pre-cotyledonary opaque (ES2) and cotyledonary (ES3) stages. Genes included in the analysis were those encoding glutamine synthetase (GS1b), asparaginase (ASPG), argininosuccinate synthase (ASS), arginase (ARG) and pyrroline-5-carboxylate dehydrogenase (P5CDH). Each value is the mean ± standard deviation of three biological replicates. Asterisks indicate significant differences between zygotic and somatic samples calculated using Student’s t-test (P < 0.01). The expression level for all genes was normalized to that of Actin-7 as a reference gene.

FIGURE 8

Differential expression of genes involved in phenylalanine biosynthesis and utilization. Quantitative analysis of transcript levels for genes that are differentially expressed in zygotic embryos (Zy) at early cotyledonary (EC) and cotyledonary stages (C) versus four somatic lines, PN519, PP7, PP8, and PP10, at pre-cotyledonary opaque (ES2) and cotyledonary (ES3) stages. Genes included in the analysis were those encoding arogenate dehydratase (ADT-F and ADT-G isoforms) ADH, arogenate dehydrogenase; PAL, phenylalanine ammonia-lyase; 4CL, 4-coumaroyl-CoA ligase; CSE, caffeoyl shikimate esterase; and CAD, cinnamyl alcohol dehydrogenase. Each value is the mean ± standard deviation of three biological replicates. Asterisks indicate significant differences between zygotic and somatic samples calculated using Student’s t-test (P < 0.01). The expression level for all genes was normalized to that of Actin-7 as a reference gene.