Although both ammonium and nitrate can be used by plants, nitrate is the major form of nitrogen in many soils (Chrispeels et al., 1999). We focused initially on 15 NRTs and two other genes that encode channels that transport nitrate (Figure 1), each of which has some experimental evidence of its function (Wang Y. -Y. et al., 2012). We first explored differential expression, since those genes are believed to be regulated to respond to certain signals and hence they may play a role during the studied biological process (Tarca et al., 2006). Among the 371 published Arabidopsis expression series datasets we collected from GEO, 50 datasets are root-specific and 49 datasets are leaf-specific (see Table S2). For each leaf- and root-specific dataset, we identified the DEG and tallied the number of experiments in which each gene was differentially expressed (see Section Materials and Methods). Comparing the differential expression events between roots and leaves was suggestive of the function of some genes. For example, NPF2.7 (formerly NAXT1) was differentially expressed in more than 30% of root-specific datasets but in only about 12% of leaf-specific datasets (Figure 1), suggesting a context-specific function of NPF2.7 in roots. Previous studies have demonstrated that NPF2.7 is involved in the excretion of nitrate from roots of Arabidopsis (Segonzac et al., 2007). Also, our results show that NPF2.13 (formerly NRT1.7) was differentially expressed in only about 20% of roots but in about 40% of leaves (Figure 1). It has been reported that NPF2.13 is involved in translocation of nitrate from old leaves to young leaves (Fan et al., 2009). Furthermore, NRT2.1 and NRT2.2 are differentially expressed in roots twice as much as in leaves (Figure 1), which corresponds to their roles in the uptake of nitrate from soil (Wang Y. -Y. et al., 2012). We must caution that there are caveats to using this approach as an indicator of functionality. For example, NRT2.4 has much higher differential expression in roots than in leaves, which is consistent with its function in root nitrate uptake (Kiba et al., 2012). However, NRT2.4 has a second role, relating to the loading of nitrate into the phloem in shoots, which would have been missed by the differential expression approach alone. These examples suggest that comparing the relative level of responsiveness of genes in particular tissues is one approach that could be combined with other approaches to focus functional genomics studies of gene networks, especially for large gene families like NRT (e.g., cytochrome P450s, glycosyltransferases, glycoside hydrolases, etc.). Additionally, the measure of plasticity in expression provided by this analysis is suggestive of the degree to which a gene's function within a tissue is dependent on context and conditions.

When building a coexpression network, Pearson Correlation Coefficient (PCC) is often used to measure the weight of correlation between two expression profiles. The challenge in selecting a cutoff to define what elements to include is that the minimum value of PCC that is significantly different from zero (i.e., no correlation) heavily depends on the sample size. For example, at the significance level 0.05, the minimum PCC is 0.6 when the sample size is 10, and 0.2 when the sample size is 100. Generally speaking, there is no standardization of the cutoff amongst studies, and it may vary dramatically (Jordan et al., 2004; Van Noort et al., 2004; Wang S. et al., 2012). Although the p-value of correlation can be used as cutoff between datasets of different sample size (Ponomarenko et al., 2013), it is tricky to calculate an average value using p-value. Since coexpression networks are intended to make complex systems understandable to the human intellect, others have included only a small number of most highly correlated genes (e.g., the top 20 genes or top 0.1%) (Kim et al., 2001; Bergmann et al., 2004). We utilized a similar strategy to focus attention on the most highly correlated network members (See Section Materials and Methods). A network was constructed using the top 20 coexpressed partners for our 17 focal genes (i.e., 15 NRT and 2 channels) (Figure 2). The weight (i.e., PCC) for each GEO dataset was calculated independently and the average value of all datasets was used to measure the strength of coexpression between a pair of genes. Unlike using a combined meta-dataset where the transient/context-specific signals may be swamped, those relationships are more likely to be captured by our method (Usadel et al., 2009).

