A novel approach to zipper fractal interpolation theory for functions of several variables is proposed. We develop multivariate zipper fractal functions in a constructive manner. We then perturb a multivariate function to construct its zipper α-fractal varieties through free choices of base functions, scaling functions, and a binary matrix called signature.In particular, we propose a multivariate Bernstein zipper fractal function and study its approximation properties such as shape preserving aspects, non-negativity, and coordinate-wise monotonicity. In addition, we derive bounds for the graph of multivariate zipper fractal functions by imposing conditions on the scaling factors and the Hölder exponent of the associated germ function and base function. The Lipschitz continuity of multivariate Bernstein functions is also studied in order to obtain estimates for the box dimension of multivariate Bernstein zipper fractal functions.