We present the first ab initio calculations of neutrinoless double beta decay matrix elements in A=6-12 nuclei using Variational Monte Carlo wave functions obtained from the Argonne v18 two-nucleon potential and Illinois-7 three-nucleon interaction. We study both light Majorana neutrino exchange and potentials arising from a large class of multi-TeV mechanisms of lepton number violation. Our results provide benchmarks to be used in testing many-body methods that can be extended to the heavy nuclei of experimental interest. In light nuclei we have also studied the impact of two-body short range correlations and the use of different forms for the transition operators, such as those corresponding to different orders in chiral effective theory.