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 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 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¹ (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²). 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³ were used for mapping, and the reference genome was Arabidopsis_thaliana.TAIR10.28.dna.toplevel.fa. Genome annotation was based on Ensembl⁴ and KEGG (Kyoto Encyclopedia of Genes and Genomes⁵). Gene expression analysis was performed using HTSeq⁶ and DESeq⁷ at different expression levels (Supplementary Table S6). 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 S4 and Supplementary Tables S7, S8).

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 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₂SO₄ and H₂O₂. Nitrate concentration was measured at OD₂₁₀. KNO₃ was used for drawing the standard curve. The formula was A = 7.5857C_NO3- – 0.0466. r = 0.9970.

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.

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

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

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 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 2B,C and Supplementary 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.

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 down-regulated 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).

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

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

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

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