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}_{n\times \infty}$, $n\geq 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(\det\left(I-s^{n}\tau_{n}\right)\right)^{-1}$. The zeta function $\zeta=\prod_{n=1}^{\infty} \left(\det\left(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.