Here's a novel take on the problem and a really elegant way to see what the
"natural" chain rule is for the discrete derivative.
First, we'll assume all of our functions are polynomials. It isn't that hard to extend to analytic functions, being essentially "infinite-degree polynomials", but we'll just look at polynomials now to keep it simple.
Then the chain rule of the usual derivative is really derived from the product and sum rules. Given only those two things, we can determine the chain rule.
One really elegant way to do this is to use the algebra of dual numbers, which are of the form $a + b \epsilon$, with $\epsilon^2 = 0$. Dual numbers enable a kind of "algebraization of calculus," which is typically called "automatic differentiation." Most functions which can be extended to the complex plane can be extended to the dual plane in a similarly natural way, and in general, for any polynomial (or analytic function) $f(x)$, we have
$$
f(x + k \epsilon) = f(x) + f'(x) k \epsilon
$$
simply due to the way that dual numbers work. If you haven't tried this, play around with it, because it is almost magic. You can plug $x + \epsilon$ into basically any polynomial or Taylor series, multiply things through and cancel out the $\epsilon^2$ terms, and it will automagically compute the derivatives of huge chains of compositions for you.
The reason this works so well is that the notion of $\epsilon^2$ being zero causes the product rule to naturally flow from the properties of the algebra. For instance, assuming you have two dual numbers of the form $f(x) + f'(x) \epsilon$ and $g(x) + g'(x) \epsilon$, their product will be
$$
(f(x) + f'(x) \epsilon)(g(x) + g'(x) \epsilon) \\
= f(x)g(x) + (f'(x)g(x) + f(x)g'(x))\epsilon + (f'(x)g'(x))\epsilon^2\\
= f(x)g(x) + (f'(x)g(x) + f(x)g'(x))\epsilon
$$
So you can see that the dual part of the product is indeed the derivative of the product of the two original functions. As a result of this, it's pretty easy to verify that this property, along with the linearity of derivative addition, gives the correct "chain rule" when composing functions (and it isn't that much harder to extend to analytic functions).
The discrete derivative has a slightly different product rule, and thus a different chain rule. To see this, we'll define a "generalized discrete derivative" with step size $h$ as the following:
$$
f'_h(x) = \frac{f(x+h) - f(x)}{h}
$$
Both the forward and backward difference are instances of this derivative, as is the usual true derivative:
$$
f'_1(x) = f(x+1) - f(x) \\
f'_0(x) = f'(x) \\
f'_{-1}(x) = f(x) - f(x-1)
$$
The "correct" product rule for this, although it's a little tedious to derive, happens to be:
$$
[fg]'_h(x) = f(x)g'_h(x) + f'_h(x)g(x) + f'_h(x)g'_h(x) h
$$
which is the same as the usual product rule, but with an extra $f'_h(x)g'_h(x) h$ term.
But now we have something particularly interesting. Let's try to use this to make a new "discrete dual algebra" which can perform "automatic discrete differentiation". Instead of an element $\epsilon$, we'll have an element $\epsilon_h$. So basically, let's assume that we already have two "discrete dual numbers" of the form $f(x) + f'_h(x) \epsilon_h$ and $g(x) + g'_h(x) \epsilon_h$, and we want to use this information to compute the discrete directive of their product as the "discrete dual" part of their product. How would the algebra need to work to make this happen?
$$
(f(x) + f'(x) \epsilon_h)(g(x) + g'(x) \epsilon_h) \\
= f(x)g(x) + (f'(x)g(x) + f(x)g'(x))\epsilon_h + (f'(x)g'(x))\epsilon_h^2
$$
It turns out the magic thing we want is to declare the relation $\epsilon_h^2 = h \epsilon$, which gives the correct expression for the product rule as the discrete dual part:
$$
f(x)g(x) + (f'(x)g(x) + f(x)g'(x) + f'(x)g'(x) h)\epsilon_h
$$
and this really does work. You can go play around with this for any polynomial you want and it will always give the right result: $p(x+\epsilon_h) = p(x) + p'_h(x)\epsilon$.
And now that we have linearity and our product rule, we can compute the correct chain rule for any polynomial. This post has gotten pretty long already, and deriving this turns out to involve a bunch of tedious algebra, so I will just skip to the end. As a first step, you have this result, for any polynomial (or again, analytic) function $f(x)$
$$
f(x + k\epsilon_h) = f(x) + f'_{kh}(x) \, k \, \epsilon_h
$$
And thus, we have our generalized chain rule for any generalized discrete derivative:
$$
f(g(x + \epsilon_h)) = f(g(x) + g'_h(x) \epsilon_h) \\
= f(g(x)) + \left(f'_{g'_h(x)\cdot h}(g(x)) \cdot g'_h(x) \right) \epsilon_h
$$
which is fairly elegant. It's very similar to the original chain rule, and indeed as $h \to 0$ you get the original chain rule. The interesting thing to note here is that the step size changes for what would normally be the $f'(g(x))$ term: instead of having a step size of $h$, it now has a step size of $g'_h(x) \cdot h$. But no problem.
Lastly, it turns out that all of these "discrete dual number" algebras are really just the same algebra, which is the algebra of "split complex numbers." The split-complex numbers are isomorphic to $\Bbb R^2$ with pointwise product and sum, and to represent our "discrete dual numbers" in this algebra, we send $a + b \epsilon_h$ to $[a,a] + b[h,0] = [a+bh, a]$. (Note that the one exception is with $h=0$, where you get the dual numbers instead.)