I've noticed that the equality relation in Peano arithmetic in wikipedia and many other sites, is formalized simply as an equivalence relation closed on the naturals.
Is there a proof that such formulation of PA would prove both substitution schemata of first order logic with equality?
The short answer is that it doesn't really matter.
At a syntactic, proof-theoretic level if you work in FOL without equality, add a binary predicate symbol to your signature, state the PA axioms in terms of this binary predicate symbol, and then prove that this binary predicate symbol is reflexive and satisfies the substitution property, then you know it's consistent to instead take reflexivity and the substitution property as axiom( scheme)s. This is all FOL with equality is doing at the syntactic level.
The main differences then are 1) in one case a failure of reflexivity/substitutivity is an inconsistency whereas in the other case it just means our binary predicate symbol doesn't behave like equality, and 2) in one case a use of the substitution property is a single proof step and in the other case it is a reference to a meta-theorem and corresponds to a potentially quite large derivation. For the way most mathematicians work, these differences don't make any difference.
At a semantic level, the story is only slightly more interesting. The semantics of FOL with equality interprets the equality symbol as the equality of the meta-theory in which we're doing the semantics, ZFC say. A binary predicate symbol will just be interpreted as a binary relation. If that binary predicate symbol satisfies reflexivity and the substitution property, then we can show it is an equivalence relation and every definable subset of the domain is closed under this relation. That is, if $S$ is some definable subset of the domain and $\sim$ is the induced semantic equivalence relation, then $s\in S$ and $s\sim s'$ implies $s'\in S$. The upshot is if $D$ is the semantic domain of a model of the theory where equality is treated as an arbitrary relation with $\sim$ its induced semantic relation, then the quotient set $D/{\sim}$ is a model of the theory where equality is primitive. Therefore the extra models of the former don't change the validity relation.
However, the class of models are distinct in the semantic meta-theory. How these different classes of models relate to other objects of the semantic meta-theory or non-definable subsets is quite different. These are often unnatural questions to ask though. Most of the time it's simpler to consider the smaller class where the equivalence relation is required to be the meta-theory's equality.
You should prove the claims I'm making here yourself.