RNA-Seq Profiling Shows Divergent Gene Expression Patterns in Arabidopsis Grown under Different Densities.

Plants growing under high-density (HD) conditions experience increased competition for water, nutrients, and light, possibly leading to changes in size, biomass, morphology, and productivity. However, no research has focused on the relationship between whole-genome expression patterns and growth density. Here, we performed whole-genome RNA sequencing to examine the gene expression patterns in Arabidopsis grown under low and high densities. Of the 20,660 detected genes, the expression levels of 98 were enhanced and 107 were repressed under HD growth. Further analysis revealed that changes in density influenced metabolism- and stimulus-related genes the most. Furthermore, HD growth led to a shade avoidance phenotype, represented by upward growth and a reduction in rosette leaves. Moreover, a cluster of glutaredoxin genes, GRXS3, 4, 5, 7, and 8, were significantly down-regulated under high density, suggesting that high density affects plant growth mainly by nitrate limitation.

Introduction

Most plant studies have been conducted using individual potted plants under controlled conditions (Hecht et al., 2016), and much physiological and molecular data have been derived from this. Sowing density is important as plants grown under high-density (HD) conditions compete with each other for water, nutrition, and light, which often leads to changes in plant size, biomass, morphology, and productivity (Hecht et al., 2016). Research about plant density refers to productivity, organ development, nutrition absorption, water heterogeneity, shade avoidance, and flowering (Lemaire et al., 2005; Hagiwara et al., 2010; Li et al., 2011, 2016; El-Zaeddi et al., 2016; Roig-Villanova and Martinez-Garcia, 2016; Song et al., 2016).

Studies on plant density and productivity have mainly focused on crops, vegetables, and medical plants, including wheat, maize, potato, Brassica, and Salvia miltiorrhiza. People measure the organ size and biomass to determine the optimal sowing density (Li et al., 2011, 2016; Kuai et al., 2015; Liu et al., 2016; Song et al., 2016; Zheng et al., 2016). Besides productivity, plant density also affects organ development (Hecht et al., 2016; Song et al., 2016). The roots of spring barley increased with higher sowing density, particularly in the topsoil. In contrast, the root mass decreased while the stem mass increased with higher sowing density (Hecht et al., 2016). The effects of plant density on maize canopy structure indicated that both lamina width and internode diameter were reduced under high density as a result of a lower growth rate (Song et al., 2016).

Plants compete for nutrition, water and light under HD conditions. Shoot nitrogen (N) dynamics analysis in lucerne showed that N accumulation was less rapid when plant density was high, which was caused by a decrease in leaf biomass (Lemaire et al., 2005). Additionally, plant density influences the volatile compositions of dill, suggesting differences in nutrient absorption at various densities (El-Zaeddi et al., 2016). The effects of water supply heterogeneity on Perilla frutescens were greater at higher density than at lower density, indicating that competition for water makes plants more sensitive to water heterogeneity (Hagiwara et al., 2010). To respond the light competition, plants employ two opposing strategies: shade tolerance and shade avoidance. Shade avoidance is associated with a set of responses known as the shade avoidance syndrome (SAS), including stem-like organ elongation, apical dominance, flowering acceleration, branch reduction, and decreasing leaf expansion and yield. Phytochromes, parotid isoelectric focusing (PIF) proteins, HD-ZIP II transcription factors, and auxins were found to be involved in shade avoidance (Morelli and Ruberti, 2002; Ruberti et al., 2012; Gonzalez-Grandio et al., 2013; Roig-Villanova and Martinez-Garcia, 2016).

Flowering is a pivotal event for plants and represents the switch from vegetative to reproductive development. Botto and Coluccio (2007) analyzed the effects of plant density on flowering time in Arabidopsis using recombinant inbred lines, and found a variation in flowering time across seasonal and density environments.

Since plant density affects many traits, we want to characterize gene expression under different growth densities. Here, we use RNA sequencing to examine the whole-genome expression patterns of Arabidopsis under high and low growth densities to identify density-regulated genes.

Materials and Methods

Plant Materials and Growth Conditions