Many of the 17 N transporters were coexpressed with several other N transporters, showing the potential functional association among those genes. In order to test whether this coexpression network makes biological sense, GO enrichment analysis was performed to detect over-represented functional categories after removing all of the 17 genes from the network (Figure 2). Interestingly, “water channel activity” was over-represented (p = 6.4 × 10⁻⁸), which indicates that transport of water and transport of nitrate may be under coordinated transcriptional control. Other than those 17 genes, there are various transporters, metabolic enzymes, and transcriptional regulators in this network (Figure 2), some of which appear likely to have a relationship with the nitrogen transport system, based on their known functions. For example, NPF6.3 (formerly NRT1.1) is highly coexpressed with H⁺-ATPase 2, AT4G30190 (Table 1 and Table S3). NPF6.3 is a nitrate/proton symporter, requiring a proton gradient for the uptake of nitrate from the soil (Parker and Newstead, 2014). The H⁺-ATPases that maintain the proton gradient comprise a large superfamily (Axelsen and Palmgren, 2001; Palmgren, 2001), but our coexpression analysis suggests that H⁺-ATPase 2, specifically, may contribute to the proton gradient needed for the transport of nitrate from the soil into the root by NPF6.3. Furthermore, H⁺-ATPase 2 is strongly expressed in the root pericycle, cortex, epidermis, and root cap according to the Arabidopsis eFP browser, particularly after nitrate addition (Figure S1A; Winter et al., 2007).

We further included 171 genes potentially involved in nitrogen transport in a similar analysis as above. A network of top 20 coexpressed partners for each nitrogen transporter can be visualized in Figure S2. In total, 2047 other genes are in this network, many of which are connected with more than one nitrogen transporter (Table S4). Interestingly, other transporter genes are enriched among those 2047 genes, such as genes from the GO categories “ion transport” and “carbohydrate transport” (Table S5), indicating those biological processes might be regulated similarly in Arabidopsis (Koprivova et al., 2000; Scheible et al., 2004). These other transporters could possibly be involved in transport of counter-ions to help maintain charge balance across membranes during sustained NO3- transport, or might be involved in the uptake or homeostasis of other essential nutrients that would be needed for growth and development at the same time as N uptake. Multiple other GO biological process categories were also significantly over-represented in the network (Table S5), such as those relating to phenolics. It is well documented that phenolic compound biosynthesis is upregulated when N is limited relative to C (Scheible et al., 2004; Cross et al., 2006). Our analysis suggests that there may be coordinated regulation of N transporter genes and phenylpropanoid biosynthetic genes. Additionally, there were various carbohydrate (C) metabolism and transport categories that were over-represented in the N transporter network (Table S5). This is consistent with previous evidence for extensive coordination to balance C and N metabolism (Palenchar et al., 2004) but may also indicate the need for increased carbohydrates in tissues where N uptake is strong to provide energy to maintain the proton gradient needed for N uptake, and to provide energy and organic building blocks for the lateral root proliferation that is common in high N regions of soil (Hodge, 2004). Finally, the network also included responses to numerous stimuli, such as water deficit, abscisic acid, and wounding. These “response” categories represent a rich resource for hypothesis generation, as they may reflect the importance of coordinating N utilization with other aspects of plant physiology in response to different environmental conditions. For example, N uptake may need to be altered if water uptake declines during drought, since N delivery to the shoot requires transport with water through the xylem. In addition to these over-arching insights, similar networks that focus on particular aspects of N transport (e.g., N export during leaf senescence) may be useful to identify a more focused set of processes that are associated with particular aspects of N utilization.

