I think you have two main alternatives:

1) Use a MatrixSymbol object with explicit symmetry in its two levels:

In[1]:= A = MatrixSymbol["A", {n, n}, Symmetric[{1, 2}]]

In[2]:= A + Transpose[A] // TensorReduce
Out[2]= 2 MatrixSymbol["A", {n, n}, Complexes, Symmetric[{1, 2}]]

In[3]:= A - Transpose[A] // TensorReduce
Out[3]= SymbolicZerosArray[{n, n}]

In[4]:= Inverse[A] . Transpose[A] // TensorReduce
Out[4]= MatrixPower[MatrixSymbol["A", {n, n}, Complexes, Symmetric[{1, 2}]], 0]

2) or declare assumptions for a symbol B (I use a separate symbol here to avoid Element resolving to True immediately):

In[5]:= $Assumptions = Element[B, Arrays[{n, n}, Symmetric[{1, 2}]]];

In[6]:= B + Transpose[B] // TensorReduce
Out[6]= 2 B

In[7]:= B - Transpose[B] // TensorReduce
Out[7]= SymbolicZerosArray[{n, n}]

In[8]:= Inverse[B] . Transpose[B] // TensorReduce
Out[8]= MatrixPower[B, 0]

Though correct, MatrixPower[squarematrix, 0] should be converted into IdentityMatrix[n].

Thanks Jose -- this is very helpful. A follow-up question, extending the above... Would it be possible using Mathematica 14.1 capabilities to automatically carry out the calculation on the 4th line below, eliminating b, y, A, and Inv[A] in favor of x and c? Here A is an n x n matrix, b, c, x and y are n-vectors. It seems that to do this one needs (Full)Simplify/Eliminate/Reduce type commands that work directly on matrices. Thanks!