Let A :=(A1, A2) be a pair of expansive dilations and φ : Rn×Rm×[0, ∞) → [0, ∞) an anisotropic product Musielak-Orlicz function. In this article, we introduce the anisotropic product Musielak-Orlicz Hardy space HφA(Rn× Rm) via the anisotropic Lusin-area function and establish its atomic characterization, the g-function characterization, the gλ*-function characterization and the discrete wavelet characterization via first giving out an anisotropic product Peetre inequality of Musielak-Orlicz type. Moreover, we prove that finite atomic decomposition norm on a dense subspace of HφA(Rn× Rm) is equivalent to the standard infinite atomic decomposition norm. As an application, we show that, for a given admissible triplet(φ, q, s), if T is a sublinear operator and maps all(φ, q, s)-atoms into uniformly bounded elements of some quasi-Banach spaces B, then T uniquely extends to a bounded sublinear operator from HφA(Rn× Rm) to B. Another application is that we obtain the boundedness of anisotropic product singular integral operators from HφA(Rn× Rm) to Lφ(Rn× Rm)and from HφA(Rn×Rm) to itself, whose kernels are adapted to the action of A. The results of this article essentially extend the existing results for weighted product Hardy spaces on Rn× Rmand are new even for classical product Orlicz-Hardy spaces.