Transcriptome Analysis of Coffee Coffea arabica L. Drought Tolerant Somatic Embryogenic Line, Mediated by Antisense Trehalase

1. Introduction

Coffee is a perennial tropical crop originated in Ethiopia, characterized by abundantly distributed rainfall [1]. Coffee depends on the environment and an increase of a few degrees of average temperature and/or short periods of drought can substantially decrease yields of quality. This could result in environmental, economic, and social problems [2]. In Brazil drought and frost decreased 25% the yield of coffee beans in 2021 and is expected reductions up to 60% [3]. This has triggered considerable tension on international markets, leading to a two-fold increase in coffee prices. According to a recent study, a 2–2.5 °C temperature increase would considerably reduce the available coffee growing area [4].

The molecular mechanisms that impact drought in coffee physiology and yield has been reported [5,6,7,8,9,10,11,12]. However, the process of identifying and utilizing these traits is lengthy. An alternative tool to overcome these problems is the plant biotechnology by means of genetic manipulation by the somatic embryogenesis (SE) process. SE occurred when an embryonic stem cell is induced from a somatic cell that differentiates into a somatic embryo (SEs), with the capacity to develop a plant containing the same genetic information as its precursor and sets the template for post-embryonic development and sculpt the adult body pattern.

SE requires the transcriptional regulation of a set of genes in response to stress mediated by plant hormones, osmotic, heavy metals, salinity, signaling elements that triggers cellular reprogramming and transformation of somatic cells into somatic embryos. The application of transcriptomic analysis has revealed a great number of differentially expressed genes (DEGs) during SE in several crops and Arabidopsis thaliana for review see: [13,14,15,16,17,18,19,20,21].

Several genes have been identified as markers in SE, SERK1, LEC1, FUS3, BBM, AGL15, WUS. In coffee, C. canephora and C. arabica L, transcriptomic and proteomic analysis have been reported to understand the regulatory networks involved in the SE process [22,23,24,25,26,27]. One molecular approach is the role of trehalose in plant surviving in extreme drought/salt environments. Trehalose accumulation in plant improves abiotic stress tolerance. Trehalose is a non-reducing disaccharide of glucose (a-D-glucopyran-osyl-a-D-glucopyranoside) that serves as a reserve metabolite in yeast and fungi [28]. This disaccharide stabilizes proteins and lipid membranes.

Trehalose is synthesized in a two-reaction process, in which trehalose-6- phosphate (T6P) is first synthesized from glucose-6-P and UDP-glucose by the enzyme trehalose phosphate synthase (TPS) and subsequently de-phosphorylated by trehalose-6P phosphatase (TPP) [29]. In plants, trehalase activity hydrolyses trehalose and maintains its concentration at low levels, to prevent detrimental effects. Trehalase is present in all organs of higher plants, with the highest activities in flowers [30].

Inhibition of trehalase activity by validamicin-A in A. thaliana, led to changes in trehalose and sucrose contents in different parts of the plant [30]. In 2004, Gámez-Escobedo [31] produced an increase in tobacco plant regeneration in osmotic medium using the alfalfa trehalase antisense coupled with the RD29 promoter of A. thaliana. The same cassette was used to produce drought tolerant maize B73 inbred line [32,33].

In the present work we used a new cassette containing the antisense C. arabica L. trehalase gene under the C. arabica RD29 promoter and the NOS terminator to genetically modified an embryogenic coffee cv typica line under osmotic stress supplemented with mannitol (0.3M) and sorbitol (0.3M) containing medium. Drought resistant SE lines were recovered and were able to regenerate plants. Resistant SE lines showed a higher rate of competence to induce secondary SE. A transcriptome analysis was performed from high competent SE line to understand the acquired molecular mechanisms to developed SE. The transcriptome analysis revealed 1,549 up-regulated (Log2 [fold change (FC)] ≥2.0) and 2,301 down-regulated (Log2 [FC] ≤ −2.0) genes. A PPI network mediated by STRING database v12.0 with high confidence was performed to understand the molecular mechanisms involved in the SE process.

2. Results

2.1. SE-RD29TAS transcriptomic-wide analysis

