Comparative Transcriptome Analysis of Genes Involved in GA-GID1-DELLA Regulatory Module in Symbiotic and Asymbiotic Seed Germination of Anoectochilus roxburghii (Wall.) Lindl. (Orchidaceae)

Anoectochilus roxburghii (Wall.) Lindl. (Orchidaceae) is an endangered medicinal plant in China, also called “King Medicine”. Due to lacking of sufficient nutrients in dust-like seeds, orchid species depend on mycorrhizal fungi for seed germination in the wild. As part of a conservation plan for the species, research on seed germination is necessary. However, the molecular mechanism of seed germination and underlying orchid-fungus interactions during symbiotic germination are poorly understood. In this study, Illumina HiSeq 4000 transcriptome sequencing was performed to generate a substantial sequence dataset of germinating A. roxburghii seed. A mean of 44,214,845 clean reads were obtained from each sample. 173,781 unigenes with a mean length of 653 nt were obtained. A total of 51,514 (29.64%) sequences were annotated, among these, 49 unigenes encoding proteins involved in GA-GID1-DELLA regulatory module, including 31 unigenes involved in GA metabolism pathway, 5 unigenes encoding GID1, 11 unigenes for DELLA and 2 unigenes for GID2. A total of 11,881 genes showed significant differential expression in the symbiotic germinating seed sample compared with the asymbiotic germinating seed sample, of which six were involved in the GA-GID1-DELLA regulatory module, and suggested that they might be induced or suppressed by fungi. These results will help us understand better the molecular mechanism of orchid seed germination and orchid-fungus symbiosis.


Introduction
Bioactive gibberellins (GAs) are diterpene phytohormones that modulate plant growth and development processes, including seed germination, stem elongation, leaf expansion, and reproductive development [1,2]. GA is perceived by GIBBERELLIN INSENSITIVE DWARF1 (GID1) receptors, the GA-GID1 complex is able to interact with DELLA proteins (which are negative regulators of the GA action) when binding to the receptor. The formation of the GA-GID1-DELLA complex triggers enhanced recognition of DELLA by the F-box protein SLY1/GID2, that is part of the SCF SLY1/GID2 E3 ubiquitin ligase complex, targets the DELLA repressors for degradation via ubiquitin/26S proteasome [3,4]. The GA-GID1-DELLA signaling module regulates plant growth and development by integrating multiple environmental cues and endogenous signals [5,6].

Transcriptome Profile of A. roxburghii Dry and Germinating Seeds
Total RNAs were isolated from DS, AGS and SGS (Figure 1), respectively. 3 µg of total RNA with high quality was equally pooled from each sample for sequencing library construction and they were subjected to Illumina sequencing by the HiSeq 4000 platform. In this study, all samples were sequenced in two biological replicates which were found to be highly reproducible ( Figure S1). A total of 53,140,236-56,491,050 raw reads with the length of 150 bp were generated and 41,747,128-48,108,916 clean reads were obtained (Table 1). The raw reads are available at the NCBI SRA database under the accession number PRJNA299493. The cleaned reads were used to assemble the transcripts and a total of 245,063 transcripts were generated with a mean length of 793 nt, an N50 of 1389 nt. Among the transcripts representing the same gene, the longest one was considered as unigene. Finally, we obtained 173,781 unigenes with a mean length of 653 nt and an N50 of 1060 nt. The range in unigene length was from 201 to 16,676 nt, and the length distribution of unigene is illustrated in Figure 2.
Int. J. Mol. Sci. 2015, 16, page-page 3 1389 nt. Among the transcripts representing the same gene, the longest one was considered as unigene. Finally, we obtained 173,781 unigenes with a mean length of 653 nt and an N50 of 1060 nt. The range in unigene length was from 201 to 16,676 nt, and the length distribution of unigene is illustrated in Figure 2.

