Zygote may be running into an issue that can also be understood by considering the analytical gradient
\frac{\partial}{\partial x} x \tanh(e^x) = \tanh(e^x) + \frac{x e^x}{\cosh(e^x)^2}
The NaN will come from the second term because e^x will overflow at an x large enough for \cosh(x) to overflow, resulting in Inf/Inf which gives NaN. To fix this, notice that the second term is always supposed to be small, so we may write it as a small number divided by a big number if we expand the hyperbolic cosine in terms of exponentials, \cosh(x) = (e^x + e^{-x})/2, and use standard exponent rules e^ae^b =e^{a+b}
\frac{\partial}{\partial x} x \tanh(e^x) = \tanh(e^x) + \frac{x }{(e^{e^x-x/2}+e^{-e^x-x/2})/2)^2}
You could register the following custom gradient with Zygote:
julia> dtelu(x) = tanh(exp(x)) + x/((exp(exp(x)-x/2) + exp(-exp(x)-x/2))/2)^2
dtelu (generic function with 1 method)
julia> dtelu(300f0)
1.0f0
julia> dtelu(-300f0)
0.0f0