Gene co-expression networks have been used successfully for discovering biological relationships among genes on a wholegenome scale, such as predicting gene functional modules and cis-regulatory elements. However, those networks are often constructed in an ad hoc manner, and various methods for network construction and analysis have not been fully evaluated and compared. In this study, we propose a method for constructing gene co-expression networks based on mutual k-nearest-neighbor graphs (mKNN), and compare it with two widely used approaches: threshold-based approach and asymmetric k-nearest-neighbor graph approach (aKNN). We show that mKNN is more robust with respect to the presence of experimental noise and scatter genes, and is less sensitive to parameter variations. Furthermore, we propose a topology-based criterion to guide the selection of the optimal parameter for mKNN, and combine the method with a modularity-based community discovery algorithm to predict functional modules. We evaluate the method on both synthetic and real microarray data. On synthetic data, our method, which does not require any user-tuned parameters, is superior to several popular methods in recovering the embedded modules. Using the yeast stress-response microarray data, we show that the overall functional coherence of the modules predicted by our method using the automatically determined parameters is close to optimal. Finally, we apply the method to study a large collection of gene expression microarray data in Arabidopsis thaliana. Remarkably, with our simple method, we have found many functional modules that are much more significant than those reported by previous studies on the same data set. In addition, we are able to predict cis-regulatory elements for the majority of the functional modules, and the association between the cis-regulatory elements and the functional modules can often be confirmed by existing knowledge.