What does the Pythagorean Theorem really prove?
Let $V$ be a real or complex vector space equipped with an inner product $\langle -, - \rangle$, and let $\| - \|$ be the corresponding norm. One version of the Pythagorean theorem states that if $v_1, \dots v_n$ are mutually orthogonal vectors in $V$ (meaning that $\langle v_i, v_j \rangle = 0$ for $i \neq j$), then
$$\left\| \sum v_i \right\|^2 = \sum_i \| v_i \|^2.$$
This theorem has content even in the case that $V = \mathbb{R}^d$ and $\langle -, - \rangle$ is the dot product, the point being that $v_1, \dots v_n$ are arbitrary and not necessarily aligned with the coordinate axes. It has various uses in other more exotic cases; for example, for random variables of mean zero equipped with the covariance inner product, this statement is linearity of variance. And of course there is the version with infinitely many $v_i$ for use in Hilbert spaces, e.g. to prove Parseval's identity.
From this coordinate-free point of view you can't even define rotations and reflections without a choice of inner product, so it's unclear how to interpret your last question.
The basic reference for the significance of the Pythagorian theorem is the illuminating article by Givental:
Givental, Alexander. "The Pythagorean theorem: What is it about?" Amer. Math. Monthly 113 (2006), no. 3, 261–265.
I rather love the question you linked concerning the behavior of triangles in non-Euclidean geometry. This may not be exactly what you are looking for, but I think it really forms the basis for why the Pythagorean theorem is special. The other answer gives a cool algebraic answer in some sense, so here is a geometric one.
The basic idea is this:
Euclidean space is the "sweet point" where "straight lines" (geodesics) satisfy the Pythagorean theorem. This is because it is the "transition point" between positively and negatively curved space (i.e. where there is no curvature).
Suppose we have a point $x$ and start two geodesics moving at constant unit speed from each other with angle $\theta$. If you look in Villani, Optimal transport, old and new, eq 14.1 (also see this question), for 2-manifolds, we get the distance $d(t)$ between the two moving points at time $t$ to be: $$ d(t) = t\sqrt{2(1-\cos(\theta))} \left[ 1 - \frac{\kappa(x)\,\cos^2(\theta/2)}{6}t^2 + O(t^4) \right] $$ where $\kappa(x)$ is the Gaussian curvature at $x$.
Let's construct a right triangle, with two sides of equal length $t$ and measure the hypotenuse as the distance between them. Assume $t<<1$. Take $\theta=\pi/2$. Then: $$ d(t) = t\sqrt{2}\left[ 1 - \frac{\kappa(x)\,t^2}{8} \right] $$ is the hypotenuse length (i.e. $c=d(t)$).
Recall that the other two sides of our triangle are of length $a=b=t$. Notice that if $\kappa=0$, we get: $$ a^2+ b^2 = t^2+t^2 = 2t^2=d(t)^2 = c^2 $$ so on an uncurved manifold, Pythagorean theorem is correct!
Note that everything depends on the curvature of the space. If $\kappa(x)<0$, then $a^2+b^2 < c^2$; if $\kappa(x)>0$, then $a^2+b^2 > c^2$. Indeed, the degree to which Pythagoras is accurate decreases as space becomes more curved.
So, my answer to what is the underlying truth of the Pythagorean Theorem, in modern mathematical terms, is that I would consider the Pythagorean theorem to be the special case result of geodesic distances in uncurved space.
This holds somewhat more generally for Riemannian manifolds: negative Ricci scalar curvature causes geodesics to diverge, while positive does the opposite. Hence, if we construct "triangles" by treating "straight lines" to be geodesics and measure the distance between them as they grow, in the former case, the distance will "overshoot" the Pythagorean theorem; in the latter case, it will undershoot. (There's a cool little graph of that in here, albeit in a probabilistic sense). Only in the flat Euclidean case (at least locally), will they be exactly equal all the time.
Maybe some people will also want to note that, in geodesic normal coordinates, the metric is locally euclidean (i.e. $g_{ij} \approx \delta_{ij}$, so the Pythagorean theorem will hold infinitesimally (since the square geodesic distance will be $s^2=g_{ij}x^ix^j \approx \delta_{ij}x^ix^j=(x^1)^2 + (x^2)^2$ in 2D).