Increasing the coexpressed partners in a network beyond the top 20, as above, may be meaningful but there is a risk of increasing the false-positive rate. One strategy to detect those broader relationships when individual gene relationships are relatively weak is to compute the correlation between a gene and meaningful pathways, such as GO Biological Processes (Huang et al., 2006; Tegge et al., 2012; Bateman et al., 2014). Since a pathway is a pre-defined group of genes, taking all those genes into account may boost the power to detect the real signal (Lee et al., 2011). Furthermore, the presence or absence of specific genes or networks may be tissue- or cell-type dependent (Anglani et al., 2014). We calculated the correlation between expression of the 17 N transporter genes and GO Biological Process pathways within each GEO dataset for samples from the same tissue type, and used the top 10 correlations to construct a network of N transporter-pathways that displays the tissue-specificity of each edge (See Section Materials and Methods; Figure 3). Each connection between a N transporter and a GO pathway represents a statistically significant correlation in a tissue type, which is represented by the color of the edge. Some N transporters are connected by network edges of a single tissue type, such as NPF2.12 (NRT1.6) and NRT2.2, while others are connected by network edges of multiple tissue types, such as NPF6.2 (NRT1.4) and NPF4.6 (NRT1.2) (Figure 3).

These correlations may provide hints as to the function of uncharacterized N transporter genes or additional functions of previously characterized genes. NPF2.12 (formerly NRT1.6) is connected with several pathways in our N transporter-pathway network and all those relationships are based on seed-specific datasets (Figure 3). This is consistent with previous evidence, which suggests that NPF2.12 is involved in the delivery of nitrate from the maternal plant to the developing embryo, particularly the transfer of nitrate from the vascular tissue into the seed (Almagro et al., 2008). Knockout of NPF2.12 has profound impacts such as reduced nitrate content of seeds, and substantially increased incidence of seed abortion. Our network suggests that, in addition to carpel development, NPF2.12 may also be strongly linked with anther and pollen development, vacuolar protein localization, and phenolic metabolism. In recent years, intact phenolic metabolism has been linked with proper pollen development and pollen fertilization of embryos (Matsuno et al., 2009; Fellenberg et al., 2012; Fellenberg and Vogt, 2015).

All of the edges connected to NRT2.2 are based on root-specific datasets, which is consistent with its role in the uptake of nitrate from soil (Li et al., 2007). One of the pathways connected to NRT2.2 is “specification of organ identity,” which might reflect the tight relationship between nitrate uptake and cellular differentiation or between nitrate consumption and root growth (Walch-Liu et al., 2006). NRT2.2 and NRT2.1 shared a strong link with “imidazole-containing compound metabolic process,” which includes multiple genes associated with histidine biosynthesis. Although the relevance of this co-linkage to histidine biosynthesis is not immediately clear, it is interesting that NRT2.2 and NRT2.1 are linked in our network since the two genes reportedly have some overlap of function in inducible high affinity nitrate uptake by roots (Li et al., 2007).

Table 2 shows the top correlated pathway for each of the 17 N transporters. Detailed information about all those tissue-specific correlations can be found in Table S6. We believe that data such as these will provide potential candidate genes and interesting hypotheses for further studies. For example, leaf-specific “sulfate (S) assimilation” is the best correlated pathway with NPF6.3 (formerly NRT1.1) and NPF7.3 (formerly NRT1.5). It is not surprising that N uptake and S assimilation genes correlate, since plant processes that require a lot of N also tend to need S, for example for the biosynthesis of cysteine, methionine, and several important cofactors (Koprivova et al., 2000). Similarly, multiple genes were best associated with, “photosynthesis, light harvesting,” including NPF4.6, NPF6.2 (formerly NRT1.2 and NRT1.4, respectively), NRT2.7, and CLCA, which appear as a cluster in the network (Figure 3). This may reflect the importance of tight co-regulation of N and C metabolism and also the fact that the biosynthesis of light harvesting proteins and pigments is highly dependent on the availability of N (Scheible et al., 2004).

Another example is NPF2.13 (formerly NRT1.7) which is coexpressed with “negative regulation of cell death” pathway in leaf. As reported previously, NPF2.13 plays a role in recycling of nitrogen within the plant (Fan et al., 2009) and cell death is a prominent event during senescence (Lim et al., 2007). It is probably crucial that cell death is slowed or delayed until most of the N is exported from the leaves through the phloem. Several other N transporter genes are also strongly correlated with the “negative regulation of cell death” pathway (NPF4.6/NRT1.2, NPF7.2/NRT1.8, NRT2.6, and NRT2.5), which might suggest that multiple N transporter genes are involved in N remobilization during senescence. NRT2.5 has been linked previously with N remobilization (Lezhneva et al., 2014). Alternatively, the “negative regulation of cell death” hub might indicate a general role of abundant N in delaying senescence. Indeed nitrogen status and senescence are known to be closely linked (Cooke et al., 2005; Diaz et al., 2008).

