In this work, we extend the hybrid Cherno tau-leap method to the multilevel Monte Carlo (MLMC) setting. Inspired by the work of Anderson and Higham on the tau-leap MLMC method with uniform time steps, we develop a novel algorithm that is able to couple two hybrid Cherno tau-leap paths at di erent levels. Using dual-weighted residual expansion techniques, we also develop a new way to estimate the variance of the di erence of two consecutive levels. This is crucial because the computational work required to stabilize the coecient of variation of the sample variance estimator of the di erence between two consecutive levels is often una ordable for the deepest levels of the MLMC hierarchy. Our algorithm enforces the total error to be below a prescribed tolerance,TOL, with high probability. This is achieved with nearly optimal computational work. Indeed, the computational complexity of our method is of order O (TOL^-2), the same as with an exact method, but with a smaller constant. Our numerical examples show substantial gains with respect to the previous single-level approach and the Stochastic Simulation Algorithm.