$\def\ms{\boldsymbol{\sigma}}
\def\s{\sigma}
\def\r{\rho}
\def\f{\varphi}
\def\o{\cdot}
\def\ra{\rightarrow}
\def\ve{{\bf e}}
\def\vv{{\bf v}}
\def\MM{{\bf M}}
\def\id{\mathbb{I}}
\newcommand\dv[1]{\nabla\cdot #1}
\newcommand\mt[4]{\left[\begin{array}{cc}
#1 & #2 \\
#3 & #4
\end{array}\right]}
\newcommand\cvc[2]{\left[\begin{array}{cc} #1 & #2 \end{array}\right]}
\newcommand\vc[2]{\left[\begin{array}{c} #1 \\ #2 \end{array}\right]}$Consider $\ms$ in the natural basis,
$$\ms = \mt{\s_{xx}(x,y)}{\s_{xy}(x,y)}{\s_{yx}(x,y)}{\s_{yy}(x,y)}.$$
We can write this as
$$\ms = \s_{ij}\ve_i \ve_j^T,$$
where $\ve_x=\cvc{1}{0}^T$ and
$\ve_y=\cvc{0}{1}^T$.
We can relate the components of $\ms$ in different orthonormal bases in the following way.
We have
\begin{align*}
\ms=\s_{ij}\ve_i \ve_j^T = \s_{i'j'}\ve_{i'}\ve_{j'}^T, \tag{1}
\end{align*}
and so
$$\s_{ij} = \ve_i^T \ve_{i'} \s_{i'j'} \ve_{j'}^T \ve_j$$
or
$$\ms = \MM^T\ms' \MM,$$
where $\ms'=[\s_{i'j'}]$ is the matrix of components in the primed basis and where
$$M_{i'i} = \ve_{i'}^T \ve_i.$$
By assumption the bases are orthonormal, so $\MM^T = \MM^{-1}$.
(Note that, by (1) and the fact that the Kronecker delta is unchanged by coordinate transformation,
$\id = \ve_i \ve_i^T = \ve_{i'}\ve_{i'}^T$.
Thus,
$[\MM^T\MM]_{ij}
= \ve_{i}^T \ve_{i'} \ve_{i'}^T\ve_j
= \ve_i^T \ve_j
= [\id]_{ij}.$)
A similar argument can be made to show that the components of a vector $\vv$ in the different bases are related by
\begin{align*}
\vv' &= \MM \vv, \tag{2}
\end{align*}
where $\vv$ and $\vv'$ are the vectors whose components are in the natural and primed bases, respectively.
(We assume that $\MM=\MM(x',y')$.)
We write the divergence in the primed basis as $[\dv\ms]'$.
Note that this quantity is a vector, and so transforms as indicated in (2).
Thus,
\begin{align*}
[\dv\ms]' &= \MM[\dv\ms]_{(x,y)\ra(x',y')} \\
&= \MM[\dv(\MM^T\ms' \MM)_{(x',y')\ra(x,y)}]_{(x,y)\ra(x',y')}.
\end{align*}
For this problem we have
\begin{align*}
[\dv\ms]' &= \MM[\dv(\MM^T\ms' \MM)_{(\r,\f)\ra(x,y)}]_{(x,y)\ra(\r,\f)}.\tag{3}
\end{align*}
Note that
\begin{align*}
\ve_\r &= \cos\f \,\ve_x + \sin\f \,\ve_y \\
\ve_\f &= -\sin\f \,\ve_x + \cos\f \,\ve_y.
\end{align*}
(This basis is orthonormal by inspection.)
This implies, for example, that
$$M_{\r x} = \ve_\r^T \ve_x
= \cvc{\cos\f}{\sin\f}
\vc{1}{0} = \cos\f.$$
Calculating the other components, one finds
$$\MM = \mt{\cos\f}{\sin\f}{-\sin\f}{\cos\f}$$
or
$$\MM_{(\r,\f)\ra(x,y)} = \frac{1}{\sqrt{x^2+y^2}}
\mt{x}{y}{-y}{x}.$$
Note also that
$$[\ms'(\r,\f)]_{(\r,\f)\ra(x,y)}
= \ms'(\sqrt{x^2+y^2},\arctan y/x).$$
It is now a straightforward, if tedious, task to work out the correct form for $[\dv\ms]'$.
Addendum
(* mma code to check (3) *)
fs[foo_] := FullSimplify[foo, {r > 0, -Pi < f < Pi}];
rpc = {r -> Sqrt[x^2 + y^2], f -> ArcTan[x, y]};
rcp = {x -> r Cos[f], y -> r Sin[f]};
Apolar = {{srr[r, f], srf[r, f]}, {sfr[r, f], sff[r, f]}};
M = {{Cos[f], Sin[f]}, {-Sin[f], Cos[f]}};
Mt = Transpose[M];
M //. rpc // fs
Apolar[[1]][[1]] //. rpc
(* out: [a check on some relations above] *)
ansmma = Div[Apolar, {r, f}, "Polar"] // fs
(* out: [divergence according to mma] *)
ans = M.Div[Mt.Apolar.M //. rpc, {x, y}, "Cartesian"] //. rcp // fs
(* out: [divergence using (3)] *)
ans == ansmma // fs
(* out: True *)
