Trait evolution models allow for the estimation of evolutionary correlations between a set of traits observed in a sample of related organisms. By directly modeling the evolution of the traits on a phylogenetic tree in a Bayesian framework, the model's structure allows us to control for shared evolutionary history between the organisms in the sample and avoid spurious inference. However, relevant correlations are assessed through a heuristic procedure based on correlation high posterior density intervals. In order to employ a model based method to identify the underlying association structure between the variables we explore the use of Gaussian graphical models (GGM) for covariance selection in this context. This association structure can be translated by the conditional independence embedded in the inverse covariance, the precision matrix K, such that zero off-diagonal entry (i,j), implies that the corresponding variables are conditionally independent given all the other variables in the model. In GGM's the zero entries in the off-diagonal of the precision also correspond to the non-existing edges in the undirected graph G, that models the dependence structure directly. We model K|G with a G-Wishart conjugate prior which results in sparse precision matrix estimates. The sparsity pattern on the across-trait precision matrix reflects on the corresponding correlation matrix potentially shrinking some of its off-diagonals to zero as desired. We evaluate our approach through Monte Carlo simulation studies and compare the results to the standard method. We find smaller precision and correlation mean square errors (MSE) in sparse models for entries between variables simulated as conditionally independent, and similar MSE's between sparse and standard models for variables simulated as conditionally dependent. We also test our approach to examine the association structure of prokaryotic genomic and phenotypic traits. Our approach provides a systematic solution for elimination of spurious correlations and better inference for the across-trait precision and correlation matrices, especially on conditionally independent variables which are the target for sparsity. Although keeping the Gibbs sampling efficiency from G-Wishart conjugacy on precision matrix, the computational burden of including GGM’s in the model increases at the cost of G-Wishart normalizing constant evaluations performed during graph estimation.