Morphofunctional interactions between galanin and GnRH-containing neurones in the diencephalon of the ewe. The effect of oestradiol

Online ISSN 2092-9293 Print ISSN 1976-9571

Genes & Genomics https://doi.org/10.1007/s13258-018-0687-z

RESEARCH ARTICLE

Transcriptome analysis for identifying possible gene regulations during maize root emergence and formation at the initial growth stage Sun‑Goo Hwang1 · Kyung‑Hee Kim2 · Byung‑Moo Lee2 · Jun‑Cheol Moon3  Received: 28 December 2017 / Accepted: 15 March 2018 © The Genetics Society of Korea and Springer Science+Business Media B.V., part of Springer Nature 2018

Abstract The root plays an important role during plant development and growth, i.e., the plant body maintenance, nutrient storage, absorption of water, oxygen and nutrient from the soil, and storage of water and carbohydrates, etc. The objective of this study was attempted to determine root-specific genes at the initial developmental stages of maize by using network-based transcriptome analysis. The raw data obtained using RNA-seq were filtered for quality control of the reads with the FASTQC tool, and the filtered reads were pre-proceed using the TRIMMOMATIC tool. The enriched BINs of the DEGs were detected using PageMan analysis with the ORA_FISHER statistical test, and genes were assigned to metabolic pathways by using the MapMan tool, which was also used for detecting transcription factors (TFs). For reconstruction of the co-expression network, we used the algorithm for the reconstruction of accurate cellular networks (ARACNE) in the R package, and then the reconstructed co-expression network was visualized using the Cytoscape tool. RNA-seq. was performed using maize shoots and roots at different developmental stages of root emergence (6–10 days after planting, VE) and 1 week after plant emergence (V2). A total of 1286 differentially expressed genes (DEGs) were detected in both tissues. Many DEGs involved in metabolic pathways exhibited altered mRNA levels between VE and V2. In addition, we observed gene expression changes for 113 transcription factors and found five enriched cis-regulatory elements in the 1-kb upstream regions of both DEGs. The network-based transcriptome analysis showed two modules as co-expressed gene clusters differentially expressed between the shoots and roots during plant development. The DEGs of one module exhibited gene expressional coherence in the maize root tips, suggesting that their functional relationships are associated with the initial developmental stage of the maize root. Finally, we confirmed reliable mRNA levels of the hub genes in the potential sub-network related to initial root development at the different developmental stages of VE, V2, and 2 weeks after plant emergence. Keywords  Maize · Root · Initial developmental stage · Differentially expressed genes · RNA-Seq