To identify the set of genes involved in the development of SE-RD29TAS under osmotic conditions producing high competent embryogenic lines, a transcriptomic-wide analysis was performed. We sequenced cDNA libraries constructed from two treatments: SE-RD29TAS producing high competent embryogenic lines and SE-WT, embryogenic line induced by conventional protocols [34], using the Illumina HiSeq 2000 platform. This produced a total of 409.947.036 sequence reads, encompassing 122.984 Mbases from all four cDNA libraries, SE of RD29TAS and SE-WT. On average 90.6% of the quality filter passed reads generated for all three somatic embryo samples were mapped uniquely to the reference genome and Q-40 mean of 37.9. The overall distribution was 24,081 (73%) in SE-RD29TAS and 23,675 (72%) in SE-WT annotated genes were transcriptionally activated. A total of 4,879 genes (3,850 protein coding genes and 999 TEs) are expressed in SE-RD29TAS.

To determine differentially expressed genes (DEGs), transcript levels (FPKM). A total of 3,850 DEGs were found between SE-RD29TAS and SE-WT in vitro culture with 1,549 up-regulated (Log2 [fold change (FC)] ≥2.0) and 2,301 down-regulated (Log2 [FC] ≤ −2.0) genes. It was found that 50 (3.23 %) Transcription factors (FT) encoding genes are differentially expressed between these two treatments and the majority were from TF families such as MYB (6 genes), basic helix-loop-helix (bHLH) (6), APETALA2 (AP2)/ethylene responsive element binding proteins (EREBP) (6), B3-Domain (6), Nuclear factor (14) were the main TF families observed. Additionally, several embryogenesis related genes such as the embryo lethal (CDC48A (A0A068UGL5), FUS3 (A0A068V7Y1), ABI3 (A0A068U8A0), WOX2 (A0A068UL49), NFYA8 (A0A068VIL0), NFYB9 (A0A068UXD0); development (AGL15 (A0A068V010), AHK5 (A0A068V1M6), ARR4 (A0A068TWC5), AUX1 (A0A068UD59), BBM (A0A068U6P3), EM1 (A0A068UD22), EM6 (A0A068V633), F22M8.6 (A0A068VAP1), GA3OX2 (A0A068VEF1), LEA46 (A0A068UPG4); nuclear factors; NF-YC13 (A0A068V0C8), NFYA2 (A0A068UPW6), NFYA5 (A0A068TPA3), NFYA6 (A0A068UH68), NFYA7 (A0A068V9V1), NFYA9 (A0A068VE79), NFYB6 (A0A068U7K3), others (PI4KG4 (A0A068UM62), SERK1 (A0A068TXX7), UFD1 (A0A068TNP5), VAL2 (A0A068UCW0), UBP14 (A0A068UBB6), RLK5 (A0A068TVW5), RUB1 (A0A068V111).

The highly enriched GO terms (p < 10−9) found in each main functional category is shown in Figure 1. For instance, The Gene Ontology (GO) of our transcriptome revealed 227 exclusive functions of up-regulated gene products, distributed into Molecular Function (MF) (93), Biological Process (BP) (109) and Cellular Component (CC) (25). The main products include, in MF (GO:0003777) microtubule motor activity, (GO:0005544) calcium-dependent phospholipid binding, (GO:0004799) thymidylate synthase activity, (GO:0004146) dihydrofolate reductase activity, (GO:0017056) structural constituent of nuclear pore, (GO:0016491) oxidoreductase activity, (GO:0098662) cation transmembrane transporter activity, (GO:0003899) DNA-directed 5'-3' RNA polymerase activity, (GO:0004518) nuclease activity, (GO:0004721) phosphoprotein phosphatase activity, (GO:0008233) peptidase activity, (GO:0003950) NAD+ ADP-ribosyltransferase activity, (GO:0008934) inositol monophosphate 1-phosphatase activity, (GO:0047216) inositol 3-alpha-galactosyltransferase activity. In the BP (GO:0006270) DNA replication initiation, (GO:0006260) DNA replication, (GO:0007018) microtubule-base movement, (GO:0006310) DNA recombination, (GO:0046654) tetrahydrofolate biosynthetic process, (GO:0006730) one-carbon metabolic process, (GO:0006231) dTMP biosynthetic process, (GO:0043069) negative regulation of programmed cell death, (GO:0036297) interstrand cross-link repair, (GO:0070588) calcium ion transmembrane transport, (GO:0007076) mitotic chromosome condensation, (GO:0015630) microtubule cytoskeleton, (GO:0006269) DNA replication, synthesis of RNA primer, (GO:0005992) trehalose biosynthetic process. In Cellular component (GO:0000786) nucleosome, (GO:0042555) MCM complex, (GO:0043240) fanconi anaemia nuclear complex, (GO:0005643) nuclear pore, (GO:0000015) phosphopyruvate hydratase complex, (GO:0031011) Ino80 complex, (GO:0005876) spindle microtubule, (GO:0000942) condensed chromosome outer kinetochore, (GO:0005819) spindle, (GO:0048500) signal recognition particle, (GO:0080008) Cul4-RING E3 ubiquitin ligase complex, (GO:0000796) condensin complex, (GO:0005694) chromosome.These were the main 15 functional sub-groups that showed the highest significance and gave rise to 22 main biological pathways that were up-regulated, as identified by KEGG pathways (see Supplementary Table S1).

The KEGG pathways found in this work are: metabolic pathways (77), biosynthesis of secondary metabolites (56), DNA replication (15), homologous recombination (13), glycerolipid metabolism (10), phenylpropanoid biosynthesis (15), mismatch repair (7), glycerophospholipid metabolism (8), ubiquinone and other terpenoid-quinone biosynthesis (5), biosynthesis of amino acids (12), flavonoid biosynthesis (5), stilbenoid, diarylheptanoid and gingerol biosynthesis (4), galactose metabolism (5), starch and sucrose metabolism (8), indole alkaloid biosynthesis (2), glucosinolate biosynthesis (3), phenylalanine metabolism (3), glycolysis/Gluconeogenesis (5), arginine biosynthesis (3), nucleotide excision repair (7) (see Supplementary Table S1).

The SE-RD29TAS produced three times more somatic embryos than the SE-WT (Figure 2A, B). A principal component analysis (PCA) was performed to examine the data obtained (Figure 2C). The replicates of each treatment were grouped together with to 0.92 for lines SE-WT and 0.04 of variance for lines SE-RD29TAS but separated and segregated between the embryogenic lines SE-WT vs SE-RD29TAS. Volcano plots visualize distribution of DEGs (p-values and fold changes). MA plots show log2 fold changes (y-axis) and the mean of normalized counts (x-axis) on scatter plots showing up- and down-regulated genes in the embryogenic lines SE-WT vs SE-RD29TAS (Figure 2D).

2.2. Up-regulated genes in SE-RD29TAS

From 1,549 up-regulated genes in SE-RD29TAS, 230 were tightly associated in 10 modules when using a PPI network devised in a STRING database v12.0 with high confidence (0.700) (Figure 3). The network revealed the presence of the 10 modules, Cell cycle, chromatin remodeling, somatic embryogenesis, oxidative stress, generic transcription pathway, carbon metabolism, phenylpropanoid biosynthesis, trehalose biosynthesis, proline biosynthesis and glycerolipid metabolism (Figure 3, supplementary Table S2 to S11). Out of the 230 up-regulated genes, 31 are embryo lethal, distributed across six modules. 7 of these genes are found in the somatic embryogenesis, 17 in the cell cycle, 2 in chromatin remodeling, 2 in carbon metabolism, 2 in generic transcription pathway and one in the proline biosynthesis module (Supplementary Table S12).

2.3. Down-regulated genes in SE-RD29TAS

Overall, 373 out of 2,301 down-regulated genes were present in the PPI network with medium-to-high confidence (0.600). The PPI network revealed 11 modules related to: Hormonal regulation of stem cell proliferation in the root cell niche, AUX/IAA family and auxin binding, response to auxin, abscisic acid (ABA) metabolic process, transmembrane transporter activity, phenylpropanoid biosynthesis, cell wall biogenesis, glutathione metabolism, beta-glucosidase activity, amino sugar, nucleotide sugar metabolism, jasmonic acid signaling, cytokinin signaling and ribosomal proteins (Figure 4, Figure 5, supplementary Table S13 to S23).

2.4. Up-regulated genes related to SE in SE-RD29TAS

The SE module consists in 23 genes, from them 10 have been involved during SE development; NFYA5, NFYA6, NFYA9, FUS3, NFYB9(LEC1), AGL15, SERK1, BBM, ABI3, WOX2, 13 in post-embryonic development and embryo development ending and seed dormancy; CDC48, RLK5, EM1, EM6, AUX1, LEA46, ARR4, AHK5, GA3OX2, NFYA2, NFYA7, NFYA8, NFYB6 and VAL2. The SE module interact with the cell cycle, chromatin remodeling, carbon metabolism and oxidative stress modules. 7 embryo lethal genes are present in this module (Figure 3, Figure 6, supplementary Table S2 to S11), the levels of regulation are listed in Figure 7.

2.5. Up-regulated genes related to CC in SE-RD29TAS

We found 55 up-regulated transcript genes directly related to the cell cycle module derived from SE-RD29TAS under osmotic medium. 37 are involved in DNA repair, 30 in DNA replication, 25 in check point, 9 in microtubule cytoskeleton, 4 in G0/early G1, 17 in G2/M, 13 in cell division, 6 in cell cycle mitotic, 6 in generic transcription pathway, 17 are embryo lethal (Figure 3, Figure 8, supplementary Table S2 to S11).

2.6. Up-regulated genes related to chromatin remodeling in SE-RD29TAS

The chromatin remodeling module consist in 25 genes, related to the nucleosome, chromatin organization, and epigenetic regulation of development. Two are embryo lethal [35], (Figure 3, Figure 9, supplementary Table S2 to S11).

2.7. Up-regulated genes related to trehalose biosynthesis in SE-RD29TAS

The trehalose biosynthesis module contains 9 genes, 4 are trehalose phosphate phosphatases TPPs (Figure 3, Figure 10, supplementary Table S2 to S11).

2.8. Up-regulated genes related to carbon metabolism in SE-RD29TAS

The carbon metabolism module contains 9 genes, two are embryo lethal [35]. (Figure 3, Figure 11, supplementary Table S2 to S11).

2.9. Up-regulated genes related to oxidative stress in SE-RD29TAS

The oxidative stress module contains 14 genes, 3 are heat shock proteins, 2 ascorbate peroxidase, 3 galactinol synthases, 1 catalase, (Figure 3, Figure 12, supplementary Table S2 to S11).

2.10. Up-regulated genes related to phenylpropanoid biosynthesis, generic transcription pathway, glycerol-lipid metabolism, proline module in SE-RD29TAS.

19 genes were found to be involved in secondary metabolites, from which 15 are in phenylpropanoid biosynthetic process, 5 in flavonoid biosynthesis, 3 in ubiquinone and terpenoid-quinone biosynthesis (Figure 3, supplementary Table S2 to S11). 30 genes are involved in generic transcription pathway (Figure 3, supplementary Table S2 to S11). 10 genes in glycerolipid metabolism (Figure 3, supplementary Table S2 to S11). 12 genes were found in the proline biosynthesis module, 3 in glucosinolate biosynthesis, 4 in strictosidine biosynthesis (Figure 3, supplementary Table S2 to S11).

2.11. Validation of the Transcriptome-Wide Analysis

The gene validation of the transcriptomic-wide analysis was made by selecting 19 DEG and analyzing their regulation by quantitative reverse transcription PCR (qRT-PCR), using the primers described in supplementary table S24. The results are shown in Figure 13, indicating that the values are consistent with those obtained in the transcriptomic-wide analysis.

Figure 13. Validation of the transcriptomic-wide analysis by quantitative reverse transcription PCR (qRT-PCR) of 16 DEG up-regulated genes involved in somatic embryogenesis, trehalose biosynthesis, carbon metabolism, cell cycle, and the down-regulated trehalase gene.

Figure 13. Validation of the transcriptomic-wide analysis by quantitative reverse transcription PCR (qRT-PCR) of 16 DEG up-regulated genes involved in somatic embryogenesis, trehalose biosynthesis, carbon metabolism, cell cycle, and the down-regulated trehalase gene.

Figure 14. General model to explain the development of drought tolerant SE line of coffee C. arabica mediated by the antisense trehalase gene under the RD29A promoter of C. arabica and NOS terminator.

Figure 14. General model to explain the development of drought tolerant SE line of coffee C. arabica mediated by the antisense trehalase gene under the RD29A promoter of C. arabica and NOS terminator.

3. Discussion

Considerable drought coffee economic losses due to climate change on prolonged drought periods in coffee trees areas affecting the vegetative growth, flowering and bean development have been reported [36]. In coffee C. arabica, there are few studies concerning about carbon and nitrogen metabolism, defense mechanisms against oxidative stress, hormonal regulation under drought conditions [7,37,38].

In the present work we generate a highly competent and drought tolerant SE line of C. arabica L. mediated by antisense trehalase gene. A transcriptome analysis was performed to understand the molecular mechanisms involved in the SE enhanced capability and drought tolerance.

3.1. Up-regulated genes in SE-RD29TAS

A PPI analyses carried out by the STRING data base (v12.0) with high confidence (0.700) revealed a set of 230 out of 1,549 up-regulated genes, clustered in several modules: cell cycle, chromatin remodeling, embryogenesis, oxidative stress, generic transcription pathway, carbon metabolism, phenylpropanoid biosynthesis, trehalose biosynthesis, proline biosynthesis and glycerolipid metabolism. We found a set of 23 up-regulated genes involved in the SE process that fulfill the requirements to explain the molecular mechanisms related to the activation of somatic embryos interacting with protection mechanisms, including those of oxidative and osmotic stress protection.

The overall interpretation underlying the SE-RD29TAS molecular mechanisms is as follows: The C. arabica trehalase (A0A068UUM1) was down-regulated (-6.87-Log2) through antisense technology, which positively affected the all-trehalose phosphate phosphatase (TPPA (A0A068U9V1), TPPB (A0A068TXJ5), TPPD (A0A068TYT9), and TPPG (A0A068U944)) genes activity in C. arabica for their consistent up-regulation, of 1.22-Log2. Those genes are directly involved in trehalose biosynthesis, and trehalose accumulation in plants may improve abiotic stress tolerance (Figure 3, Figure 6, Figure 10).

All mentioned TPPs interact with SUS6 (A0A068TQS5) (sucrose synthase 6). SUS6 is actively involved in callose synthesis at the site of phloem sieve-elements, and in tuber biomass [39,40]. Furthermore, SUS6 interacts with T22P22.110 (A0A068UW10, a glycosyl hydrolase family 31 protein member), and in turn T22P22.110 interacts with PHS1-3 (A0A068U3V8, an alpha-glucan phosphorylase 1). Phosphorylases are important allosteric enzymatic regulators in carbohydrate metabolism, and in cellular osmotic regulation and mutants can cause embryo arrest (Figure 3, Figure 6, Figure 10, Figure 13) [41,42].

PHS1-3 interacts with ENO2 (A0A068V643) and ENO3 (A0A068ULH0), which are enolase proteins. Enolases are involved in carbon metabolism, they act as a positive regulator of cold stress, are directly involved in senescence, reproductive, vegetative, vascular and are embryo lethal [43,44]. ENO2 interacts with F16L1.10 (A0A068TYM2, a phosphoglycerate mutase family protein member). F16L1.10 interacts with PGDH1 (A0A068TM02, (3-phosphoglycerate dehydrogenase 1 chloroplastic isoform). It is involved in the plastidial phosphorylated pathway of serine biosynthesis (PPSB), required for mature pollen development (Figure 3, Figure 6, Figure 11, Figure 13).

3.1.1. Interaction with cell cycle

It is remarkable that proteins that have central roles in the carbon metabolism have established an intricate and tight interaction network in SE-RD29TAS somatic embryos. This PPI interaction network has a last member PGDH1. The protein PGDH1 interacts with THY-1 (A0A068UIN3, a bifunctional dihydrofolate reductase-thymidylate synthase 1). It is a key enzyme in folate metabolism playing two different roles: de novo synthesis of tetrahydrofolate or recycling of the dihydrofolate released, depending on the source of dihydrofolate. THY-1 interacts with PCNA2 (A0A068UNF1, proliferating cell nuclear antigen 2). This protein is directly involved in the control of eukaryotic DNA replication (Figure 3, Figure 6, Figure 8, Figure 13).

3.1.2. PCNA2 interacts with the somatic embryogenesis module

PCNA2 interacts with NF-YC13 (A0A068V0C8). NF-Y transcription factors play crucial roles in embryogenesis, seed maturation, and SE induction [45,46,47]. Another NF-Y transcription factor directly involved in embryogenesis is NF-YB9, which was identified as LEAFY COTYLEDON1 LEC1 [48,49,50]. NF-YC13 interacts with NF-YA5 (A0A068TPA3), which in turn interacts with NF-YB6 (A0A068U7K3) and NF-YA2 (A0A068UPW6). Those transcriptions factors also interact with NF-YB9/LEC1 (A0A068UXD0). NFYA5 (A0A068TPA3), it is involved in the blue light and abscisic acid (ABA) signaling pathways. Overexpression of NF-YA5 and NF-YB6 overexpression in Arabidopsis activates somatic embryogenesis [47]. NF-YB9/LEC1 interacts with NFYA8. NFYA8 (A0A068VIL0 Nuclear transcription factor Y subunit A-8) is a transcription factor directly involved in embryo development [51], (Figure 3, Figure 6, Figure 7, Figure 13). NF-YB9/LEC1 (A0A068UXD0) interacts with AGL15 (A0A068V010, an Agamous-like MADS-box protein). AGL15 is a transcription factor involved in the negative regulation of flowering, prevents premature perianth senescence and abscission, fruits development and seed desiccation, induces the expression of DTA4, LEC2, FUS3, ABI3, AT4G38680/CSP2 and GRP2B/CSP4, promotes somatic embryo development, and stimulates SE reprogramming via histone acetylation-related mechanisms [52,53], (Figure 3, Figure 6, Figure 7, Figure 13). AGL15 interacts with VAL2. VAL2 is a transcriptional repressor of the sugar-inducible genes, is also involved in seed maturation, regulates the expression of LEC1, ABI3, and FUS3, which in turn directly impacts embryonic pathways, and regulates the transition from seed maturation to seedling growth, SE, and germination [54,55,56], (Figure 3, Figure 6, Figure 7, Figure 13). VAL2 interacts with HTR2, a Histone H3.2, a core component of nucleosome. In coffee C. arabica L. cv. Catuaí Amarelo IAC 62, it was found that the CaABI3 activity correlates with the embryogenic potential with highly expressed in embryogenic masses and expression of the VAL2 gene is increased at the end of the embryogenic process [57]. Moreover, AGL15 interacts with BBM. BBM is an AP2-like ethylene-responsive transcription factor, regulates the expression of LEC1, LEC2, FUS3, and ABI3, promotes cell proliferation, cellular differentiation, morphogenesis, embryogenesis, and somatic embryogenesis induction [16,58]. And AGL15 also interacts with SERK1 (A0A068TXX7, a somatic embryogenesis receptor kinase 1) (Figure 3, Figure 6, Figure 7, Figure 13). SERK1 regulates cell proliferation and embryogenic competence, is a central regulator of gametophyte production, regulates the brassinesteroid signaling pathway, and is highly expressed during early embryogenesis stages [59,60,61].

WOX2, AGL15 and NF-YB9/LEC1 altogether interact with FUS3 (A0A068V7Y1, a transcription regulator). FUS3 regulates late embryogenesis and embryo development, controls foliar organ identity, positively regulates the abscisic acid (ABA) synthesis, and negatively regulates gibberellin production, is positive regulator of ABI3 expression and its protein accumulation in the seed, actively regulates developmental phase transitions and lateral organ development, is an active regulator during germination [61,62,63,64], (Figure 3, Figure 6, Figure 7, Figure 13). Mutations in LEC1 and FUS3 genes caused embryo lethality due to the loss of desiccation tolerance during late seed development [62,63].

LEC1/NFYB9, NFYB6 and AGL15 interact with ABI3 (A0A068U8A0, a B3 domain-containing transcription factor) (Figure 3, Figure 6, Figure 7, Figure 13). ABI3 participates in abscisic acid-regulated gene expression during seed development and embryo development, involved in leaf and embryo degreening, regulates the transition between embryo maturation and early seedling development, rather than simply a transducer of the abscisic acid signal [64,65]. ABI3 interacts with AUX1 (A0A068UD59, an auxin transporter protein 1), GEA6 (is an Em-like protein), and GA3OX2. GEA6 is stress induced, also involved in ABA response, required for normal seed development and seed maturation processes. GEA6 interacts with LEA46 (a late embryogenesis abundant protein 46). LEA46 is involved in dehydration tolerance and in the adaptive response to water deficit, is also involved in somatic embryogenesis [66,67]. AUX1 regulates auxin delivering from the mature phloem to the root meristem via the protophloem cell files. AUX1 interacts with ARR4 (A0A068TWC5, a two-component response regulator (RR) (Figure 3, Figure 6, Figure 7, Figure 13). ARR4 actively participates in the phosphorelay signal transduction system modulating the red light signaling, is directly involved in embryogenesis through CK signaling and SAM establishment during maturation of SE [68,69]. ARR4 interacts with AHK5 (A0A068V1M6, histidine kinase 5). AHK5 transmits the stress signal through the MAPK signaling cascade, is a negative regulator of the ABA and ethylene signaling pathway inhibiting root elongation and regulates stomatal activity. GA3OX2 (Gibberellin 3-beta-dioxygenase 2) participates actively in the gibberellin synthesis, regulates vegetative growth and development, and an active regulator of embryogenesis [70,71].

SERK1 interacts with WOX2 (A0A068UL49), RLK5 (A0A068TVW5, Receptor-like protein kinase 5) and CDC48 (A0A068UGL5, cell division control protein 48 homolog A) (Figure 3, Figure 6, Figure 7, Figure 13). WOX2 is involved in embryonic development and patterning, is highly expressed early during somatic embryo development [72]. CDC48 regulates cell division, development and growth processes, is actively involved in seedling, pollen and embryo development, mutants are seedling lethal [73].

CDC48 interacts with HSP70-4 (A0A068UKG5 Heat shock 70 kDa protein 4) (Figure 3, Figure 6, Figure 7, Figure 13). HSP70-4 is a key component in the de novo synthesized proteins folding process, assist precursor proteins translocation into organelles, and are responsible for damaged protein degradation under stressful conditions, involved in seed maturation processes [74]. It is worth noting the relevant functions carried out by HSP proteins, such as appropriate protein folding, translocation assistance for protein precursors into cellular organelles, involved in leaf and siliques differentiation and proper development, in seed maturation processes, flowering, regulation of cytokinins, brassinosteroid, and ABA signaling, and regulates plant cell overall transcriptional activity [75,76,77].

3.1.3. Proline biosynthesis module

SUS6 gene interact with PAL1 (A0A068VM15, phenylalanine ammonia-lyase 1) (Figure 3, Figure 6, Figure 13). PAL1 is a key enzyme of proline metabolism catalyzing the first reaction in the biosynthesis from L-phenylalanine of a wide variety of natural products based on the phenylpropane skeleton. Also is involved directly in the SA acid biosynthesis and in the response against microbial pathogens. PAL1 (A0A068VM15) interacts with ARGAH1 (A0A068TNZ, Arginase 1, mitochondrial). ARGAH1 catalyzes the hydrolysis of L-arginine to urea and L-ornithine and regulates the urea cycle and the proline and polyamines synthesis. ARGAH1 (A0A068TNZ) interacts with P5CSB (A0A068TXS1, Delta-1-pyrroline-5-carboxylate synthase B) (Figure 3, Figure 6, Figure 13). P5CSB plays a key role in proline synthesis, is directly involved in the osmoregulation process in plants, embryo development and floral transition [78,79]. The proline has a central role for plant cell wall composition, signal transduction cascades, plant development, stem elongation, root and shoot growth, inflorescence architecture, seed development and germination, stress tolerance, modulates the cyclin genes expression, embryo formation, gametophyte development [80].

3.2. Down-regulated genes in SE-RDTAS

A PPI analyses carried out by the STRING data base (v12.0) with medium to high confidence (0.600) revealed a set of 373 out of 2,301 down-regulated genes, clustered in 11 modules (Figure 4, Figure 5, supplementary Table S13 to S23). The high competent and drought tolerant SE-RDTAS line down-regulate 59 genes related to hormones, which 16 response to auxin, 13 auxin-related interacting with ethylene, six are auxin biosynthesis-related, two auxin-responsive to YABB1-5, three auxin-responsive to SAUR genes, 11 peroxidases involved in auxin catabolism, two auxin transporters and four WUSCHEL-RELATED HOMEOBOX (WOX5, WOX11, WOX12, WOX13) induced by auxins, the TCTP1 (translationally-controlled tumor protein 1), several peroxidases, glutathione S-transferases, APXs, malate dehydrogenases, thioredoxins, and light signaling regulators [11,81]. Altogether, those proteins are directly implicated during SE induction, and auxins play an essential role in apical-basal patterning and in the embryo axis establishment (Figure 4, supplementary Table S13 to S23). Furthermore, we found five proteins related to cytokinins regulation, five related to ethylene biosynthesis, 14 related to jasmonic acid regulation, and eight to gibberellin regulation. This is in accordance with the fact that different transcription factors work in a coordinated manner to orchestrate the stem cell differentiation and maintenance for meristem development.

In summary, the higher competence in SE developed for the SE-RDTAS line might be supported in the regulation model that we propose and is shown in the Figure 14. The C. arabica antisense trehalase sequence, triggering two related processes: trehalose accumulation and trehalose biosynthesis, that in turn modulates the carbon metabolism, oxidative stress, and the cell cycle by interacting with the set of genes involved in somatic embryogenesis regulation and establishment.

4. Materials and Methods

4.1. Generation of marker-free SE-RD29TAS of coffee C. arabica L. var Typica.

Genetically modified SE-RD29TAS was developed using SE lines of coffee C. arabica var. Typica [82]. The particle bombardment PDS-1000/He device was used to deliver the purified 1.5 Kb cassette containing the antisense trehalase gene under the RD29 promoter and the NOS terminator, according to [34,83]. Selection of drought tolerant SE lines were done in the CP2 medium [34]. supplemented with mannitol (0.3M) and sorbitol (0.3M) (1,602 MPa) and 3.0 g/L gelrite pH 8. 20 embryogenic masses (1 cm², 50 mg fresh weight) containing globular to early torpedo stage were placed in the center of plates for bombardment. Bombarded SE plates were incubated at 25±2°C, under a 12/12 h photoperiod at 50 µmol/m²-s irradiance provided by fluorescent lamps T8 Phillips P32T8/TL850 combined with natural light increasing red/far red light in the spectrum. Three subcultures every two weeks onto fresh medium were applied to get resistant SE lines while non bombarded SE lines deceased growth and development and became necrotic.

4.2. RNA isolation and qPCR Analysis

SE-RD29TAS and control SE-WT, in globular-early torpedo stage were used to isolate RNA using Trizol (Invitrogen, Carlsbad, CA, USA). RNA concentration was measured by its absorbance at 260 nm, the ratio 260 nm/280 nm was assessed, and its integrity was verified by electrophoresis in agarose 2% (w/v) gels. Samples of cDNA for validations were amplified by PCR using SYBR Green qPCR (Bio-Rad, Hercules, CA, USA) in Real-Time PCR Systems (CFX Bio-Rad, Hercules, CA, USA). The expression of actin, RP29, and S24 was used as reference for calculating the relative amount of target gene expression using the 2^−ΔΔCt method [84]. qPCR analysis was based on at least three biological replicates for each sample with three technical replicates. Oligonucleotides were designed to qPCR (supplementary Table S24). In parallel the sequencing of cDNA was made in GENEWIZ, Plainfield, NJ, USA. To sequence the Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) was applied. Quality of sequence reads was assessed by the software package FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and to remove sequence adapters and low-quality bases using the software Trimmomatics [85].

4.3. Aligned and Analysis of DEG and GO in the SE-RD29TAS line of coffee C. arabica

Adapter removal was performed using the Trimmomatic v0.3.6 program [85]. RNA-seq reads were aligned RNA-seq reads to the coffee C. arabica reference genome available in Phytozome v12.1.(https://phytozome.jgi.doe.gov/pz/portal.html) with the STAR aligner v.2.5.2b. [86]. In this step, the BAM (Binary Alignment/Map) files were generated. Subsequently, a count and set of transcripts were made using the featureCounts program of the Subread v.1.5.2 package [87]. A quantification and differential analysis of the transcripts was performed using the DESeq2 v1.12.4 program. Finally, an ontology analysis was performed using Blast2GO.

4.4. PPI analysis of SE-RD29TAS

A gene network with high confidence (0.700) was performed with the STRING database v12 [88], based on C. canephora, and homologous in Arabidopsis thaliana genome present in Phytozome, NCBI, the gene identifier (Id) was made according to the UNIPROT [89], NCBI [90] database. Homologous in C. canephora greater than 60% in protein sequence with A. thaliana were considered.

5. Conclusions

Trehalase silencing in somatic embryos of coffee C. arabica allowed the expression of 1,549 genes. A PPI network made by STRING v12.0 with high confidence (0.700) revealed that 230 interact tightly in 10 gene modules.

The SE module consists in 23 genes and fulfill the requirements to develop somatic embryos interacting with the cell cycle and with the trehalose biosynthesis, carbon metabolism, oxidative stress, secondary metabolites to provide drought tolerance.

SE drought tolerant lines allowed us to get plants with enhance capability of resistant to water deprivation.