The Unigenes Encoding the Proteins Involved in GA-GID1-DELLA Regulatory Module
GA-GID1-DELLA signaling module in plants is shaped by phytohormone GA, receptor GID1 and repressor DELLA, and it plays a pivotal role in regulating seed germination [18,19]. Therefore, we searched the unigenes encoding proteins involved in GA-GID1-DELLA module by searching the annotation with gene name as key words. As a result, 20 unigenes involved in GA biosynthesis pathway were found, including two unigenes encoding ent-copalyl diphosphatesynthase (CPS) and one unigene encoding ent-kaurene synthase (KS), both which convert precursor trans-geranylgeranyl diphosphate (GGPP) to ent-kaurene [20], two unigenes encoding ent-kaurene oxidase (KO), and one unigene encoding ent-kaurenoic acid oxidase (KAO) which covert ent-kaurene to intermediate GA12,  Table 2. Unigenes involved in GA-GID1-DELLA regulatory module.

Pathway
Gene Name

Number of Unigenes in Transcriptome
GA biosynthesis GA signaling downstream GA receptor GID1 5 DELLA family protein 11 F-box protein GID2 2

Differential Expression Genes at Seed Germination with or without Mycorrhiza Fungi
In order to obtain the differential expression genes (DEGs), the gene expression profiles of AGS and SGS were compared to those of DS, respectively. Meanwhile, the gene expression profiles of SGS was compared to those of AGS. A total of 8012, 16,331, and 11,881 DEGs were observed in AGS vs. DS, SGS vs. DS, and SGS vs. AGS, respectively ( Figure 4 and Figure S2). The venn diagram showed that a total of 3018, 1560, and 939 DEGs were only changed in AGS vs. DS, SGS vs. DS, and SGS vs. AGS, respectively ( Figure 4). In order to analyze their functions, these DEGs were annotated in GO. Under the GO classification, the unigenes involved in biological process, cellular component and molecular function in SGS vs. DS were very different from those of AGS vs. DS, but interestingly, the number of unigenes associated with catalytic activity is very high in both ( Figure 5A,B), indicating that these unigenes annotated to catalytic activity are important for A. roxburghii seed germination. For the GO classification of DEGs in SGS vs. AGS, the number of unigenes associated with metabolic process and catalytic activity are very high ( Figure 5C), implying that mycorrhizal fungi might influence the metabolic process by playing a role in regulating catalytic activity. When DEGs were searched against KEGG pathway, 3057 (AGS vs. DS), 7964 (SGS vs. DS), and 6586 (SGS vs. AGS) DEGs with significant hits were returned, respectively. Among these pathways, the top 20 pathways for the three groups of DEGs were summarized (Table 3), and it is noteworthy that "ubiquitin mediated proteolysis" is listed in SGS vs. DS and SGS vs. AGS, but it is not in AGS vs. DS. It is well known that the DELLA proteins are polyubiquitinated by SCF SLY1/GID2 E3 ubiquitin ligase and then proteolyzed via the 26S proteasome [22], and regulate arbuscule formation in arbuscular mycorrhizal symbiosis [23], so the results suggested that DELLA proteins also may be required in orchid symbiosis and fungi may influence on GA-GID1-DELLA module by mediating DELLA proteolysis.

Differential Expression Genes at Seed Germination with or without Mycorrhiza Fungi
In order to obtain the differential expression genes (DEGs), the gene expression profiles of AGS and SGS were compared to those of DS, respectively. Meanwhile, the gene expression profiles of SGS was compared to those of AGS. A total of 8012, 16,331, and 11,881 DEGs were observed in AGS vs. DS, SGS vs. DS, and SGS vs. AGS, respectively (Figure 4 and Figure S2). The venn diagram showed that a total of 3018, 1560, and 939 DEGs were only changed in AGS vs. DS, SGS vs. DS, and SGS vs. AGS, respectively (Figure 4). In order to analyze their functions, these DEGs were annotated in GO. Under the GO classification, the unigenes involved in biological process, cellular component and molecular function in SGS vs. DS were very different from those of AGS vs. DS, but interestingly, the number of unigenes associated with catalytic activity is very high in both ( Figure 5A,B), indicating that these unigenes annotated to catalytic activity are important for A. roxburghii seed germination. For the GO classification of DEGs in SGS vs. AGS, the number of unigenes associated with metabolic process and catalytic activity are very high ( Figure 5C), implying that mycorrhizal fungi might influence the metabolic process by playing a role in regulating catalytic activity. When DEGs were searched against KEGG pathway, 3057 (AGS vs. DS), 7964 (SGS vs. DS), and 6586 (SGS vs. AGS) DEGs with significant hits were returned, respectively. Among these pathways, the top 20 pathways for the three groups of DEGs were summarized (Table 3), and it is noteworthy that "ubiquitin mediated proteolysis" is listed in SGS vs. DS and SGS vs. AGS, but it is not in AGS vs. DS. It is well known that the DELLA proteins are polyubiquitinated by SCF SLY1/GID2 E3 ubiquitin ligase and then proteolyzed via the 26S proteasome [22], and regulate arbuscule formation in arbuscular mycorrhizal symbiosis [23], so the results suggested that DELLA proteins also may be required in orchid symbiosis and fungi may influence on GA-GID1-DELLA module by mediating DELLA proteolysis.

