Joint probabilistic inversions of magnetotelluric (MT) and seismic data has great potential for imaging the thermochemical structure of the lithosphere as well as mapping fluid/melt pathways and regions of mantle metasomatism. In this contribution we present a novel probabilistic (Bayesian) joint inversion scheme for 3D MT and surface-wave dispersion data particularly designed for large-scale lithospheric studies. The approach makes use of a recently developed strategy for fast solutions of the 3D MT forward problem (Manassero et al., 2020) and combines it with adaptive Markov chain Monte Carlo (MCMC) algorithms and parallel-in-parallel strategies to achieve extremely efficient simulations. To demonstrate the feasibility, benefits and performance of our joint inversion method to image the conductivity, temperature and velocity structures of the lithosphere, we apply it to two numerical examples of increasing complexity. The inversion approach presented here is timely and will be useful in the joint analysis of MT and surface wave data that are being collected in many parts of the world. This approach also opens up new avenues for the study of translithospheric and transcrustal magmatic systems, the detection of metasomatised mantle and the incorporation of MT into multi-observable inversions for the physical state of the Earth’s interior.