High-throughput genomics allows genome-wide quantification of gene expression levels in tissues and cell types and, when combined with sequence variation data, permits the identification of genetic control points of expression (expression QTL or eQTL). Clusters of eQTL influenced by single genetic polymorphisms can inform on hotspots of regulation of pathways and networks, although very few hotspots have been robustly detected, replicated, or experimentally verified. Here we present a novel modeling strategy to estimate the propensity of a genetic marker to influence several expression traits at the same time, based on a hierarchical formulation of related regressions. We implement this hierarchical regression model in a Bayesian framework using a stochastic search algorithm, HESS, that efficiently probes sparse subsets of genetic markers in a high-dimensional data matrix to identify hotspots and to pinpoint the individual genetic effects (eQTL). Simulating complex regulatory scenarios, we demonstrate that our method outperforms current state-of-the-art approaches, in particular when the number of transcripts is large. We also illustrate the applicability of HESS to diverse real-case data sets, in mouse and human genetic settings, and show that it provides new insights into regulatory hotspots that were not detected by conventional methods. The results suggest that the combination of our modeling strategy and algorithmic implementation provides significant advantages for the identification of functional eQTL hotspots, revealing key regulators underlying pathways.T HE current focus of biological research has turned to high-throughput genomics, which encompasses largescale data generation and a variety of integrated approaches that combine two or more -omics of data sets. An important example of integrative genomics analysis is the investigation of the genetic regulation of transcription, also called expression quantitative trait locus (eQTL) or "genetical genomics" studies (Cookson et al. 2009;Majewski and Pastinen 2011). A typical eQTL analysis follows a natural structure of parallel regressions between the large set of q responses (i.e., expression phenotypes), and that of p explanatory variables (i.e., genetic markers, often single nucleotide polymorphism, SNPs), where p is typically much larger than the number of observations n.From a statistical point of view, the size and the complex multidimensional structure of eQTL data sets pose a significant challenge. Not only does one wish to detect the set of important genetic control points for each response (expression phenotype), including cis-and trans-acting control for the same transcript, but, ideally, one would wish to exploit the dependence between multiple expression phenotypes. This will facilitate the discovery of key regulatory markers, so-called hotspots (Breitling et al. 2008), i.e., genetic loci or single polymorphisms that influence a large number of transcripts. Identification of hotspots can inform on network and pathways, which are likel...