You're putting the cart before the horse here. First you need to decide what your transformation does to both the independent and dependent variables of your PDE, and then you can decide whether this transformation is a symmetry or not.
For equations like the diffusion equation, the definition of the time-reversal transformation is $x\mapsto x$, $t\mapsto -t$, and $\rho \mapsto \rho$. This is fairly easy to motivate - if you're watching a video of some concentration $\rho$ evolve with time and then you pause and watch the movie backward, then nothing happens to the concentration (e.g. it doesn't suddenly acquire a minus sign, like the velocities of the particles do). With this definition, the diffusion equation is time-reversal symmetric if, for all solutions $\rho(x,t)$, we have that $\rho(x,-t)$ is also a solution. Since this is clearly false, the diffusion equation is not time-reversal symmetric, while equations like the elementary wave equation are.
On the other hand, in quantum mechanics we define time-reversal symmetry to be implemented by an antiunitary operator, which implies that the wave function be mapped to its complex conjugate. The reason we do this is ultimately because we want to preserve the canonical commutation relations, and we cannot do this if the time-reversal operator does not include a complex conjugation.
All that being said, you are asking what $f$ would have to be in order for the diffusion equation to have $(x,t,\rho)\mapsto(x,-t,f\circ\rho)$ as a symmetry. This is not an uninteresting question, but it isn't about time-reversal. Incidentally, the answer is only the trivial map $f\circ \rho =0$. You can see this intuitively by noting that e.g. in free space, all solutions to the diffusion equation asymptotically approach the same steady state - namely the uniform distribution - as $t\rightarrow \infty$. This would serve as the $t\rightarrow -\infty$ initial condition for $\sigma(x,t):=(f\circ \rho)(x,-t)$, which would approach a localized distribution as $t\rightarrow 0^-$.
However, if $\sigma$ is also a solution to the diffusion equation, then this doesn't make sense. If $\sigma$ is initially uniform, then it will stay uniform forever. Ultimately, this is a result of the fact that the diffusion equation is first-order in time, which implies that you only need to supply the initial value of the function to specify its evolution. Since you don't need to supply its initial time-derivative, the solution has no "inertia", and the time-evolution of the solutions all follow the same generic pattern (i.e. they spread out, rather than become localized).
I think I don't seem to get what you mean by "This would serve as the $t\rightarrow -\infty$ initial condition...". Can you please explain?
Pick an arbitrary solution $\rho(x,t)$, with specified initial condition $\rho(x,0)$, and evolve it forward in time to find $\rho(x,T)$. Now define $\sigma(x,t):=f\big(\rho(x,-t)\big)$, which must also be a solution to the diffusion equation with initial condition $\sigma(x,-T) = f\big(\rho(x,T)\big)$.
This must hold for arbitrary solutions, so let's pick a specific one - $\rho(x,0) = c \delta(x)$ (or if you don't like delta functions, pick an arbitrarily well-localized Gaussian). At time $t$, we would have
$$\rho(x,t) = \frac{1}{\sqrt{4\pi t}} e^{-x^2/4t}$$
Our transformed function is then
$$\sigma(x,t) = f\big(\rho(x,-t)\big) = f\left( \frac{1}{\sqrt{-4\pi t}} e^{x^2/4t}\right), \qquad t\in [-T,0]$$
It should be clear that $\rho$ evolves from a highly localized distribution to a delocalized distribution. Since the spatial dependence of $\sigma$ comes through $\rho$, it should also be clear that no matter what $f$ is (unless it's trivial), $\sigma$ evolves the other way - from a delocalized distribution to a localized one.
But that's not what solutions to the diffusion equation do. They exhibit a very specific trend - localized to delocalized - as time goes forward. As I mentioned in a comment, this one-way time evolution means that a snapshot of $\rho$ at any time $t$ is sufficient to determine its future evolution, which is mathematically reflected in the fact that the diffusion equation is first-order in time.
Contrast this with solutions to the wave equation, for which you need to know the solution and its time derivative as initial conditions. These solutions have inertia in a certain sense - if you flip the initial time derivative, the wave will just propagate in the opposite direction while still satisfying the wave equation.
A rigorous proof could be drawn up for this, but I find this argument sufficiently convincing.
