Summary: We develop the multilevel Markov chain Monte Carlo finite element method (MLMCMC-FEM) to sample from the posterior density of the Bayesian inverse problems. The unknown is the diffusion coefficient of a linear, second-order divergence form, elliptic equation in a bounded, polytopal subdomain of \(\mathbb{R}^d\). We provide a convergence analysis with absolute mean convergence rate estimates for the proposed modified MLMCMC-FEM showing in particular error versus work bounds, which are explicit in the discretization parameters. This work generalizes the MLMCMC-FEM algorithm and the error versus work analysis for the uniform prior measure from V. H. Hoang et al. [Inverse Probl. 29, No. 8, Article ID 085010, 37 p. (2013; Zbl 1288.65004)], which we also review here, to linear, elliptic, divergence-form PDEs with a log-Gaussian uncertain coefficient and Gaussian prior measure. In comparison to [loc. cit.], we show by mathematical proofs and numerical examples that the unboundedness of the parameter range under Gaussian prior and the non-uniform ellipticity of the forward model require essential modifications in the MCMC sampling algorithm and in the error analysis. The proposed novel multilevel MCMC sampler applies to general Bayesian inverse problems for linear, second order elliptic divergence-form PDEs with log-Gaussian coefficients. It only requires a numerical forward solver with essentially optimal complexity for producing an approximation of the posterior expectation of a quantity of interest within a prescribed accuracy. Numerical examples using independence and pCN samplers are in agreement with our error versus work analysis.