  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.
