This work is concerned with zeta functions of two-dimensional shifts of finite type. A two-dimensional zeta function $zeta^{0}(s)$, which generalizes the Artin-Mazur zeta function, was given by Lind for $mathbb{Z}^{2}$-action $phi$. In this paper, the $n$th-order zeta function $zeta_{n}$ of $phi$ on $mathbb{Z}_{ntimes infty}$, $ngeq 1$, is studied first. The trace operator $mathbf{T}_{n}$, which is the transition matrix for $x$-periodic patterns with period $n$ and height $2$, is rotationally symmetric. The rotational symmetry of $mathbf{T}_{n}$ induces the reduced trace operator $tau_{n}$ and $zeta_{n}=left(detleft(I-s^{n}tau_{n}right)right)^{-1}$. The zeta function $zeta=prod_{n=1}^{infty} left(detleft(I-s^{n}tau_{n}right)right)^{-1}$ in the $x$-direction is now a reciprocal of an infinite product of polynomials. The zeta function can be presented in the $y$-direction and in the coordinates of any unimodular transformation in $GL_{2}(mathbb{Z})$. Therefore, there exists a family of zeta functions that are meromorphic extensions of the same analytic function $zeta^{0}(s)$. The natural boundary of zeta functions is studied. The Taylor series for these zeta functions at the origin are equal with integer coefficients, yielding a family of identities, which are of interest in number theory. The method applies to thermodynamic zeta functions for the Ising model with finite range interactions.