We consider a fully-implicit formulation for two-phase flow in a porous medium with capillarity, gravity, and compressibility in three dimensions. The method is implicit in time and uses the multiscale mortar mixed finite element method for spatial discretization in a non-overlapping domain decomposition context. The interface conditions between subdomains are enforced in terms of Lagrange multiplier variables defined on a mortar space. There are two novelties in our approach: first, we linearize the coupled system of subdomain and mortar variables simultaneously to form a global Jacobian and eliminate variables by taking Schur complements; and second, we adapt a two-stage preconditioning strategy to solve the resulting formulation. The new formulation exploits the local invertibility of mortar block matrices to obtain a system with subdomain variables as unknowns. The two-stage preconditioner uses the Householder transformation to decouple the pressure and saturation variables for each grid block. This allows the elliptic equation for the pressure and hyperbolic equation for the saturation to be solved in a decoupled manner. The algorithm is fully parallel and numerical tests show the low computational costs of eliminating the mortar variables to obtain the final Jacobian matrices. We have demonstrated parallel scalability using our two-stage preconditioner on large-scale heterogeneous two-phase reservoir simulation problems with over 10 million elements and over 1000 processors.