Let X be a metric space with doubling measure and L a non-negative self-adjoint operator on L 2 (X) whose heat kernels satisfy the Gaussian upper bound estimates. Assume that the growth function ϕ : X×[0, ∞) → [0, ∞) satisfies that ϕ(x, ·) is an Orlicz function and ϕ(·, t) ∈ A ∞ (X) (the class of uniformly Muckenhoupt weights). Let H ϕ, L (X) be the Musielak-Orlicz-Hardy space defined via the Lusin area function associated with the heat semigroup of L. In this article, the authors characterize the space H ϕ, L (X) by means of atoms, non-tangential and radial maximal functions associated with L. In particular, when µ(X) < ∞, the local non-tangential and radial maximal function characterizations of H ϕ, L (X) are obtained. As applications, the authors obtain various maximal function and the atomic characterizations of the "geometric" Musielak-Orlicz-Hardy spaces H ϕ, r (Ω) and H ϕ, z (Ω) on the strongly Lipschitz domain Ω in R n associated with second-order self-adjoint elliptic operators with the Dirichlet and the Neumann boundary conditions; even when ϕ(x, t) := t for any x ∈ R n and t ∈ [0, ∞), the equivalent characterizations of H ϕ, z (Ω) given in this article improve the known results via removing the assumption that Ω is unbounded.(1.2)V(x, λr) ≤ Cλ n V (x, r) 2010 Mathematics Subject Classification. Primary 42B25; Secondary 42B35, 46E30, 30L99. Key words and phrases. Musielak-Orlicz-Hardy space, atom, maximal function, non-negative self-adjoint operator, Gaussian upper bound estimate, space of homogeneous type, strongly Lipschitz domain.