Some of the other GO Biological Process pathways appear to represent informative hubs. For example, the “nitrate transport,” and “response to nitrate,” groups are coexpressed with multiple genes in root tissues, including NPF6.3, NPF4.6, NPF2.9, NPF2.7, and SLAH3 (Figure 3). Several of these genes are known to be involved in nitrate uptake and redistribution in roots (Segonzac et al., 2007; Wang and Tsay, 2011; Glass and Kotur, 2013) and this association suggests that the others might play other roles in these processes. For example, SLAH3 functions in nitrate release from guard cells (Geiger et al., 2011), but based on the fact that SLAH3 is strongly expressed in the pericycle of roots (Figure S1B) combined with this coexpression relationship, one could hypothesize that SLAH3 might facilitate nitrate loading into the xylem or phloem in the roots via an apoplastic route. The role of NPF2.9 (formerly NRT1.9) was described as mediating nitrate distribution between shoot and root (Wang and Tsay, 2011). Based on expression of NPF2.9 in companion cells in roots and the relationship of NPF2.9 with NPF6.3 and NPF4.6 (formerly NRT1.1 and 1.2, respectively) in our coexpression network, perhaps NPF2.9 might have a more direct link with nitrate uptake, such as delivery of nitrate from the maturation zone of the root, where much of the water and nitrate uptake occurs, toward the developing root tip via the phloem. By identifying clusters or hubs such as this, our analysis provides guidance for further experimentation. For example, in order to understand nitrate uptake, we need to understand the functions of these genes and how these functions are integrated as a system.

Three hundred and seventy one expression series datasets based on platform GPL198 were collected from GEO (Table S7). Each dataset contains at least 12 samples. Three hundred and twenty datasets which contain CEL files were used in our analysis. Robust Multiarray Average (RMA) was used to normalize the microarray data for each dataset (Irizarry et al., 2003). The IDs of probesets were converted into gene locus ID based on the annotation file for GPL198. The replicate group and tissue types were manually curated.

For each expression dataset, we applied ANOVA to identify gene expression which has larger variation between two replicate groups than within a replicate group. Only those with false discovery rate (FDR) < 0.001 were considered as DEG.

We used the following equation to measure the strength of coexpression between a nitrogen transporter and another gene on the array:

where r_k is the coexpression weight (i.e., PCC) between gene i and nitrogen transporter x in the GEO dataset n. R_{x, i} is the average value of 320 weights between gene i and nitrogen transporter x. We used the following equation to measure the strength of tissue-specific coexpression between a nitrogen transporter and another gene on the array:

where T is the subset of GEO dataset of a specific tissue type. Rx,i′ is the strength of tissue-specific coexpression between a nitrogen transporter x and gene i. r_k in Equation (2) represents the weights from a specific tissue type, T. We used the “biological process” in the GeneOntology system to define a pathway and the following equation was used to measure the tissue-specific coexpression between a nitrogen transporter x and a pathway:

where m is the number of genes in a pathway p. Rx,k′ is the tissue-specific coexpression between nitrogen transporter x and another gene k. And k is a gene in the pathway p. When we constructed the coexpression network between nitrogen transporters and pathways, only the datasets where a nitrogen transporter is differentially expressed were used. In order to determine the statistical significance of Rx,p′, the genes of a pathway were replaced by randomly selected genes in the genome and the Rx,p′ was calculated. We repeated this process 100 times for each pathway. A empirical p < 0.01 was assigned if none of the resulted Rx,p′ is higher than the real Rx,p′. For more detail, see Supplemental Note.

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.