Introduction Byung-Moo Lee and Jun-Cheol Moon equally contributed to this article. Electronic supplementary material  The online version of this article (https​://doi.org/10.1007/s1325​8-018-0687-z) contains supplementary material, which is available to authorized users.

During plant development, the root system plays an important role in the anchorage of the plant body, absorption of water, oxygen and nutrient from the soil, storage of water and carbohydrates, etc. Thus, the root is a very important

* Byung‑Moo Lee [email protected]

2



Department of Life Science, Dongguk University-Seoul, Seoul 04620, Republic of Korea

* Jun‑Cheol Moon [email protected]

3



Agriculture and Life Sciences Research Institute, Kangwon National University, Chuncheon 24341, Republic of Korea

1



Plant Genomics Laboratory, Department of Applied Plant Sciences, Kangwon National University, Chuncheon 24341, Republic of Korea

13

Vol.:(0123456789)



organ for the physical and physiological maintenance of vascular plants. The growth and development of different root zones are affected by plant hormones and environmental stimuli. For example, auxins are essential to plant hormones that positively affect root initiation through cell elongation generated by the altering plasticity of the cell wall (Tomic et al. 1998), and they are known to play a central role in the organized cellular structure of roots. Jasmonic acid (JA) inhibits the elongation of the main root axis during seedling growth and decreases adventitious root development (Ravid et al. 1975; Staswick et al. 1992). Among the environmental factors, the effect of gravity has been mostly explored for root growth, and roots generally show positive gravitropism. In addition, gravity is generally known to affect the root cap cells (Sack 1991). In soil, chemical and/or physical properties such as soil structure, pH, salinity, and sodicity affect the growth and development of roots. For example, the optimized composition of soil structure, which is the arrangement of distinct soil particles on the basis of size, shape, and grade, leads to improved plant growth through the sufficient storage capacity of water and adequate drainage (Russell 1978). In addition, root growth is affected by soil pH that is changed by soil chemical properties such as alkalinity and aluminum concentration (Raun et al. 2010; Ryan and Kochian 1993). Therefore, understanding the molecular aspects of maize root formation is necessary to provide an available food resource for an expanding world population affected by serious threats such as environmental changes. Molecular genetic approaches have been successfully used for detecting reliable genes involved in the growth and development of roots. Especially, the genetic dissection of complex root traits has been mostly studied in maize under different environment conditions through quantitative trait locus (QTL) analysis, e.g., root formation (Mano et al. 2005), number (Cai et al. 2012), angle (Guingo et al. 1998), and length (Liu et al. 2008). In addition, QTLs for primary root length (PRL), lateral root number, lateral root density, and total length of the lateral root system in Arabidopsis have been detected in the Bay-0 × Shahdara population, and QTLs for PRL were confirmed in a heterogeneous inbred family, which is a type of near-isogenic line population (Loudet et al. 2005). The development of next-generation sequencing (NGS), which is a cost-effective and rapid molecular technique, has helped in the investigation of the genetic diversity, structural variation, and altered gene expression levels at a genome-wide scale. Recently, this technique has been widely used to determine potential genes associated with the development and growth of roots in plants. For example, the candidate gene for the QTL associated with root development was identified on chromosome 4 by using significant SNP markers generated from genotyping-by-sequencing of 384 maize inbred lines (Pace et al.

13

Genes & Genomics

2015). In Raphanus sativus L., the SNP markers derived with NGS were used for the identification of the potential gene controlling 4-methylthio-3-butenyl glucosinolate contents in the roots (Zou et al. 2013). In addition, integrated analysis using NGS such as RNA-seq has suggested that functional interpretation of large-scale data for Arabidopsis root hairs showed a high correlation between mRNA and protein abundance related to the central dogma of molecular biology (Lan et al. 2013). To identify complex transcriptional regulations, coexpression network analysis has often been used to determine the functional interactions of genes on the basis of large-scale gene expression profiling data detected using high-throughput techniques such as microarray and/or RNA-seq. for example, 8 co-expressed modules have been identified using gene regulatory network analysis, suggesting that the LBD gene family is the transcriptional repressor of genes related to maize root development in response to nitrogen stress (He et al. 2016). In addition, transcriptomic changes between maize inbred lines Mo17 and B73 have been studied in maize seedlings exposed to cadmium during root development by using a consensus network analysis, and 12 co-expressed modules exhibited specific biological processes for each line (Peng et al. 2015). To identify gene regulations during root development, comparative transcriptome analysis using the co-expression network approach has been effectively applied for plant species other than maize, such as Arabidopsis (Lin et al. 2011; Salazar-Henao et al. 2016), rice (Coneva et al. 2014), and Populus (Wei et al. 2013). Here, we focused on detecting specific genes related to the emergence and formation of maize roots at the initial developmental stages by using comparative transcriptome analysis. For this purpose, we performed RNA-seq by using shoot and root samples obtained during different developmental stages (VE and V2) and attempted to analyze significantly changed gene functions in each tissue between the VE and V2 stages via network-based transcriptome analysis. In addition, we selected a candidate module closely related to root development and then observed the reliable expression patterns of several hub genes, which showed relatively high numbers of functionally interacting genes.

Materials and methods Plant materials and RNA extraction In this study, we used the B73 inbred line, which was donated by USDA-ARS-PIRU, USA, in 2014. The maize seeds were grown in the growth room (16-h light/8-h dark cycle) at 25 °C until 2 weeks after plant emergence. All samples of the leaf and root tissues were collected during three

Genes & Genomics

developmental stages: emergence (6 to 10 days after planting; VE), 1 week after plant emergence (V2), and 2 weeks after plant emergence (V3) (O’Keeffe 2009). The samples of the VE and V2 stages were used for RNA-seq. For quantitative real-time (qRT)-PCR, all samples of the three developmental stages were used. Total RNA was isolated from the collected samples by using the TRIzol reagent (GibcoBRL, Cleveland, OH, USA), according to the manufacturer’s instructions, and then stored at − 80 °C until use.

RNA‑sequencing Total RNA of samples from the VE and V2 stages were used to construct cDNA libraries by using the Illumina TruSeq kit for RNA-seq analysis. RNA-seq was conducted using the Illumina Hiseq2000 system for paired-end sequencing with two biological replicates of each sample. The raw data obtained using RNA-seq were filtered for quality control of the reads with the FASTQC tool (http://www.bioin​forma​ tics.babra​ham.ac.uk/proje​cts/fastq​c/), and the filtered reads were pre-proceed using the TRIMMOMATIC tool (Bolger et al. 2014). The FASTQC files were aligned to a reference genome (AGPv3.29) from the Ensembl Plants database by using the TOPHAT tool (Trapnell et al. 2009), and the aligned reads were assembled to detect novel and alternative splicing transcripts by using the Cufflinks tool (Trapnell et al. 2012). Finally, we determined the differentially expressed genes (DEGs) for each sample (shoot and root) of the V2 stage in comparison with the VE stage by using the fragments per kilobase of transcript per million mapped reads of genes generated using the Cufflinks tool (Trapnell et al. 2012).

In silico analysis of DEGs The enriched BINs of the DEGs were detected using PageMan analysis with the ORA_FISHER statistical test (ORA cutoff value = 1.0) (Usadel et al. 2006), and genes were assigned to metabolic pathways by using the MapMan tool, which was also used for detecting transcription factors (TFs) (Thimm et al. 2004). To identify the enriched motifs of cisregulatory elements, 1-kb upstream promoter sequences of both DEGs were obtained from the Grassius database (Yilmaz et al. 2009), and local motif enrichment analysis was performed using the CentriMo tool with default parameters (Bailey and Machanick 2012).

Weighted gene co‑expression network analysis of DEGs To construct an integrated microarray database, we collected 451 Affymetrix CEL files for the microarray by using the NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/),

and transcript levels of the genes among different CEL files were normalized using the R package affy (see Supplementary Table 1) (Gautier et al. 2004). The normalized transcript levels were retrieved from the constructed microarray dataset, and gene expression profiling of the DEGs was used to conduct co-expression analysis with the dynamic branch cutting method of the R package WGCNA (Langfelder and Horvath 2008). For reconstruction of the co-expression network, we used the algorithm for the reconstruction of accurate cellular networks (ARACNE) in the R package (Margolin et al. 2006), and then the reconstructed co-expression network was visualized using the Cytoscape tool (Shannon et al. 2003). The enriched GO functions of each module were determined using the DAVID bioinformatics resources with a significant cutoff value (p < 0.05) (Huang et al. 2009), and the transcript levels of the DEGs detected at different developmental stages and tissues were retrieved using the GENEVESTIGATOR tool (Zimmermann et al. 2004).

QRT‑PCR of hub genes For qRT-PCR, total RNA from the samples of the three developmental stages was used to construct first-strand cDNA with the ReverTra Ace® qPCR RT Master Mix and gDNA Remover (Toyobo Co., Japan), according to the manufacturer’s instructions. qRT-PCR was performed in triplicate by using the Real-Time System (BIO-RAD™, USA), and each reaction mixture contained 5 ng of cDNA, 10 µl of SYBR® Green TOP real qPCR 2 × PreMix (Enzynomics™, Korea), and 5 pmol of each candidate gene primer at a final volume of 20  µl. The qRT-PCR program was as follows: denaturation at 95 °C for 10 min, followed by 40 cycles of denaturation at 95 °C for 10 s, annealing at 56 °C for 30 s, and extension at 72 °C for 20 s; the melting curve was checked using default settings. qRT-PCR was performed using specific primers for three candidate genes (GRMZM2G477658, GRMZM2G15536, and GRMZM2G072121) and one reference gene (Table  1). Maize 18S rRNA was used as the reference gene (Manoli et al. 2012).

Results RNA‑seq of leaf and root tissues To observe gene regulations during root emergence and formation, we conducted RNA-seq analysis by using Illumina Hiseq2000 with two independent replicated experiments for each sample collected from the VE and V2 stages of maize. The raw RNA-seq data were filtered for the redundant reads of low quality, and then the short reads were mapped to a reference genome sequence (Zea

13



Genes & Genomics

Table 1  List of primers used for qPCR Gene ID

Forward primer (5′–3′)

Reverse primer (5′–3′)

GRMZM2G477658 GRMZM2G15536 GRMZM2G072121 Reference gene 18S rRNA

5′-CCT​TAG​CCT​CAA​CCT​CAT​CC-3′ 5′-CCT​TCT​CAC​AGA​CCA​CGA​GA-3′ 5′-CTC​CTA​CCC​ATC​CAG​TCG​TT-3′

5′-GGG​CAA​GTG​AAG​TCC​TCT​G-3′ 5′-CGC​CTG​GGT​TCT​TTA​CTA​GC-3′ 5′-GGA​GAA​CCT​GAC​ACC​ACC​TT-3′

5′-CCA​TCC​CTC​CGT​AGT​TAG​CTTCT-3′

5′-CCT​GTC​GGC​CAA​GGC​TAT​ATAC-3′

mays ‘B73’; see Supplementary Table 2). We harbored numerous surviving reads generated with an average mapping rate of 93.7%. The gene expression profiles of each sample were normalized using the quantile normalization method, and similarities among the different samples were identified using the hierarchical clustering analysis with complete linkage of Euclidean distance (see Supplementary Fig. 1). The results suggest that the two independently replicated experiments for each sample were appropriately performed to study complex gene regulations during root emergence (in VE) and formation (in V2). To determine specific genes that control root development, we compared the mRNA levels of the genes in the shoot and root tissues at different stages (VE and V2). A total of 1286 genes, composed of 675 upregulated genes and 611 downregulated genes, were differentially expressed in each tissue between VE and V2 (Fig. 1a). In addition, 52 genes commonly exhibited altered gene expression levels, including 29 upregulations and 23 downregulations, in both tissues at the different developmental stages. These findings suggest the functional differences of genes involved in the plant developmental mechanisms of shoots and roots. The statistically significant BINs as over-represented functional pathways of DEGs were detected using Fisher’s exact test with a cut-off value of 1.0 (Fig. 1b). The enriched BINs of the V2 stage upregulated genes detected in the shoot tissue were Photosynthesis (PS). Light reaction, major CHO metabolismrelated functions, RNA regulation of the NAC TF, and transport of peptides and oligopeptides and several BINs such as Rubisco-related PS, cell wall degradation, metal handling, secondary and nucleotide metabolism, protein synthesis of ribosomal proteins, and protein degradation and folding were over-represented among the different V2 stage downregulated genes. In the upregulated genes detected in the root tissue during the V2 stage, BINs such as stress-related function and signaling of receptor kinases were over-represented, and PS-related DEGs involved in photosystem II exhibited a significantly small percentage of genes. The V2 stage downregulated genes detected in the root tissue exhibited enriched BINs such as amino acid and secondary metabolism and signaling of receptor kinases. The enriched BINs of DEGs detected in the

13

shoots and roots were different, indicating that the initial development of shoots and roots was promoted by a distinct mechanism. To identify reliable transcript levels of genes involved in metabolic pathways, DEGs were assigned to the corresponding BINs by using the MAPMAN tool (Thimm et al. 2004) (Fig. 2). Overall, the transcript levels of these genes were different between the shoot and root. For example, many genes involved in the light reaction, nucleotide breakdown, flavonoid synthesis (in secondary metabolism), and starch and sucrose synthesis pathways exhibited strong changes in transcript levels in the shoot during the V2 stage. In contrast, the mRNA levels of genes associated with ­SO4 assimilation, TCA, and biosynthesis of terpenes, phenylpropanoids, and phenolics (in secondary metabolism) were strongly changed in the root during the V2 stage. These results support the possible existence of a distinct mechanism underlying shoot and root development.

TFs associated with shoot and root development In all, 113 TFs were differentially expressed in both maize tissues between V2 and VE (Table 1). The highest frequency of differentially expressed TFs was exhibited in the upregulated genes detected in the root tissues between V2 and VE, while the lowest frequency of differentially expressed TFs was observed in the downregulated genes detected in the shoot tissues. In addition, we identified higher numbers of TFs in the bHLH group than in the other groups. Six TFs (GRMZM2G035465, GRMZM2G098986, GRMZM2G125653, GRMZM2G163200, GRMZM2G174834, and GRMZM2G390641) were commonly upregulated in the shoot and root tissues during the V2 stage, whereas the gene expression patterns of the GRMZM5G894568 gene were different between the shoot and root tissues during the V2 stage. To identify the cisregulatory elements that control the genes, we detected the enriched motifs in the 1-kb upstream promoter region of DEGs by using local motif enrichment analysis of the CentriMo tool with default parameters (Bailey and Machanick 2012) (Table 2). We identified five enriched motifs, namely, NAC043, DOF5.3, ATHB-6, ATHB-15, and CCA1. We suggest that the differentially expressed TFs

Genes & Genomics

Fig. 1  Characterization of differentially expressed genes (DEGs) detected in the shoot and root tissues during different developmental stages (V2 and VE). a The number of upregulated or downregulated genes between the shoot and root tissues. b Enriched BINs of DEGs

detected using the PageMan analysis. The box colors represent the average signals of both genes involved in each BIN category. (Color figure online)

binding to the five enriched motifs of the cis-regulatory elements in the promoter region may play a possible role as regulators of DEGs involved in shoot and root development (Table 3).

Functional interaction of DEGs To identify the functional interactions of developmental stage-responsive DEGs detected in the shoot and root

13



Genes & Genomics

Fig. 2  Transcript levels of DEGs involved in metabolic pathways. Dot colors represent the altered transcript levels of the genes (red, upregulated genes; blue, downregulated genes). (Color figure online)

tissues, we conducted the weighted gene co-expression network analysis by using proper scale-free topology with a soft-thresholding value of 5 (­ R2 = 0.94; slope = − 0.87; see Supplementary Fig. 2A). A total of 600 DEGs were used to construct the co-expression network and divided into five different modules as highly coexpressed gene groups, except the module Eigen (ME) for grey color as a non-significant module (see Supplementary Fig. 2B). The gene expression correlations between MEyellow and MEturquoise were higher than the others. We identified the enriched GO biological process (BP) of each module by using functional annotation clustering (Fig. 3). For example, MEturquoise as the largest module exhibited six enriched BP functions such as cellular polysaccharide and glucan biosynthetic process, hexose and glucose metabolic process, fatty acid and lipid biosynthetic process, homeostatic process, regulation of transcription, and phosphorylation, while the BP function related to the regulation of transcription was over-represented in the DEGs of MEgreen as the smallest group. The DEGs of MEblue and MEturquoise exhibited high frequencies of upregulation in both plant tissues, when compared with the number of downregulated genes. However, the module genes for MEyellow were mostly downregulated in the shoot and root during the V2 stage. Interestingly, the module genes for MEgreen and MEbrown exhibited different frequencies of altered gene regulation between the shoot and root tissues during the V2 stage, in which these genes were mostly downregulated in the roots, but not in the shoots. We hypothesize that the DEGs included in MEgreen and MEbrown may play a possible role as distinct regulators of root or shoot development during different growth stages such as VE and V2.

13

Module genes as potential regulators during initial root development To demonstrate the functional relationship of the genes related to initial root development, we surveyed the gene expression patterns of each module on the basis of numerous gene expression profiles detected from the different maize tissues by using the GENEVESTIGATOR tool (Zimmermann et al. 2004) (Fig. 4, see Supplementary Fig. 3). One of the modules, MEgreen, exhibited coherent patterns of DEGs in two gene expression datasets detected from the germination stage and root tip samples, suggesting the possible role of MEgreen genes in the initial stage, such as the germination of maize root tips (Fig. 4a). However, the other modules exhibited no obvious coherence of gene expression patterns in each dataset (see Supplementary Fig. 3). We also identified reliable transcript levels of MEgreen genes between the shoot and root, and the subnetwork exhibited distinct differences in the altered expression levels of the genes, including three potential hub genes (GRMZM2G477658, GRMZM2G15536, and GRMZM2G072121) evaluated via the numbers of functional interactions of each gene (Fig. 4b). We hypothesized that the three hub genes within the MEgreen module may play a possible role in the developmental acceleration of plant growth, such as root formation and emergence through the division and/or elongation of cells. To study the gene expression patterns of the hub genes, we further observed reliable transcript levels of the three genes (GRMZM2G477658, GRMZM2G15536, and GRMZM2G072121) at various developmental stages (VE; emergence, V2, and V3) by using qRT-PCR of the different maize tissues, which were the shoots and roots obtained

2 5 8

C2C2(Zn) CO-like C2C2(Zn) GATA​ C2H2 zinc finger

1

1 5 1

2G004996

1 1

2G062885 5G813804

2G063522

2

1 1

2G467640, 2G072274

2G175265 2G003715

2

1

1

1

Nucleosome/chromatin assembly factor Orphan family Putative transcription regulator SNF7

1 4

2G070034, 2G098986* 2G362949, 2G081557

2 6

2 2

MYB-related TF NAG domain

2G098986*, 2G059102 2G011422, 2G419239, 2G119693 2G023557 2G123667, 2G146380, 2G179049, 2G011598

3

2

1 1

2 3

2G018254 2G314064, 2G396114, 2G005624

5 6

1 3

1 1

1 2

2 1 4

1

2G161435, 2G094241

1 2

3

2G163200* 2G126946, 2G180699

2G035465* 2G061906, 2G088443, 2G062541, 2G350312, 2G009478, AC193786.3 2G174284, 2G029979, 2G000842

1 6

2

2G124524, 5G842961, 2G174834*, 2G366434 2G390641*, 2G475263

4

1

2

2G039586 2G115388, 2G035103, 2G403590 2G117007

2G378653, 2G159456, 2G030744, 2G074438, 5G649600

2G020054

2G113640, 2G024005

2G125648

2G035092 2G171994

2G355773, 2G056231, 5G864735 2G484880

2G119999, 2G076272

2G004988 5G806358

2G062218 AC202864.3, 2G044576

2G179827, 2G096171 2G366373 2G035156, 2G072820, 2G057260, 2G173534

Number Gene

Histone acetyltransferases MADS box MYBdomain

2G368909 2G055243, 2G011588, 2G002915

1 3

1

1 3

4

1

Number Gene

Down-regulation from shoot-V2NE

3

2G038783 2G163200* 2G120563, 2G400281, 2G121151

2G098904

2G174834*, 2G057386, 2G421033 2G390641* 2G060485, 2G319187 2G035465* 2G001724

1 1 3

1

1 2 1 1

3

Number Gene

Down-regulation from shoot-V2NE Up-regulation from shoot-V2NE

Histone

1 1 1 2 10

4

bZIP

C3H zinc finger General Transcription Global TF group GRAS HB

3 4 3 15

8

Number Gene

Total number Up-regulation from shoot-V2NE

ARF ARR​ Allx/IAA bHLH

AP2/EREBP

TFs

Table 2  Transcription factors in DEGs detected between V2 and VE stage

Genes & Genomics

13

13

*Genes represented the overlapping genes within differentially expressed TFs detected in root and shoot between V2 and VE

17 Total

113

34

5G894568*, 5G871347, 2G004060, 2G125653* 4

All maize gene names were described by removing “GRMZM” labelled on the front of gene ID

35

2G143765, 2G003551, 2G125653* 3

27

2G077755 2G314660 2G018721, 5G894568* 1 1 2 AC205574.3 1

2G078077 2G037128 1 1 2G142751 1

4 2 9 TCP Trihelix WRKY domain

Number Gene Number Gene Number Gene Number Gene

TFs

Table 2  (continued)

Total number Up-regulation from shoot-V2NE

Down-regulation from shoot-V2NE

Genes & Genomics

Down-regulation from shoot-V2NE Up-regulation from shoot-V2NE



from the zones of cell division and elongation (Fig. 4c). Compared to the transcript levels of the genes detected in the shoot tissues at the VE stage, the three hub genes were more strongly induced in the shoot tissues and the division zone (DZ) of the root at the VE stage than in at the V2 and V3 stages. The transcript levels of two genes (GRMZM2G477685 and GRMZM2G135536) were higher in the shoot and DZ samples of the VE stage than in the others, while the GRMZM2G072121 gene exhibited a relatively high gene expression level in the elongation zone (EZ) root tissues at the V2 stage. In addition, the three hub genes exhibited the lowest gene expression levels in both tissues at the V3 stage and seemed to decrease the expression levels over a period. The transcript levels of the three hub genes reached a peak in the DZ tissues of the roots at the VE stage, suggesting that these genes are closely related to the formation and elongation of roots through cell division in maize during the earliest developmental stage (VE). Similarly, two genes (GRMZM2G477685 and GRMZM2G135536) exhibited strong gene expression levels in the primary root tissue at the initial stage of 6 days after seeding (DAS) on the basis of the expression data detected using MaizeGDB (http://maize​gdb.org/), while the GRMZM2G135536 gene was highly induced in the leaf-related tissues at the vegetative stage from V5 to V9 (Fig. S3). These results suggest that two of the hub genes, GRMZM2G477685 and GRMZM2G135536, are especially involved in the primary root growth mechanism of maize. These findings support the hypothesis for the functional association of MEgreen module genes related to root formation and emergence during the initial stage of plant development.

Discussion During the developmental processes of plants, the root tip is morphologically defined to have three zones, the meristematic, elongation, and specialization zones, and the lateral root formation region (Schiefelbein and Benfey 1991). The meristematic zone consists of mitotic cells that confer division, polarity establishment, and determination of root cell, and the elongation zone plays an important role in root cell expansion. The root cells undergo further differentiation in the specialization zone. We attempted to identify specific regulations of genes between shoot and root tissues. Many DEGs exhibited tissue-specific gene expressions in each sample during plant development, suggesting that the distinct metabolism of shoots and roots is associated with initial developmental stages such as VE and V2. Understandably, functional categories such as PS, light reaction, and starch synthesis (major CHO metabolism) involved in photosynthetic metabolism were over-represented in the upregulated genes detected in the shoot tissues of the V2 stage

Genes & Genomics Table 3  Enriched motifs within 1 kb upstream promoter region of both DEGs ID

Name

Class

E-value

adj_p-value

Region width

Region mathches

MA1045.1

NAC043

NAC/NAM

2.30E+00

1.00E-02

835

419

MA1071.1

DOF5.3

C2H2-zinc finger factors

6.60E+00

2.90E-02

488

334

MA0953.1

ATHB-6

Homeo domain factors

7.00E+00

3.10E-02

722

463

MA1026.1

ATHB-15

Homeo domain factors

8.10E+00

3.60E-02

841

510

MA0972.1

CCA1

Tryptophan cluster factors

8.60E+00

3.80E-02

465

253

Sequence logo

Fig. 3  Co-expression network of both module genes reconstructed using ARACNE. The enriched biological process (BP) function was detected using the DAVID bioinformatics resources (p < 0.05). Circle

colors indicate the module colors, and edges represent the functional interactions among different genes. The bar colors of the histogram represent the different types of DEGs. (Color figure online)

because maize exhibits the emergence of the second leaf at the V2 stage (Fig. 1). Similarly, the altered mRNA levels of the genes associated with light reaction and sucrose-tostarch conversion during photosynthetic metabolism were confirmed using the MapMan visualization of the DEGs detected in the shoot tissues at the V2 stage (Fig. 2). With respect to the DEGs detected in the root tissues, we found enriched BINs of the NAC domain TF family only in the shoot tissues at the V2 stage. In a previous study, many NAC TFs were highly expressed to a greater extent in shoot tissues

such as stem and leaf than in root and/or flower tissues, suggesting a distinct functional relationship of those in the shoot development process (Fan et al. 2014). In addition, we found strongly changed gene expression patterns of the functional categories assigned to MapMan BINs, as compared with shoot and root tissues, e.g., lipids and secondary metabolism. In sweet marjoram varieties under normal conditions, the flavonoid content of the plants at different growth stages (early vegetative stage, late vegetative stage, and flowering stage) is different, suggesting a distinct regulation of the

13



Genes & Genomics

Fig. 4  Gene expression patterns of MEgreen module genes. a Transcript levels of DEGs detected at different developmental stages and tissues. Box colors represent the different transcript levels of the genes. b Reliable transcript levels of the MEgreen module genes in the shoot and root tissues between V2 and VE. Circle colors represent the altered transcript levels of the genes (red, upregulated genes;

blue, downregulated genes). Edges indicate the functional interactions among different genes. c Transcript levels of the three genes (GRMZM2G477658, GRMZM2G15536, and GRMZM2G072121) at various developmental stages (VE; emergence, V2, and V3) by using qRT-PCR of the different maize tissues (shoots, DZ; cell division zone and EZ; elongation zone). (Color figure online)

genes involved in the biosynthetic pathways of flavonoids (Baâtour et al. 2012). Our findings support the distinct transcriptional regulation of flavonoid biosynthetic genes at the different growth stages such as VE and V2, and the mRNA levels of these genes are more variable in the shoots than in the roots. The bHLH domain TFs exhibited the highest number (6 upregulated and 4 downregulated genes) of DEGs detected in the maize root tissues at different growth stages, suggesting their functional relationship during root development (Table 2). A bHLH TF, SPATULA (SPT), affects root growth in Arabidopsis, resulting in the size regulation of the root meristem (Makkena and Lamb 2013). Here, many WRKY TFs were more strongly induced in the shoots at the V2 stage than in the VE stage (Table 1). The functional roles of WRKY TFs have been reported in Arabidopsis, and WRKY71/EXB1-regulating RAX genes play an important role in the control of shoot branching, which affects plant yield (Guo et al. 2015). However, the distinct functions of

many TFs at the different growth stages remain unknown, and further research is required to determine the complex tissue-specific regulations of TFs at each growth stage. We found five different cis-regulatory elements over-represented in the 1-kb upstream regions of the DEGs, suggesting that the transcriptionally changed TFs at the different growth stages may regulate the target genes by using the enriched cis-regulatory elements as TF-binding sites. In previous studies, the NAC genes for NAC043 (ID: MA1045.1) have been reported in the formation of the shoot apical meristem in Arabidopsis (Kim et al. 2007), and Medicago truncatula HB1 (MtHB1) encoding homeo domains such as ATHB-6 (ID: MA0953.1) and ATHB-15 (ID: MA1026.1) acts as the regulator of genes associated with lateral root emergence; the MtHB1-overexpressing transgenic plants exhibited increased primary root length and reduced root emergence (Ariel et al. 2010).

13

Genes & Genomics

In the co-expression network, we identified two modules (MEgreen and MEbrown) that showed a clear distinction of DEG numbers between the shoots and roots at the different growth stages (VE and V2), and the DEGs of these modules commonly exhibited higher mRNA levels of the genes in the roots at the VE stage than in those at the V2 stage, suggesting that MEgreen and MEbrown may related to root growth in the initial developmental stage (VE). In particular, the induced mRNA levels of the MEgreen module genes exhibited transcriptional coherence during the germination stage and in root-tip tissues. In addition, we identified strongly increased mRNA levels of the three hub genes (GRMZM2G477685 for the lipid-transfer protein family, GRMZM2G135536 for the cytochrome P450 family, and GRMZM2G072121 for the expansin precursor protein) for MEgreen in the division zone of the maize roots. Similarly, relatively high mRNA levels of GRMZM2G072121 (lipid-transfer protein) and GRMZM2G135536 (cytochrome P450 protein) were observed in primary root tissues at 6 DAS (Sekhon et al. 2011). Wei and Zhong (2014) have reported that four Zea mays lipid-transfer proteins (ZmLTPs: ZmLTP2.6, ZmLTP2.7, ZmLTPd1, and ZmLTPd4) may participate in root development. In addition, a gain-of-function mutant for Arabidopsis LTP5 encoding the lipid transfer protein domain exhibited primary root growth defects (Chae et al. 2010). Mutation of CYP82C2, an Arabidopsis cytochrome P450 protein, led to inhibited root growth and resulted in altered responses to JA. These findings support the hypothesis for the functional relationship of MEgreen associated with root emergence and formation during initial plant growth. In this study, we successfully detected the tissue-specific genes that participate in root emergence and formation at the initial root growth stage by using network-based transcriptome analysis. Our results might help in understanding the complex and distinct regulations of genes during initial root development. Acknowledgements  This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No. NRF2015R1D1A4A01015702) as well as a grant from the Korea Institute of Planning and Evaluation for Technology in Food, Agriculture, and Forestry (IPET) through Golden Seed Project, funded by Ministry of Agriculture, Food and Rural Affairs (MAFRA) (213009-05-2-SB710), and the “Cooperative Research Program for Agriculture Science & Technology Development” (Project No. PJ012649022018) of Rural Development Administration, Republic of Korea.

Compliance with ethical standards  Conflict of interest  The authors declare that they have no conflict of interest.

Ethical standards  This article does not contain any studies with human participants or animals performed by any of the authors.

References Ariel F, Diet A, Verdenaud M, Gruber V, Frugier F, Chan R, Crespi M (2010) Environmental regulation of lateral root emergence in Medicago truncatula requires the HD-Zip I transcription factor ­HB1. Plant Cell 22:2171–2183 Baâtour O et al (2012) Effect of growth stages on phenolics content and antioxidant activities of shoots in sweet marjoram (Origanum majorana L.) varieties under salt stress. Afr J Biotechnol 11:16486–16493 Bailey TL, Machanick P (2012) Inferring direct DNA binding from ChIP-sEq. Nucleic Acids Res 40:e128 Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120 Cai HG et al (2012) Mapping QTLs for root system architecture of maize (Zea mays L.) in the field at different developmental stages. Theor Appl Genet 125:1313–1324 Chae K, Gonong BJ, Kim SC, Kieslich CA, Morikis D, Balasubramanian S, Lord EM (2010) A multifaceted study of stigma/style cysteine-rich adhesin (SCA)-like Arabidopsis lipid transfer proteins (LTPs) suggests diversified roles for these LTPs in plant growth and reproduction. J Exp Bot 61:4277–4290 Coneva V et al (2014) Metabolic and co-expression network-based analyses associated with nitrate response in rice. Bmc Genomics 15:1056 Fan K et al (2014) Molecular evolution and expansion analysis of the NAC transcription factor in Zea mays. PLoS ONE 9:e111837 Gautier L, Cope L, Bolstad BM, Irizarry RA (2004) affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 20:307–315 Guingo E, Hébert Y, Charcosset A (1998) Genetic analysis of root traits in maize. Agronomie 18:225–235 Guo DS et al (2015) The WRKY transcription factor WRKY71/EXB1 controls shoot branching by transcriptionally regulating RAX genes in Arabidopsis. Plant Cell 27:3112–3127 He XJ et al (2016) Comparative RNA-seq analysis reveals that regulatory network of maize root development controls the expression of genes in response to N stress. PLoS ONE 11:e0151697 Huang DW, Sherman BT, Lempicki RA (2009) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4:44–57 Kim SG, Kim SY, Park CM (2007) A membrane-associated NAC transcription factor regulates salt-responsive flowering via FLOWERING LOCUS T in Arabidopsis. Planta 226:647–654 Lan P, Li WF, Lin WD, Santi S, Schmidt W (2013) Mapping gene activity of Arabidopsis root hairs. Genome Biol 14:R67 Langfelder P, Horvath S (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinform 9:559 Lin WD, Liao YY, Yang TJW, Pan CY, Buckhout TJ, Schmidt W (2011) Coexpression-based clustering of Arabidopsis root genes predicts functional modules in early phosphate deficiency signaling. Plant Physiol 155:1383–1402 Liu J, Li J, Chen F, Zhang F, Ren T, Zhuang Z, Mi G (2008) Mapping QTLs for root traits under different nitrate levels at the seedling stage in maize (Zea mays L.). Plant Soil 305:253–265 Loudet O, Gaudon V, Trubuil A, Daniel-Vedele F (2005) Quantitative trait loci controlling root growth and architecture in Arabidopsis thaliana confirmed by heterogeneous inbred family. Theor Appl Genet 110:742–753

13

Makkena S, Lamb RS (2013) The bHLH transcription factor SPATULA regulates root growth by controlling the size of the root meristem. Bmc Plant Biol 13:1 Mano Y, Muraki M, Fujimori M, Takamizo T, Kindiger B (2005) Identification of QTL controlling adventitious root formation during flooding conditions in teosinte (Zea mays ssp. huehuetenangensis) seedlings. Euphytica 142:33–42 Manoli A, Sturaro A, Trevisan S, Quaggiotti S, Nonis A (2012) Evaluation of candidate reference genes for qPCR in maize. J Plant Physiol 169:807–815 Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Dalla Favera R, Califano A (2006) ARACNE: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinform 7:S7 O’Keeffe K (2009) Maize growth & development. NSW Department of Primary Industries, State of New South Wales Pace J, Gardner C, Romay C, Ganapathysubramanian B, Lubberstedt T (2015) Genome-wide association analysis of seedling root development in maize (Zea mays L.). BMC Genomics 16:47 Peng H et al (2015) Transcriptomic changes during maize roots development responsive to Cadmium (Cd) pollution using comparative RNAseq-based approach. Biochem Biophys Res Commun 464:1040–1047 Raun AL, Borum J, Sand-Jensen K (2010) Influence of sediment organic enrichment and water alkalinity on growth of aquatic isoetid and elodeid plants. Freshwater Biol 55:1891–1904 Ravid U, Ikan R, Sachs RM (1975) Structures related to jasmonic acid and their effect on lettuce seedling growth. J Agr Food Chem 23:835–838 Russell RS (1978) Plant root systems: their function and interaction with the soil. Soil Sci 125:272 Ryan PR, Kochian LV (1993) Interaction between aluminum toxicity and calcium uptake at the root apex in near-isogenic lines of wheat (Triticum aestivum L.) differing in aluminum tolerance. Plant Physiol 102:975–982 Sack FD (1991) Plant gravity sensing. Int Rev Cytol 127:193–252 Salazar-Henao JE, Lin WD, Schmidt W (2016) Discriminative gene co-expression network analysis uncovers novel modules involved in the formation of phosphate deficiency-induced root hairs in Arabidopsis. Sci Rep 6:26820 Schiefelbein JW, Benfey PN (1991) The development of plant roots: new approaches to underground problems. Plant Cell 3:1147–1154

13

Genes & Genomics Sekhon RS, Lin HN, Childs KL, Hansey CN, Buell CR, de Leon N, Kaeppler SM (2011) Genome-wide atlas of transcription during maize development. Plant J 66:553–563 Shannon P et al (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13:2498–2504 Staswick PE, Su W, Howell SH (1992) Methyl jasmonate inhibition of root growth and induction of a leaf protein are decreased in an Arabidopsis thaliana mutant. Proc Natl Acad Sci USA 89:6837–6840 Thimm O et al (2004) MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J 37:914–939 Tomic S, Gabdoulline RR, Kojic-Prodic B, Wade RC (1998) Classification of auxin related compounds based on similarity of their interaction fields: extension to a new set of compounds. Internet J Chem 1:26 Trapnell C, Pachter L, Salzberg SL (2009) TopHat: discovering splice junctions with RNA-seq. Bioinformatics 25:1105–1111 Trapnell C et al (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7:562–578 Usadel B et al (2006) PageMan: an interactive ontology tool to generate, display, and annotate overview graphs for profiling experiments. BMC Bioinform 7:535 Wei KF, Zhong XJ (2014) Non-specific lipid transfer proteins in maize. BMC Plant Biol 14:281 Wei HR, Yordanov YS, Georgieva T, Li X, Busov V (2013) Nitrogen deprivation promotes Populus root growth through global transcriptome reprogramming and activation of hierarchical genetic networks. New Phytol 200:483–497 Yilmaz A, Nishiyama MY, Fuentes BG, Souza GM, Janies D, Gray J, Grotewold E (2009) GRASSIUS: a platform for comparative regulatory genomics across the grasses. Plant Physiol 149:171–180 Zimmermann P, Hirsch-Hoffmann M, Hennig L, Gruissem W (2004) GENEVESTIGATOR. Arabidopsis microarray database and analysis toolbox. Plant Physiol 136:2621–2632 Zou ZW, Ishida M, Li F, Kakizaki T, Suzuki S, Kitashiba H, Nishio T (2013) QTL analysis using SNP markers developed by nextgeneration sequencing for identification of candidate genes controlling 4-methylthio-3-butenyl glucosinolate contents in roots of radish, Raphanus sativus L. PLoS ONE 8:e53541