Parallel transport as algebra isomorphism
It is a classic theorem in linear algebra that any ($\mathbb{R}$-linear) automorphism $\phi$ of the the ring $M_n(\mathbb{R})$ is inner, i.e., of the form $\phi(x) = axa^{-1}$ for some invertible $a\in M_n(\mathbb{R})$. In particular, the group of automorphisms of the algebra is $\mathrm{PGL}(n,\mathbb{R})$, a simple group of dimension $n^2{-}1$, whose maximal compact subgroups are all conjugate to $\mathrm{PO}(n) = \mathrm{O}(n)/\{\pm I_n\}$.
To avoid trivial cases, assume that $n>1$ and, at least to start with, assume also that $M$ be simply connected. If $M^{n^2}$ carried a Riemannian metric $g$ whose Levi-Civita connection $\nabla$ preserved an $M_n(\mathbb{R})$ algebra structure on $TM$, then, because its holonomy would be compact, there would be an underlying $\nabla$-parallel $\mathrm{PO}(n)$-structure on $M$.
Now, the representation of $\mathrm{PO}(n)$ on $M_n(\mathbb{R})$ is reducible, being the sum of three terms:
$$
M_n(\mathbb{R}) = \mathbb{R}{\cdot}I_n\ \oplus\ A_n(\mathbb{R})\ \oplus\ S_n(\mathbb{R}),
$$
where $A_n(\mathbb{R})\subset M_n(\mathbb{R})$ consists of the skew-symmetric $n$-by-$n$ matrices, while $S_n(\mathbb{R})\subset M_n(\mathbb{R})$ consists of the traceless symmetric $n$-by-$n$ matrices.
When $n>2$, the representations of $\mathrm{PO}(n)$ on both $A_n(\mathbb{R})$ and $S_n(\mathbb{R})$ are irreducible and (almost) faithful. By the Bianchi identities, the two corresponding subbundles of $TM$ will be $\nabla$-parallel and, hence by the reducibility of the holonomy, the metric on $M$ will locally split as a product into three factors The (almost) faithfulness of two of the representations guarantees (by the second Bianchi identity), that the curvature of $\nabla$ must vanish identically. Hence, the only case that occurs is the locally flat structure. One can, of course, do things with open sets and discrete quotients, but, locally, the flat case is the only case.
When $n=2$, things are different. The group $\mathrm{PO}(2)$ acts on $A_2(\mathbb{R})\simeq \mathbb{R}$ as $\{\pm1\}\simeq \mathbb{Z}_2$, and it acts on $S_2(\mathbb{R})\simeq \mathbb{R}^2$ as $\mathrm{O}(2)$. Since $M$ is simply connected, it follows that, at least locally, $(M,g)$ can be written as a metric product $$ M = \mathbb{R}\times \mathbb{R}\times \Sigma, $$ where $(\Sigma,h)$ is an oriented Riemannian surface, which, hence, has an associated orthogonal complex structure $J$ and area form $\omega$. Then one can define an $M_2(\mathbb{R})$ structure on $T_pM$ = $T_x\mathbb{R}\times T_y\mathbb{R}\times T_z\Sigma = \mathbb{R}\oplus\mathbb{R}\oplus T_z\Sigma$ by the rule $$ (a,b,v)(a',b',v') = (aa'{-}bb'{+}h(v,v'),\ ab'{+}a'b{+}\omega(v,v'),\ av'{+}a'v{+}b\,Jv' {+} b'\,Jv). $$ Thus, the general solution essentially depends on a Riemannian metric in dimension $2$, which is one function of two variables, up to diffeomorphism.
If one is willing to consider pseudo-Riemannian metrics instead of only Riemannian ones, there are other solutions. For example, when $n=2$, $\mathrm{PGL}(2,\mathbb{R})$ acts preserving the irreducible splitting $$ M_2(\mathbb{R}) = \mathbb{R}{\cdot}I_2\ \oplus Z_2(\mathbb{R}) $$ where $Z_2(\mathbb{R})\simeq \mathbb{R}^3$ is the space of $2$-by-$2$ traceless matrices. The action of $\mathrm{PGL}(2,\mathbb{R})$ on $Z_2(\mathbb{R})$ preserves a Lorentzian inner product $\langle,\rangle$ that satisfies $\langle v,v\rangle = \det(v)$, so that $v^2 = -\langle v,v\rangle I_2$ and the 'outer' or 'skew' product $[v,w] = vw-wv$.
Consequently, if $(M^4,g)$ is a simply-connected Lorentzian $4$-manifold whose holonomy preserves a $M_2(\mathbb{R})$ algebra structure on $TM$, then $M$ can be written locally as a metric product $$ M^4 = \mathbb{R} \times \Sigma^3 $$ where $(\Sigma^3,h)$ is an oriented Lorentzian $3$-manifold and $g$-parallel algebra structure on $TM= T\mathbb{R}\oplus T\Sigma$ can be written in the form $$ (a,v)(a',v') = (aa'-h(v,v'),\ av' + a'v + v{\times}v') $$ where $v{\times}v' = \ast_h(v\wedge v')$ and where $\ast_h:\Lambda^2(T\Sigma)\to T\Sigma$ is the Hodge star associated to the Lorentzian metric $h$ and the orientation. This gives a family of (Lorentzian) solutions depending on $3$ functions of $3$ variables locally.
Finally, when $n>2$, there is the natural bi-invariant $M_n(\mathbb{R})$ structure on $M=\mathrm{GL}(n,\mathbb{R})$ itself, regarded as a pseudo-Riemannian symmetric space, together its dual symmetric space $M^* = \mathrm{GL}(n,\mathbb{C})/\mathrm{GL}(n,\mathbb{R})$. By Berger's classification of non-symmetric holonomies (or by simply computing the space of curvatures of this particular $\mathrm{PGL}(n,\mathbb{R})$-structure and seeing that the only possibility is the locally symmetric one), these are, locally, the only examples when $n>2$.
Every Riemannian metric on the sphere has holonomy $SO(n)$, since the sphere is not a product and not Kaehler unless $n=2$ and has no parallel spinors in any metric. So on the sphere, no Riemannian metric has parallel algebra structure on the tangent bundle.
I think you need to restrict your question a bit. When $M$ is $n^2$ dimensional Riemannian flat Euclidean space with Levi-Civita connection, then the parallel transport induced by this is the identity on $M_n(\mathbb{R})$. And thus is an algebra isomorphism.
Perhaps you want to look into the notion of Frobenius manifold:
https://en.wikipedia.org/wiki/Frobenius_manifold
There are non-trivial Frobenius manifolds. I think this partially answers maybe what you wanted to ask.