A proof in the parameter range $\,k\geqslant 0\,$ is proposed which relies on the known case $k=1$.
It encompasses the following steps, hereby submitted to the community's critical eye:
- Variable transformation, to simplify the subsequent
- clearing of denominators.
- Scaling the inequality and
- arguing that its $k=1$ instance yields a lower bound for other values of $k$.
- Limit case $k\to 0$
Let's send $\:a,b,c\:$ to the new variables $\:u,v,w\:$ via
$$\begin{pmatrix}u\\ v\\ w\end{pmatrix}\;=\;
\begin{pmatrix}k&1&0\\0&k&1\\1&0&k\end{pmatrix}\,
\begin{pmatrix}a\\b\\c\end{pmatrix}$$
Note that $\,u+v+w=(k+1)(a+b+c)$, and especially
$$\begin{matrix}
\sum_\text{cyc}uv & = & \left(k^2+k+1\right)\sum_\text{cyc}ab & +
& k\:\sum_\text{cyc}a^2 \\[1.5ex]
\sum_\text{cyc}u^2 & = & 2k\:\sum_\text{cyc}ab &
+ & \left(k^2+1\right)\sum_\text{cyc}a^2\,.
\end{matrix}$$
Assume $\,k>0\,$ and combine the preceding two identities to obtain
$$
\left(k+\frac 1 k +1\right)(k+1)^2\sum_\text{cyc}ab\;=\;
\left(k+\frac 1 k\right)\sum_\text{cyc}uv -\sum_\text{cyc}u^2\,.
$$
This is used when multiplying ("Case $k$"), the given inequality, with $\,(k+\frac 1 k -1)(k+1)^2\,u^2v^2w^2$, and after clearing the RHS let
$$
g(k,u,v,w)\;:=\;\left[\left(k+\frac 1 k\right)\sum_\text{cyc}uv -\sum_\text{cyc}u^2 \right]
\cdot\sum_\text{cyc}u^2v^2\; -\;9\left(k+\frac 1 k -1\right)u^2v^2w^2\,.
$$
Note that the freedom to scale the function $g$ has been exploited such that the second "$\sum$" in the brackets does not depend on $k\;\implies$ will drop out upon differentiation w.r.to $k$.
It has to be shown that $\,g(k,u,v,w)\geq 0\;\:\forall\, k,u,v,w>0\,$.
This is accomplished by one-variable analysis in $k$ with $\,u,v,w\,$ fixed.
$g$ enjoys the invariance property $\,g(1/k,\ldots)=g(k,\ldots)\,$ implying
$\,\frac 1kg_k\big(\frac 1k,\ldots\big)=-k\,g_k(k,\ldots)\,$ where $g_k\equiv\frac{\partial}{\partial k}g$. In particular $\,g_k(1,\ldots)=0\,$. Thus the known case $\,g(1,\ldots)\geqslant0\,$ is a bound for $\,g(k,\ldots)\,$, but we badly need $\,g_k(k,\ldots)\geqslant0\,$ for $\,k\geqslant 1\,$ to conclude it's a lower bound!
Now
$$
g_k(k,u,v,w)\;=\;\left(1-\frac 1{k^2}\right)
\left[\left(\sum_\text{cyc}uv\right)\left(\sum_\text{cyc}u^2v^2\right)\; -\;9u^2v^2w^2\right]
$$
and separation of terms involving all three variables from those containing two and distributing the multiples of $\,u^2v^2w^2\,$ accordingly leads to
$$
=\;\left(1-\frac 1{k^2}\right)
\left[uvw\sum_\text{cyc}u(v-w)^2\; +\;\frac{uv+vw+wu}{2}
\sum_\text{cyc}u^2(v-w)^2\right]\,.
$$
Bingo, coz the bracket is a sum of squares with positive coefficients. This completes the $\,k>0\,$ case.
Since the inequality is built on "$\le$" (and not "$<$") the case "$k=0$" holds by continuity too, just send $k$ to $0\,$.