Top kth threshold (i.e. top 300 MR genes) is well suited for designing a biological experiment based on co-expressed genes using the guide-gene PubMed ID:https://www.ncbi.nlm.nih.gov/pubmed/10956013 approach [28]. The rationale of using this threshold in the present study was supported by examining the distribution of values within the top 100 HRR for each gene (Additional file 1: Table S4). The average HRR value was between 200 and 245, while the median HRR value was between 150 and 160; both of which were well under the statistical significance of HRR values at P < 0.01. Furthermore, HRR U0126 values at the lower bound percentile (i.e. 5th and 1st) were statistically significant or very close to the P < 0.05 limit (Additional file 1: Table S4). This indicates that this threshold would be robust enough for infer meaningful co-expression relationships. Several gene co-expression studies in plants have discussed in detail the issue for defining an optimal threshold for gene co-expression, be it from raw PCC or from mutual coexpression ranks (HRR and MR), and its possible solutions [27,28]. These include defining the statistical significance of mutual co-expression ranks [27] or defining a top kth threshold [28] as described earlier. In this study (and for the first time in gene co-expression studies in plants), we have leveraged these two separate approaches by showing that the top 100 HRR for each gene would be a reasonable compromise between manageability of the co-expressed genes list combined with having the statistical power in itsunderlying gene co-expression relationships, and therefore would be suited for downstream guide GCA inferred from the citrus dataset. Furthermore, GO enrichment analysis was then applied to functionally annotate all guide-gene co-expression clusters (30,217 clusters) and assess the predictive performance of these networks (using co-expressed genes within the top 100 HRR for each guide-gene) to recover enriched GO annotations (Table 1). As an alternative approach, when there was no previous knowledge regarding the function of the target gene, identification of densely connected modules based on the graph clustering approach was performed using MCL [13]. The MCL partitions an underlying graph based on the manipulation of transition probabilities or stochastic flows between nodes of the graph. This technique has been shown to effectively identify high-quality functional clusters and is robust to noise [29,30]. Parameter optimisation of the MCL inflation score (I) is often necessary to maximise clustering performance (the quality of derived GO predictions based on specificity, sensitivity and F-measure). Using this approach, we empirically determined that a threshold of HRR30 is a reasonable compromise for MCL given that increasing the threshold to HRR50 (or more) did not improve clustering performance, while reducing to HRR10 improved clustering performance slightly but excluded a greater fraction of probesets (data not shown). Similar observations have been made while determining the optimal HRR value for obtaining biologically relevant clusters [27]. Furthermore, we show that the various HRR scores (10, 20, 30, 40 and 50) used for performance evaluation and parameter optimization of MCL clustering described above were statistically significant (P < 2.95E-04) in all conditions defined (Additional file 1: Table S5). Similar to the guide-gene approach, an evaluation of various inflation parameters on cluster characteristics and clustering performance for each weighted H.