Well, using the familiar identity you pointed out,
$$ {\rm var}(XY) = E(X^{2}Y^{2}) - E(XY)^{2} $$
Using the analogous formula for covariance,
$$ E(X^{2}Y^{2}) = {\rm cov}(X^{2}, Y^{2}) + E(X^2)E(Y^2) $$
and
$$ E(XY)^{2} = [ {\rm cov}(X,Y) + E(X)E(Y) ]^{2} $$
which implies that, in general, ${\rm var}(XY)$ can be written as
$$ {\rm cov}(X^{2}, Y^{2}) + [{\rm var}(X) + E(X)^2] \cdot[{\rm var}(Y) + E(Y)^2] - [ {\rm cov}(X,Y) + E(X)E(Y) ]^{2} $$
Note that in the independence case, ${\rm cov}(X^2,Y^2) = {\rm cov}(X,Y) = 0$ and this reduces to
$$ [{\rm var}(X) + E(X)^2] \cdot[{\rm var}(Y) + E(Y)^2] - [ E(X)E(Y) ]^{2} $$
and the two $[ E(X)E(Y) ]^{2}$ terms cancel out and you get
$$ {\rm var}(X){\rm var}(Y) + {\rm var}(X)E(Y)^{2} + {\rm var}(Y)E(X)^{2} $$
as you pointed out above.
Edit: If all you observe is $XY$ and not $X$ and $Y$ separately, then I don't think there is a way for you to estimate ${\rm cov}(X,Y)$ or ${\rm cov}(X^2,Y^2)$ except in special cases (for example, if $X,Y$ have means that are known a priori)