Here is a continuous counterexample. It has discrete analogues, some of which are given at the end. We start with a general remark, which may help understand where the counterexamples come from. (For a short version of the answer, go to the picture and ponder it, really it says everything.)
Assume that $(X,Y)$ has PDF $p$ and that there exists some measurable function $q$ such that, for every $(x,y)$, $ p(x,y)+p(y,x)=2q(x)q(y). $ One can assume without loss of generality that $q$ is a PDF. Then, for every function $f$, $ E[f(X)f(Y)]=\iint f(x)f(y)p(x,y)\mathrm dx\mathrm dy=\iint f(x)f(y)p(y,x)\mathrm dx\mathrm dy, $ hence $ E[f(X)f(Y)]=\iint f(x)f(y)q(x)q(y)\mathrm dx\mathrm dy=\left(\int f(x)q(x)\mathrm dx\right)^2. $ Thus, if, furthermore, $ q(x)=\int p(x,y)\mathrm dy=\int p(y,x)\mathrm dy, $ then indeed, $ E[f(X)f(Y)]=E[f(X)]E[f(Y)]. $ Now, a specific counterexample: assume that $(X,Y)$ is uniformly distributed on the set $ D=\{(x,y)\in[0,1]^2\,;\,\{y-x\}\leqslant\tfrac12\}, $ where $\{\ \}$ is the function fractional part. In words, $D$ (the black part in the image of the square $[0,1]^2$ below) is the union of the parts of the square $[0,1]^2$ between the lines $y=x+\frac12$ and $y=x$, and below the line $y=x-\frac12$.
$\hskip2in$
Then $(Y,X)$ is uniformly distributed on the image of $D$ by the symmetry $(x,y)\to(y,x)$ (the white part in the image of the square $[0,1]^2$ above), which happens to be the complement of $D$ in the square $[0,1]^2$ (actually, modulo some lines, which have zero Lebesgue measure). Thus, our first identity holds with $q=\mathbf 1_{[0,1]}$, that is:
If $(X,Y)$ is uniformly distributed on $D$, then $X$ and $Y$ are both uniformly distributed on $[0,1]$, $(X,Y)$ is not independent, and, for every function $f$, $E[f(X)f(Y)]=E[f(X)]E[f(Y)]$.
Note that $(X,Y)$ can be constructed as follows. Let $U$ and $V$ be i.i.d. uniform on $[0,1]$, then $(X,Y)=(U,V)$ if $(U,V)$ is in $D$, else $(X,Y)=(V,U)$.
An analogue of this in the discrete setting is to consider $(X,Y)$ with joint distribution on the set $\{a,b,c\}^2$ described by the matrix $ \frac19\begin{pmatrix}1&2&0\\0&1&2\\2&0&1\end{pmatrix}. $ Then the random set $\{X,Y\}$ is distributed like $\{U,V\}$ where $U$ and $V$ are i.i.d. uniform on $\{a,b,c\}$. For example, considering $S=\{(a,b),(b,a)\}$, one sees that $ P[(X,Y)\in S]=\tfrac29=P[(U,V)\in S], $ since $[(X,Y)\in S]=[(X,Y)=(a,b)]$, while $ P[(U,V)\in S]=P[(U,V)=(a,b)]+P[(U,V)=(b,a)]=\tfrac19+\tfrac19. $ Thus, $ E[f(X)f(Y)]=E[f(U)f(V)]=E[f(U)]E[f(V)]=E[f(X)]E[f(Y)]. $ This example can be extended to any sample space of odd size. A more general distribution on $\{a,b,c\}^3$ is, for every $|t|\leqslant1$, $ \frac19\begin{pmatrix}1&1+t&1-t\\1-t&1&1+t\\1+t&1-t&1\end{pmatrix}. $