We study the fluctuations of linear statistics with polynomial test functions for Multiple Orthogonal Polynomial Ensembles. Multiple Orthogonal Polynomial Ensembles form an important class of determinantal point processes that include random matrix models such as the GUE with external source, complex Wishart matrices, multi-matrix models and others. Our analysis is based on the recurrence matrix for the multiple orthogonal polynomials, that is constructed out of the nearest neighbor recurrences. If the coefficients for the nearest neighbor recurrences have limits, then we show that the right-limit of this recurrence matrix is a matrix that can be viewed as representation of a Toeplitz operator with respect to a non-standard basis. This will allow us to prove Central Limit Theorems for linear statistics of Multiple Orthogonal Polynomial Ensembles. A particular novelty is the use of the Baker--Campbell--Hausdorff formula to prove that the higher cumulants of the linear statistics converge to zero. We illustrate the main results by discussing Central Limit Theorems for the Gaussian Unitary Ensembles with external source, complex Wishart matrices and specializations of the Schur measure related to multiple Charlier, multiple Krawtchouk and multiple Meixner polynomials.

56 pages; Minor corrections and extended example section