Modeling genotype-phenotype relationships is a central objective in plant genetics and breeding. Commonly, variations in phenotypic traits are modeled directly in relation to variations at the DNA level, regardless of intermediate levels of biological variation. Here we present an integrative method for the simultaneous modeling of a set of multilevel phenotypic responses to variations at the DNA level. More specifically, for ripe tomato fruits, we use Gaussian graphical models and causal inference techniques to learn the dependencies of 24 sensory traits on 29 metabolites and the dependencies of those sensory and metabolic traits on 21 QTLs. The inferred dependency network which, though not essentially representing biological pathways, suggests how the effects of allele substitutions propagate through multilevel phenotypes. Such simultaneous study of the underlying genetic architecture and multifactorial interactions is expected to enhance the prediction and manipulation of complex traits.