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., 2011Li et al., , 2016El-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., 2011Kuai 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.

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-weekold 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 S1). 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 S6). Up-and downregulation 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 S4 and Supplementary Tables S7, S8).

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 S9, and were assessed via a standard curve. The amplification and melting curves are shown in Supplementary Figure S5. The values were reported using the 2 C T 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 H 2 SO 4 and H 2 O 2 . Nitrate concentration was measured at OD 210 . KNO 3 was used for drawing the standard curve. The formula was A = 7.5857C NO3 − -0.0466. r = 0.9970.

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 S1), with the GC distribution close to theoretical distribution (Supplementary Figure S2). 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 S2). 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 S3). A percentage of 90% of the reads were mapped onto the Arabidopsis genome, most of which were uniquely mapped (Supplementary Table S3). Among those genome-mapped reads, about 90% were mapped to gene region, and more than 99% were mapped onto exons (Supplementary Table S4). In total, 20,660 sequenced genes were recovered and were distributed in five nuclear chromosomes (Supplementary Table S5).

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 S6). Under high growth density, 98 genes were up-regulated and 107 genes were down-regulated in comparison with the LD conditions (Supplementary Figure  S4A and Supplementary Table S6).
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 S4B and Supplementary Table S7). 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  S4B and Supplementary Table S7). 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 S4B and Supplementary Table S7).
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 S8). 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 S4C and Supplementary  Table S8). In addition, the DNA replication and repair pathway in genetic information processing was influenced by high growth density with P < 0.05 (Supplementary  Table S8). Generally, metabolism was influenced more than genetic functioning based on the significantly smaller P-values (Supplementary Figure S4C and Supplementary Table S8). The other four categories did not exhibit much change between the high and low growth densities (Supplementary Table S8).

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 1A). After that, they grew upward to avoid the shade. The rosette leaf angle thus became smaller compared with the low-density plants (Figures 1B,C). However, the rosette leaf diameter was similar between the two treatments ( Figure 1D).
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 1). However, flowering time was similar between the high and low densities ( Table 1). 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 S6). qPCR validation demonstrated that FD was induced about 2-fold when the plants were grown at high density (Figure 2A). 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  Table S6).
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 S6). However, qPCR analysis indicated that they were induced about 1.34, 2.30, and 1.68-fold, respectively (Figure 3), indicating that shade avoidance was slightly induced.
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 S6). qPCR validation demonstrated that all these genes were downregulated at least two-fold under high density conditions in comparison to LD conditions (Figures 4A-E). 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 4F).

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 S6).

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 S6). GO enrichment analysis suggested that density obviously influences metabolic processes (Supplementary Figure S4B and Supplementary Table S7), which is consistent with the findings of the KEGG analysis whereby three metabolic pathways were significantly affected (Supplementary Figure S4C and Supplementary Table S8). Since density constitutes an environmental factor, stimulus-and stress-related genes were greatly affected (Supplementary Figure S4B and Supplementary Frontiers in Plant Science | www.frontiersin.org FIGURE 2 | 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.
FIGURE 3 | 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. Table S7), 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 S4B and Supplementary Table S7). 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 S6), 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 1A). 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 3). Another SAS response was the reduction in rosette leaves at bolting time ( Table 1). 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 1). Upon evaluating the RNA-Seq data, we found that only one flowering-related gene, FD, was up-regulated (Supplementary Table S6). qPCR validation demonstrated that the induction of FD was about 2-fold (Figure 2A). We speculated that this phenotype was caused by the low induction of FD. Since FIGURE 4 | 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.
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 2A). 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 2C and Supplementary  Table S6). 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 2B). 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 4). 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 downregulated. 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 4F). 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 S6). 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 mostused 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 downregulation 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 S6). In addition, the RNA-Seq samples of the different densities were highly correlated (Supplementary Figure S3), 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.

FUNDING
This work is supported by the Natural Science Foundation of Hebei Province (C2015209045), Science Foundation of Hebei Province Department of Education (QN2014003), Science and Technology Support Program of Tangshan City (15140202a), and the Doctoral Scientific Research Foundation of North China University of Science and Technology to DG, China National Science Foundation (3151333) to JW, and the Hebei 100 Talented Scholars project, and Tangshan Key Laboratory Project to XW.

ACKNOWLEDGMENTS
The authors thank Prof. Ning Li in HKUST for Col-0 seeds donation. They thank Personal Biotechnology for RNA sequencing and qPCR performance. They also thank LetPub (http://www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017.02001/ full#supplementary-material FIGURE S1 | Clean reads quality analysis using FastQC. Q > 28 is considered to be high quality and is labeled with a green background.