Since we have a homogeneous inequality, we may assume $xyz=1$ without loss of generality, then find the infimum of $AM(x,y,z)+HM(x,y,z)$ over the set
$$D=\left\{(x,y,z)\in\mathbb{R}_{+}^3: xyz=1\right\}.$$
So we want:
$$ \inf_{\mathbb{R}^+\times\mathbb{R}^+}\left(\frac{x+y+\frac{1}{xy}}{3}+\frac{3}{xy+\frac{1}{x}+\frac{1}{y}}\right)=\inf_{\mathbb{R}^+\times\mathbb{R}^+}\left(\frac{x^2 y+xy^2+1}{3xy}+\frac{3xy}{x^2y^2+x+y}\right).$$
Let $f(x,y)$ be the last function. It is a symmetric function, hence it is enough to study its behaviour over the set $0<y\leq x$. So, let $y=kx$ with $k\in(0,1]$. We have:
$$ g_k(x)=f(x,kx) = \frac{1}{3kx^2}+\frac{(k+1)x}{3}+\frac{3kx}{1+k+k^2 x^3} $$
and disregarding the last (positive) term, by the AM-GM inequality:
$$ g_k(x)\geq \frac{1}{3kx^2}+\frac{(k+1)x}{6}+\frac{(k+1)x}{6} \geq 3\sqrt[3]{\frac{(k+1)^2}{108k}}$$
where the last term is a decreasing function over $(0,1]$. We just need a little improvement of the last inequality to prove the original claim. So, let us study $f(x,y)$ under another perspective. Assuming $x+y=k$ we have:
$$ f(x,y) = \frac{k+\frac{1}{xy}}{3}+\frac{3}{\frac{k}{xy}+xy}=\frac{k+\frac{1}{x(k-x)}}{3}+\frac{3}{\frac{k}{x(k-x)}+x(k-x)} $$
and by putting $x=\lambda k, y=(1-\lambda)k$ we are left with:
$$ \frac{k+\frac{1}{k^2\lambda(1-\lambda)}}{3}+\frac{3}{\frac{1}{k\lambda(1-\lambda)}+k^2\lambda(1-\lambda)}=\color{red}{\frac{k+\frac{1}{k^2 L}}{3}+\frac{3}{\frac{1}{kL}+k^2 L}=h(L,k)}$$
where $L$ ranges over $\left(0,\frac{1}{4}\right)$ and $k$ ranges over $\mathbb{R}_+$. The stationary points of $h(L,k)$ for a fixed $k$ occur at
$$ L \in\left\{\pm\frac{1}{\sqrt{2} k^{3/2}},\pm\frac{1}{\sqrt{5} k^{3/2}}\right\} $$
and
$$ h\left(\frac{1}{\sqrt{2} k^{3/2}},k\right) = \frac{4}{3}\sqrt{\frac{2}{k}}+\frac{k}{3},\qquad h\left(\frac{1}{\sqrt{5} k^{3/2}},k\right) = \frac{5}{6}\sqrt{\frac{5}{k}}+\frac{k}{3},$$
so the AM-GM inequality proves the claim.
In facts, the original inequality is optimal and our approach shows that equality is attained by $$ \color{red}{(x,y,z) = \lambda\cdot (1,4,4)} $$
and cyclic shifts.