We consider a three-dimensional (3D) lattice model associated with the intertwiner of the quantized coordinate ring A q (sl 3 ), and introduce a family of layer to layer transfer matrices on m × n square lattice. By using the tetrahedron equation we derive their commutativity and bilinear relations mixing various boundary conditions. At q = 0 and m = n, they lead to a new proof of the steady state probability of the n-species totally asymmetric zero range process obtained recently by the authors, revealing the 3D integrability in the matrix product construction.