Arabidopsis Col-0 was provided by Prof. Ning Li at the Hong Kong University of Science and Technology. The soil used in this study was Klasmann-Deilmann 876. Seeds were first imbibed at 4°C for 3 days before being sowed in the soil. The plants were grown at 22 ± 1°C under continuous white light at an intensity of 100 μmol photons m-2 s-1. The pots used in this study were 32 × 48 × 9 (depth) cm. The distance between the seeds sown at low density (LD) was 10 cm, while the distance between the seeds in the HD treatment was 2 cm. Three-week-old leaf samples were collected from directly under the white light for angle and diameter measurements, RNA sequencing, and qPCR examination. Plants situated on the edge of the pots were excluded from all experiments. For RNA sequencing, the plants were randomly selected and the rosette leaves of three individual plants were separately collected. More than one leaf was used for the RNA extraction. Leaf samples from different plants were not pooled.

RNA Extraction and cDNA Preparation

RNA extraction and cDNA preparation kits were provided by Personal Biotechnology Co., Ltd., Shanghai, China. For RNA-Seq analysis, total RNA was extracted using an RNAout kit (TIANDZ, CAT#: 71203). For qPCR, total RNA was extracted using TRIZOL (Invitrogen). mRNA was purified and cDNA was prepared using the Truseq Stranded mRNA LT Sample Prep Kit (Illumina).

RNA Sequencing

RNA sequencing services were provided by Personal Biotechnology Co., Ltd., Shanghai, China. The constructed cDNA libraries were examined using an Agilent High Sensitivity DNA Kit, and the average fragment length was 250 bp. The libraries were sequenced using Illumina NextSeq500 to generate paired-end reads of 150 bp. The sequenced raw reads were processed to obtain clean reads using the following strategy: (i) remove reads with adaptor contamination; (ii) remove low quality reads whose average quality was less than Q20; and (iii) remove the reads with a final length of less than 50 bp. The quality of the clean reads was assessed using FastQC[1] (Supplementary Figure ). Three replicates of each treatment were sequenced and the raw data were uploaded to NCBI and deposited in the sequence read archive (SRA[2]). The biosample accession isSAMN07267285 and the corresponding address is: https://www.ncbi.nlm.nih.gov/sra/?term=SAMN07267285.

High-quality clean reads were used for further analysis. Bowtie2 and TopHat2[3] were used for mapping, and the reference genome was Arabidopsis_thaliana.TAIR10.28.dna.toplevel.fa. Genome annotation was based on Ensembl[4] and KEGG (Kyoto Encyclopedia of Genes and Genomes[5]). Gene expression analysis was performed using HTSeq[6] and DESeq[7] at different expression levels (Supplementary Table ). Up- and down-regulation of genes was considered to have occurred when the P-value was <0.05 and the absolute fold changes were ≥2.0. Gene ontology (GO) and KEGG orthology (KO) analysis of the genes with differential expression was also performed (Supplementary Figure and Supplementary Tables ).

qPCR Analysis

The qPCR conditions were set as follows: 95°C for 5 min, followed by 40 cycles of 95°C at 15 s and 60°C at 30 s. ACT2 was used as an internal control. Primers used in the qPCR are shown in Supplementary Table , and were assessed via a standard curve. The amplification and melting curves are shown in Supplementary Figure . The values were reported using the 2ΔΔC method. A P-value < 0.05 and absolute fold change ≥2.0 signified differential expression.

Nitrate Content Determination

Nitrate determination was performed according to an approach in a previous report (Lv et al., 2004). Three-week rosette leaves were collected, dried, and digested using the Kjeldathl method with H2SO4 and H2O2. Nitrate concentration was measured at OD210. KNO3 was used for drawing the standard curve. The formula was A = 7.5857C – 0.0466. r = 0.9970.

Gene Name and Corresponding Locus in This Study

GRXS3: At4g15700; GRXS4: At4g15680; GRXS5: At4g15690; GRXS7: At4g15670; GRXS8: At4g15660; FD: At4g35900; FT: At1g65480; FUL: At5g60910; AtHB-2: At4g16780; HFR1: At1g02340; PIL1: At2g46970; ACT2: At3g18780.

Results

RNA Sequencing of Arabidopsis under Different Growth Densities

RNA sequencing was used to examine the gene expression patterns in Arabidopsis under low and high growth densities. For each treatment, three replicates were sequenced using Illumina NextSeq500. Approximately 30 million reads of raw tags were obtained for each sample (Supplementary Table ), with the GC distribution close to theoretical distribution (Supplementary Figure ). After filtering with stringent criteria (see section Materials and Methods), we finally obtained 31.26, 32.18, and 28.70 million reads for the high growth density samples, and 31.98, 30.10, and 31.56 million reads for the low growth density samples (Supplementary Table ). All raw data were uploaded onto NCBI SRA and the biosample accession is SAMN07267285. The Pearson’s correlation coefficients of the three replicates for each density were larger than 80%, suggesting appreciable correlation between them (Supplementary Figure ). A percentage of 90% of the reads were mapped onto the Arabidopsis genome, most of which were uniquely mapped (Supplementary Table ). Among those genome-mapped reads, about 90% were mapped to gene region, and more than 99% were mapped onto exons (Supplementary Table ). In total, 20,660 sequenced genes were recovered and were distributed in five nuclear chromosomes (Supplementary Table ).

Gene Expression of Arabidopsis under Different Growth Densities

Among the 20,660 genes, the expression levels of 20,455 were not influenced by growth density, and only 205 constituted differentially expressed genes (DEGs) under the different growth densities (Supplementary Table ). Under high growth density, 98 genes were up-regulated and 107 genes were down-regulated in comparison with the LD conditions (Supplementary Figure and Supplementary Table ).

We classified the density-related genes by GO enrichment analysis. Biological processes, metabolic processes, cellular processes, death, stimulus, and stress obtained P-values < 0.01 and were thus likely related to changes in plant density (Supplementary Figure and Supplementary Table ). Cellular components, nucleus, endoplasmic reticulum, external encapsulating structure, and extracellular region were significantly related to changes in density, and the corresponding P-values for external encapsulating structure and extracellular region were <0.01 (Supplementary Figure and Supplementary Table ). With regards to molecular functions, the most changed DEGs were detected in the binding group, indicating significant change with P-values < 0.05 (Supplementary Figure and Supplementary Table ).

The KEGG enrichment analysis allowed us to divide all the genes into six categories, including metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases (Supplementary Table ). We discovered that growth density was the most influential metabolism category, among which amino acid metabolism, secondary metabolite biosynthesis, and xenobiotic biodegradation were significantly affected (Supplementary Figure and Supplementary Table ). In addition, the DNA replication and repair pathway in genetic information processing was influenced by high growth density with P < 0.05 (Supplementary Table ). Generally, metabolism was influenced more than genetic functioning based on the significantly smaller P-values (Supplementary Figure and Supplementary Table ). The other four categories did not exhibit much change between the high and low growth densities (Supplementary Table ).

High Density Causes a Shade Avoidance Response Phenotype

In our study, we found that plants grown under high density displayed a shade avoidance response, which was represented by a smaller rosette leaf angle and reduction in rosette leaves. The leaves of Col-0 began to reach the proximity of the other plants at about 2 weeks (Figure ). After that, they grew upward to avoid the shade. The rosette leaf angle thus became smaller compared with the low-density plants (Figures ). However, the rosette leaf diameter was similar between the two treatments (Figure ).

Characterization of Col-0 at high and low densities. (A) Photo of Col-0 at different stages and densities. Bars = 2 cm. (B) Side view of 3-week-old plants in different growth densities. Bars = 3 cm. (C,D) Rosette leaf angle (C) and rosette diameter (D) of 3-week-old plants. ∗∗Represents P-value < 0.01.

Another shade avoidance response is the reduction in rosette leaves at bolting time. Although there was variation between the different sets of experiments, the bolting time of the rosette leaves of the HD plants was 2 less than those grown at LD (Table ). However, flowering time was similar between the high and low densities (Table ). As rosette leaf number is an important indicator of flowering time (Abe et al., 2005; Wigge et al., 2005), we searched for flowering-related genes and found that only FD was up-regulated under high density (Supplementary Table ). qPCR validation demonstrated that FD was induced about 2-fold when the plants were grown at high density (Figure ). It was reported that the combination of FD with FT stimulates flowering, and the over-expression of FD induces FUL expression (Wigge et al., 2005; Amasino, 2010; Srikanth and Schmid, 2011). Here, we assessed the expression levels of FT and FUL under different growth densities and discovered that neither FT nor FUL showed any difference between high and low growth densities in both the RNA-Seq and qPCR analysis (Figures and Supplementary Table ).

qPCR analysis of some flowering-related genes at different densities. Three-week-old leaf samples under high and low growth densities were collected. qPCR was performed to examine the expression levels of (A–C) FD, FT, and FUL. ACT2 was used as an internal control. LD, low density; HD, high density. ∗∗Represents P-value < 0.01 and the fold changes > 2.

Since a shade avoidance response phenotype was detected, we further assessed any changes in SAS-related genes. We found that marker genes, such as ATHB-2, HFR1, and PIL1, did not show any obvious differences in the RNA-Seq analysis (Supplementary Table ). However, qPCR analysis indicated that they were induced about 1.34, 2.30, and 1.68-fold, respectively (Figure ), indicating that shade avoidance was slightly induced.

qPCR analysis of SAS-related genes at different densities. Three-week-old leaf samples grown under high and low densities were collected. qPCR was performed to examine the expression levels of (A–C) AtHB-2, HFR1, and PIL1. ACT2 was used as an internal control. LD, low density; HD, high density. ∗∗Represents P-value < 0.01 and the fold changes > 2.

High Density Repressed the GRXS3/4/5/7/8 Cluster Genes

The RNA sequencing results showed that a cluster of GRXS genes on chromosome 4, named GRXS3, 4, 5, 7, and 8, were significantly down-regulated under high density (Supplementary Table ). qPCR validation demonstrated that all these genes were down-regulated at least two-fold under high density conditions in comparison to LD conditions (Figures ). Since this GRXS cluster was previously reported to be closely associated with soil nitrates (Patterson et al., 2016; Walters and Escobar, 2016), the nitrate content of the rosette leaves was assessed. However, the nitrate content was similar between plants grown under high and low densities (Figure ).

qPCR analysis of GRXS genes and rosette leaf nitrate content at different densities. Three-week-old leaf samples under high and low growth densities were collected. (A–E) qPCR was performed to examine the expression levels of the GRXS3, 4, 5, 7, and 8 genes. ACT2 was used as an internal control. ∗∗Represents P-value < 0.01 and the fold change > 2. (F) Nitrate content in the rosette leaves was measured. LD, low density; HD, high density; mg/gdw, mg per gram dry weight.

High Density Does Not Induce an Abiotic Stress Response

Plants growing under high density compete with each other for water, light, and nutrients, which often leads to changes in size, biomass, morphology, and productivity (Hecht et al., 2016). Thus, extreme high density constitutes an abiotic stress. A previous study showed that abscisic acid (ABA) was responsible for plant defense against abiotic stresses, while ethylene (ET), salicylic acid (SA), and jasmonate (JA) mainly responded to biotic stresses (Verma et al., 2016). Here, we assessed the expression patterns of some ABA biosynthesis and response genes, such as ABA1, AAO3, NCED3, ABA3, ABF2, ABF3, and ABF4 (Verma et al., 2016). However, no significant changes were observed in the expression of all these genes between the high and low densities (Supplementary Table ).

Discussion

Growth Density Caused Variation in Global Gene Expression

We found that relatively few (1%, 205/20,660) of the considered Arabidopsis genes were differentially expressed under the different planting densities (Supplementary Table ). GO enrichment analysis suggested that density obviously influences metabolic processes (Supplementary Figure and Supplementary Table ), which is consistent with the findings of the KEGG analysis whereby three metabolic pathways were significantly affected (Supplementary Figure and Supplementary Table ). Since density constitutes an environmental factor, stimulus- and stress-related genes were greatly affected (Supplementary Figure and Supplementary Table ), and all P-values were < 0.01, implying interaction between these genes and the changing environment. By performing cellular component analysis, we found that genes in the extracellular region, particularly the external encapsulating structure, changed the most (P-values < 0.01, Supplementary Figure and Supplementary Table ). This group of genes is located in the primary cell wall, implying that density influences cell wall development.

High Density Led to a Shade Avoidance Phenotype

Plants exhibit shade avoidance when grown at high density as the proximity to other vegetation results in competition for light. However, most SAS-related genes did not show much change in this study (Supplementary Table ), although a shade avoidance phenotype was observed. This may be the result of the growth conditions set in this study. Previous studies reduced both the ratio and intensity of R/FR light to induce a shade avoidance response (Sessa et al., 2005; Ciolfi et al., 2013). This treatment ensures that the entire plant is under low light and cannot escape the shade, which often induces a quick and obvious shade avoidance response. In our study, continuous white light was used and the distance between the Arabidopsis seeds in the HD treatment was 2 cm. In the first 2 weeks, the seedlings did not crowd each other, which allowed each plant to get enough light. Therefore, no shade avoidance was induced at this time. When the plants became larger, they started touching at their proximity at about 2 weeks (Figure ). Shade avoidance was then induced and some response genes were up- or down- regulated. After that, the rosette leaves changed the growth direction to avoid the shade. When the rosette leaves were growing upward, 2 cm offered enough space for the plants to avoid the shade and obtain adequate light. Since the shade avoidance response is rapid and reversible (Morelli and Ruberti, 2002), the plants then stopped the shade avoidance response and those up- or down- regulated genes restored to normal levels. As we collected the samples at 3 weeks, the additional 1-week was sufficient for the restoration of the SAS response. This may explain why the induction of SAS-related genes was relatively low (∼2-fold) in this study (Figure ).

Another SAS response was the reduction in rosette leaves at bolting time (Table ). The rosette leaf number and bolting time are two indicators of flowering time (Pouteau et al., 2004; Abe et al., 2005; Wigge et al., 2005). In most cases, a reduction in rosette leaves is accompanied by earlier bolting (Pouteau et al., 2004). However, in our study there was no difference in bolting time, and only a reduction in rosette leaves was observed (Table ). Upon evaluating the RNA-Seq data, we found that only one flowering-related gene, FD, was up-regulated (Supplementary Table ). qPCR validation demonstrated that the induction of FD was about 2-fold (Figure ). We speculated that this phenotype was caused by the low induction of FD. Since SAS caused earlier flowering, we hypothesized that when the plants began to encroach on each other, shade avoidance occurred and FD was up-regulated. Therefore, the plants successfully avoided the shade and stopped the SAS response, resulting in a low induction of FD (∼2-fold).

It was also reported that the over-expression of FD by the 35S promoter caused obvious earlier flowering, representing by 5 fewer rosette leaves than the wild-type at bolting time (Wigge et al., 2005). Compared with the constitutive promoter 35S, the induction fold of FD was very low (∼2-fold) in our study (Figure ). This low induction is not enough to cause the obvious flowering phenotype. In our study, only 2 fewer rosette leaves were observed, and no earlier bolting was noted. In addition, this low FD induction was not enough to induce other related genes, such as FUL and FT. The overexpression of FD by the 35S promoter resulted in an enhanced expression level of FUL in the leaves (Wigge et al., 2005). FUL was only induced 1.2-fold in our results (Figure and Supplementary Table ). Although FD, in combination with FT, was considered to stimulate plant flowering (Abe et al., 2005), the expression level of FT in our study did not differ under the different densities (Figure ). Taken together, we suggest that the reduction in rosette leaves was caused by enhanced FD expression.

A similar phenotype was reported previously (Botto and Coluccio, 2007). For the Arabidopsis ‘Bay’ ecotype, bolting time was similar between the high and low densities, but three fewer rosette leaves at bolting time were observed at high density than at LD (Botto and Coluccio, 2007). In their study, 8 seeds per pot (8 cm in diameter) constituted the HD treatment, which is similar to our conditions (2 cm distance). We speculated that the flowering phenotype in their study was also caused by SAS and induction of the FD gene.

Plants under High Density Compete for Nitrate Absorption

In our study, we found that GRXS3, 4, 5, 7, and 8 were significantly repressed under high density (Figure ). Previous literature reported that the GRXS3, 4, 5, 7, and 8 genes belong to class III glutaredoxins and have high similarity in both RNA and protein sequences (Gutsche et al., 2015; Walters and Escobar, 2016). Silencing GRXS3 led to the reduced expression of GRXS4, 5, 7, and 8, and a longer primary root compared with the control group, suggesting that they acted as negative regulators of primary root growth (Patterson et al., 2016; Walters and Escobar, 2016). Therefore, we presumed that Arabidopsis grown at high density possessed longer roots compared with that at LD, since the expression of GRXS3, 4, 5, 7, and 8 was obviously down-regulated. At this stage we must note that we failed to retrieve root length data because it was difficult to separate the roots from the soil, especially in the HD treatment. However, a previous study reported that root length in barley increased with increased sowing density, especially in the topsoil (Hecht et al., 2016). Combined with our results, this increase might also be caused by the depression of GRXS genes.

GRXS genes are particularly sensitive to nitrate, but not ammonia. Nitrate treatment causes an obvious accumulation of GRXS genes and an inhibitory effect on primary roots. Silencing the GRXS3, 4, 5, 7, and 8 genes resulted in a longer primary root compared with the control group, implying that GRXS genes mediated the inhibitory effect of nitrate on primary root growth (Patterson et al., 2016; Walters and Escobar, 2016). When plant roots are grown in nitrate-rich regions of the soil, the accumulated nitrate induces glutaredoxins, which suppress primary root growth. Lateral roots are also induced by soil nitrate, which enhance the root absorption ability (Zhang and Forde, 1998; Patterson et al., 2016; Walters and Escobar, 2016). Considering our data, we hypothesized that plants at LD do not need to compete for nutrients, and the adequate nitrate would induce the GRXS genes and suppress primary root growth. In contrast, plants at high density compete with each other for soil nutrients, including nitrate, which leads to a surrounding nitrate deficiency. This deficiency will cause a reduction in GRXS genes, which results in a longer primary root. The longer root helps the plants to reach more nitrates in the soil. This explains why nitrate accumulation was similar between the different densities (Figure ). This mechanism helps plants to compete for nutrition and survive when growing under high density.

Based on our RNA-Seq analysis, we did not find significant expression changes in other nutrition-related genes between the two growth densities, including ammonia (Supplementary Table ). This result indicated that nitrate is the major limiting factor in agricultural productivity when plants are grown under high density. This is not surprising as nitrate is the most-used fertilizer in agriculture as a source of nitrogen compared with ammonium and urea (Noguero and Lacombe, 2016). In addition, some literature has reported that cabbage prefers nitrate to ammonium as a nitrogen source (Tian and Li, 2000; Zhang and Wei, 2002). Since both Arabidopsis and cabbage belong to Brassicaceae, nitrate should also be the primary nitrogen source of Arabidopsis. That explains the obvious down-regulation of glutaredoxins, which are especially sensitive to nitrate.

High Density Does Not Cause a Stress Effect When There Is Adequate Water, Nutrition, and Light

High density did not cause an abiotic stress response, and similar expression levels of ABA biosynthesis and response genes were observed (Supplementary Table ). In addition, the RNA-Seq samples of the different densities were highly correlated (Supplementary Figure ), indicating that growth density did not cause any serious global changes in gene expression. The possible reason for the phenotype is that although it was a HD treatment, the water, nutrition, and light were adequate. Therefore, we conclude that under conditions of adequate light, water, and nutrition, density itself does not significantly affect plant growth in Arabidopsis.

Author Contributions

XW conceived and led the research. DG and XS performed all the experiments and data analysis. XW, DG, and XS wrote the paper. MY, ZW, WG, LW, and JW participated in data discussion and manuscript revision.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Table 1

Bolting time and rosette leaf number of Col-0 under different growth conditions.

	Experiment 1			Experiment 2
	Bolting time	Rosette leaf	N	Bolting time	Rosette leaf	N
	(day)	number		(day)	number
Low density	26.5 ± 2.2	14.8 ± 1.6 (12–18)	28	36.8 ± 3.0	21.0 ± 2.1 (18–25)	25
High density	25.4 ± 1.7	12.8 ± 1.1 (11–15)	97	36.9 ± 1.5	18.7 ± 1.1 (17–21)	79