Inferring the underlying causal relationships between genes from expression data is a task of critical importance in systems biology. In the particular case where a mixture of observational and intervention experiments (ex: single or multiple knock-out or knock-down experiments) are available, only a few methods are currently available. The first, called Intervention calculus when the Directed graph is Absent (IDA; Maathuis et al., 2009), provides causal bounds for direct and indirect effects once a skeleton graph has been estimated using the PC algorithm (pcalg R package). The second approach instead relies on the notion of a causal ordering of genes, whose posterior distribution is inferred from the data using a probabilistic generative model which allows for single and multiple interventions. This can be conveniently done within an MCMC simulation (Rau et al., 2013), in particular for the case where the posterior distribution is efficiently approximated (Hartmann and Nuel, 2017). In the present work, we extend this second approach by introducing two novelties: 1) we use a Laplace approximation to obtain a fast approximation of the integrated likelihood of the model; and 2) we use parallel tempering combined with the classic MC3 algorithm (Barker et al., 2010) to efficiently explore the directed acyclic graph (DAG) space. This new approach proves to be both faster and more reliable than the previous causal approach based on causal node orderings. It also has the advantage of providing a collection of DAG structures that can be aggregated to provide robust estimates for each dataset. We also introduce a simple mixture model over the DAG space to help represent this posterior distribution. Finally, we illustrate the method on both simulated and real datasets, where it shows promising results.