The second order symmetry operators that commute with the Dirac operator with external vector, scalar and pseudo-scalar potentials are computed on a general two-dimensional spin-manifold. It is shown that the operator is defined in terms of Killing vectors, valence two Killing tensors and scalar fields defined on the background manifold. The commuting operator that arises from a nontrivial Killing tensor is determined with respect to the associated system of Liouville coordinates and compared to the the second order operator that arises from that obtained from the unique separation scheme associated with such operators. It shown by the study of several examples that the operators arising from these two approaches coincide.The two-dimensional Dirac equation is currently of great interest due to the connections with graphene's physics [24] and other experimental studies [12]. It is known that the existence of multiplicatively separated solutions of the Dirac equation Dψ = µψ in a given coordinate system and frame implies the additive separability in the same coordinates of the geodesic Hamilton-Jacobi equation. The logical chain of implications works as follows: if the Dirac equation admits multiplicatively separated solutions, then so does the squared Dirac equation; the highest-order terms of this equation coincide with those of the Laplace-Beltrami operator ∆ acting on each component ψ i of the solution. Indeed, it can be shown [27,34] that in this case the Helmholtz equation ∆ψ = µ 2 ψ must admit separable solutions. The multiplicative separation of the last equation is possible only if the same coordinates allow the additive separation of the geodesic Hamilton-Jacobi equation and the Ricci tensor is diagonalized in the same coordinates (Robertson condition)([4] and references therein). There exists examples of separable Schrödinger equations whose corresponding Dirac equations are not separable [16,36]. Among the relevant steps towards a general theory we mention the works by Miller [27] and Shapovalov and Ekle [29] where a theory of complete separation associated with first-order symmetry operators of the Dirac operator is developed. Several contributions to the search for exact solutions of the Dirac equation in curved spaces may be found in Bagrov and Gitman [2] and Cook [18]. In Shishkin [30,31] and Shiskin et al. [1,32,33] an algebraic procedure is developed to obtain separation relations for the solution of the Dirac equation. The definition of separation used therein coincides essentially with our definition of "naive separation" [26] in dimension two. It seems that any definition of separation of variables in two dimensions would yield the same results. In [25] it is proved that in Minkowski spaces all second-order symmetry operators of the Dirac operator (without external fields) can be factorized into products of the first-order operators. It thus seemed that the problem of separation could be reduced to separation associated only to first-order operators. Howewer, Fels and Kamran showed in [23] that...