DEGs Involved in GA-GID1-DELLA Regulatory Module
Further analysis of DEGs between symbiotic and asymbiotic seed germination (SGS vs. AGS), six genes involved in GA-GID1-DELLA regulatory module were identified, including two genes encoding GA20ox, two genes encoding GA2ox, and two genes encoding DELLA protein ( Table 4). The results suggest that fungi mediate the GA-GID1-DELLA regulatory module through changing in expression of genes encoding GA20ox, GA2ox and DELLA protein. Table 4. DEGs (SGS vs. AGS 1 ) involved in GA-GID1-DELLA regulatory module.

Pathway
Gene Name Number Unigene ID log 2 fold 2 padj 3

The Complexity of Orchid Seed Germination
Seed germination is a comprehensive process. After absorbing water, the mature seed quickly shifts from a program of maturation to germination-driven and then prepares for seedling growth [24]. Almost all orchid seeds are extremely small (0.3-14 µg), dust-like, and contain notably few nutrient reserves [16]. Further, in the majority of orchid species, seed germination and early seedling development is largely dependent on mycorrhizal fungi under natural conditions because of the lack of nutritional reserves [25,26], thus the process of orchid seed germination in nature is complex and special, involving various processes such as colonization, plant growth stimulation, plant response to fungi, etc. According to the GO classification of DEGs ( Figure 4A,B), We can clearly see that DEGs from AGS vs. DS and SGS vs. DS involved in many biological processes, cellular components, and molecular functions, and the result also suggested that orchid seed germination is a complex, multi-stage process and requires the coordinated expression of numerous different functional genes. In addition, the DEGs were more in SGS vs. DS than in AGS vs. DS, implying that symbiotic fungi promote metabolic process in A. roxburghii seed germination and it is also a reflection of the fact that orchid seed germination under natural condition is a very complex physiological and biochemical process.

The Role of Gibberellins in Seed Germination and Unigenes Related to the GA Signaling Pathway
GA and ABA play significant but antagonistic roles in the regulation of the seed germination process. Whereas GA promotes seed germination, ABA inhibits the process. Furthermore, the antagonistic relationship and the ratio between these two hormones regulate the germination process [27]. The GA-deficient Arabidopsis thaliana mutant ga1-3 is defective in seed germination because of containing greatly reduced levels of bioactive GAs [28,29]. In recent years, major progress has been made in understanding the molecular mechanism of GA action, that is GA-GID1-DELLA. In brief, a soluble receptor GID1 or GID1-like elevate their direct interaction with DELLA proteins in the presence of bioactive GAs, the GA-GID1-DELLA complex is recognized by the GA-specific F-box proteins (GID2 in Oryza sativa and SLY1 in A. thaliana) for ubiquitination of the DELLA protein, which is then degraded through the 26S proteasome pathway. The degradation of DELLAs will relieve the DELLA-mediated growth restraint [30][31][32]. In this study, 49 unigenes involved in the GA-GID1-DELLA regulatory module were identified (Table 2), and the number of unigenes identified is similar to that obtained from lettuce seed (43 unigenes involved in GA-GID1-DELLA module), which implies that the assembled transcriptome is a reliable one and might be a very useful resource for gene annotation and discovery. Among the 49 gene set, there are two unigenes with better similarity to the rice GID2 than the Arabidopsis SLY1. This is not surprising since orchids are monocotyledonous species and thus more closely related to rice than Arabidopsis. GA can promote seed germination in opposition to ABA by enhancing the proteasome-dependent degradation of repressor DELLA and many genes involved in the regulation of the GA signaling pathway, which mainly divided into three categories: genes involved in GA biosynthesis, catabolism, and GA signaling downstream (receptor and repressor).

