Abstract. We construct a family of integrable volume-preserving maps in R 3 with a bidimensional heteroclinic connection of spherical shape between two fixed points of saddle-focus type. In other contexts, such structures are called Hill's spherical vortices or spheromaks. We study the splitting of the separatrix under volume-preserving perturbations using a discrete version of the Melnikov method.Firstly, we establish several properties under general perturbations. For instance, we bound the topological complexity of the primary heteroclinic set in terms of the degree of some polynomial perturbations. We also give a sufficient condition for the splitting of the separatrix under some entire perturbations. A broad range of polynomial perturbations verify this sufficient condition. Finally, we describe the shape and bifurcations of the primary heteroclinic set for a specific perturbation.Key words. Separatrix splitting, volume-preserving maps, primary heteroclinic set, Melnikov method, bifurcations AMS subject classifications. 34C37, 34C23, 37C29, 33E201. Introduction. A fundamental question in dynamical systems is the effect that small perturbations of a dynamical system cause on its unperturbed invariant sets. The most studied unperturbed invariant sets are tori and stable/unstable invariant manifolds of hyperbolic sets. Usually, the unperturbed dynamical system is integrable and has separatrices; that is, its stable and unstable invariant manifolds overlap. After a generic perturbation, the perturbed stable and unstable invariant manifolds intersect transversely, which give rise to the onset of chaos, through the creation of Smale horseshoes. This phenomenon is known as the problem of splitting of separatrices. A widely used technique for detecting such intersections is the Melnikov method.Our goal is to apply the Melnikov method to the splitting of separatrices in the discrete volume-preserving framework. Similar questions have been considered before. However, we believe this is the first time that detailed analytical results about the structure of the primary heteroclinic set and its bifurcations are established for specific maps. This represents a step forward with respect to previous works [22,23], in which once written down a formula for the Melnikov function in terms of an infinite series, the approach becomes mainly numerical, because of the technical difficulties that obstruct the analytical one. Here, we have overcome some of these difficulties using basic tools: complex variable theory, quasielliptic functions, homology, and several algebraic tricks. Nevertheless, we have not been able to find an explicit expression of the Melnikov function in terms of elementary functions for any specific perturbation. In opposition, such explicit expressions (in terms of elliptic functions) are known for almost twenty years in the discrete area-preserving setting [17,13].This study is interesting because volume-preserving maps are the simplest, and most natural higher-dimensional versions of the much-studied class o...