Given a multivariate polynomial $h\left(X_1,\cdots ,X_n\right)$ with integral coefficients verifying an hypothesis of analytic regularity (and satisfying $h\left(\mathbf{0}\right)=1$), we determine the maximal domain of meromorphy of the Euler product ${\prod }_{p\mathrm{prime}}h\left(p^{-s_1},\cdots ,p^{-s_n}\right)$ and the natural boundary is precisely described when it exists. In this way we extend a well known result for one variable polynomials due to Estermann from 1928. As an application, we calculate the natural boundary of the multivariate Euler products associated to a family of toric varieties.