Symbiotic Fungi Influence on the GA-GID1-DELLA Regulatory Module
The GA-GID1-DELLA growth regulatory mechanism enables a flexible response to environmental variability, such as light, high salinity, temperatures, and pathogens [33]. Symbiotic fungi as a biotic factor have an effect on the GA-GID1-DELLA regulatory module. The gene expression profiles of SGS were compared to those of AGS, six DEGs involved in GA-GID1-DELLA module were identified (Table 4). Of the six DEGs, two unigenes encoded a GA20ox, two unigenes encoded a GA2ox, and two unigenes encoded a DELLA protein. It is widely accepted, and there is abundant evidence that GA20ox, GA3ox, and GA2ox catalyze key regulatory steps in GA biosynthesis [20], GA20ox, and GA3ox that catalyze the final steps in the synthesis of bioactive GAs. However, GA2ox converts bioactive GA or their precursors to inactive forms by 2β-hydroxylation [4,34,35]. Our results also proved this point and suggested that fungi have effect on GA biosynthesis through the modulation of the expression of members of two small families, GA20ox and GA2ox.
DELLA proteins function as repressors in GA signaling and are highly conserved in Arabidopsis, and in several crop plants, including barley, grape, maize, rice, and wheat. A single DELLA protein gene SLENDER RICE1 (SLR1) and SLENDER1 (SLN1) is present in rice and barley, respectively. Surprisingly, five DELLA protein genes (GA-Insensitive (GAI), Repressor of ga1-3 (RGA), RGA-like1 (RGL1), RGL2, and RGL3) have been identified in Arabidopsis, with RGL2 playing a major role in GA-dependent germination of Arabidopsis seeds. [36,37]. Three independent studies have now shown that DELLA proteins are required for AM development [23,38,39]. The results illustrate two DEGs (from SGS vs. AGS) encoding DELLA protein SLR1 and suggested DELLA proteins play an important role in orchid symbiotic seed germination, meanwhile, are also required for orchid-fungus symbiosis formation. All in all, symbiotic fungi might have effect on GA-GID1-DELLA regulatory module only through the modulation of the expression of GA20ox, GA2ox, and DELLA, but the specific mechanism still remains undefined.

Symbiotic and Asymbiotic Germination of Anoectochilus roxburghii Seeds
Mature, undehiscent A. roxburghii capsules were collected from Yongchun county in Fujian province, southeast China. The capsules were surface sterilized in 75% ethanol for 1 min and 2.5% NaClO for 5 min, then rinsed three times in sterile water, followed by transferring to sterile filter paper for drying and were split longitudinally with a sterile scalpel. Lastly, axenic seeds were cleaned from the capsule debris and stored in a 1.5-mL centrifuge tube on silica gel at 4˝C. A. roxburghii seeds were identified by DNA barcoding [40][41][42].
Co-inoculation of mycorrhizal fungi and A. roxburghii seeds was performed in a 9 cm petri dish. Sowing the axenic seeds onto a square of autoclaved nylon cloth (4 cmˆ4 cm) by an autoclaved writing brush, and the nylon cloth was positioned on water agar medium (1.2% agar). An inoculum of isolate Ar-34 (unpublished data), obtained from A. roxburghii and consisting of four 3 mmˆ3 mm portions of actively growing mycelium, was placed in each petri dish. Asymbiotic seed germination was performed on 1/2MS culture media. All treatments were maintained at 25˝C with a 12 h/12 h light-dark cycle. Dry and the third-stage seeds ( Figure 1) were collected, immediately frozen in liquid nitrogen, and stored at´80˝C prior to RNA extraction.

Total RNA Extraction
Total RNAs were extracted from each sample using RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany) and treated with an RNase-free DNase I digestion kit (Aidlab, Beijing, China) to remove residual genomic DNA. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA was quantified using Qubit ® 2.0 Flurometer (Life Technologies, San Diego, CA, USA), checked for integrity on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Finally, 3 µg of total RNA with RIN (RNA integrity number) above 7.0 was equally pooled from each sample for cDNA library construction, and each sample has two biological replicates.

cDNA Library Construction and Sequencing
Library construction, quality detection and sequencing were performed by the Novogene Bioinformatics Institute (Beijing, China). Sequencing libraries were generated using NEBNext ® Ultra ™ RNA Library Prep Kit (NEB, Ipswich, MA, USA) and the index-coded samples were clustered using TruSeq PE Cluster Kit v3-cBot-HS (Illumia, San Diego, CA, USA), according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 4000 instrument with 150 bp paired-end reads.

Data Filtering and de Novo Assembly
Following cDNA library sequencing, sequencing-received raw image data were transformed by base calling into sequence data, after quality filtering, raw reads (sequence data without fungi reads) were obtained. High-quality clean reads were generated from the raw reads of each library by removing adaptor sequences, duplicated sequences, reads with unknown nucleotides ("N") larger than 10%, and low-quality reads with more than 50% of quality value ď20 bases. All the downstream analyses were based on the clean reads.
De novo assembly was applied to construct transcripts from the clean reads because of the absence of reference genomic sequences, and transcriptome assembly was accomplished using Trinity [43] with min_kmer_cov set to 2 by default. The longest transcript from each gene were designated unigenes.

Functional Annotation and Metabolic Pathway Analysis of Unigenes
The unigenes were searched against four databases, NCBI [44] non-redundant protein database (Nr), NCBI nucleotide database (Nt), as well as SwissProt protein database [45] with E value less than 10´5, and euKaryotic Ortholog Groups of proteins (KOG) database [46] with E value threshold of 10´3 by BLASTX (version 2.2.28+, NCBI, Bethesda, MA, USA). And the protein domains were predicted by HMMER 3.0 package in protein family (Pfam) database [47]. With Nr and Pfam annotation, the unigenes were annotated according to biological process, cellular component, and molecular function by Blast2GO v2.5 in Gene Ontology (GO) database [48]. Metabolic pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes Pathway database (KEGG) [49] and related software KAAS (KEGG Automatic Annotation Server, Tokyo, Japan).

Analysis of Different Expression Genes
Dry seeds (DS) and the third-stage seeds from asymbiotic germination (AGS) or symbiotic germination (SGS) were used as materials to investigate different expression genes by digital gene expression profile analysis. Clean reads were mapped back onto the assembled transcriptome by RSEM (v1.2.15) software [50] with bowtie2 (mismatch 0), and readcount for each gene was obtained from the mapping results. Next, differential expression analysis was performed using the DESeq R package (v1.10.1) [51]. The resulting p values were adjusted by the method of Benjamini and Hochberg [52] and adjusted p-value of 0.05 was set as the threshold for selecting genes that differentially expressed.

Conclusions
A better understanding of germination and the role for mycobiont in germination is needed for conservation of orchid populations. Due to the obligate association with a mycobiont, the germination of orchid seeds is extremely complex and varied [53]. It is a multi-stage process requiring numerous genes coordinately regulated both spatially and temporally [54].
In this study, we generated the first large-scale transcriptome dataset of the A. roxburghii seed that examined both symbiotic and asymbiotic seed germination. We identified the functions and metabolic pathways of genes expressed in the A. roxburghii seed, and focused on genes involved in the GA-GID1-DELLA regulatory module. In total, we identified 49 genes related to the regulatory module, of which six genes were differential expressed in SGS vs. AGS. The generated transcriptome dataset of A. roxburghii seed lays the foundation for further studies of the molecular mechanisms related to orchid seed germination and orchid-fungus symbiosis. Furthermore, the focus on expression profile of genes involved in the GA-GID1-DELLA module provides new insights into the molecular mechanisms of A. roxburghii seed germination.