Inductive definition of Bernstein polynomials

For $$n\in \mathbb{N}$$ let $$B_n$$ be the linear operator taking a function $$f$$ on the unit interval $$I=[0,1]$$ to its $$n$$-th Bernstein polynomial $$B_nf$$, $$B_nf(x):=\sum_{k=0}^n\binom{n}{k} f\Big(\frac{k}{n}\Big)x^k(1-x)^{n-k}\label{1}\tag{1}$$ The polynomial $$B_nf(x)$$ has a natural probabilistic interpretation, namely, it is the expected value of $$f(\xi)$$, where $$\xi=\frac{1}{n}\sum_{j=1}^n \omega_j$$ is the average of $$n$$ independent random variables with identical Bernoulli distribution of parameter $$x$$, that is, $$\mathbb{P}(\omega_j=1)=x$$. In fact, this is the starting point in the beautiful Bernstein's proof of the Weierstrass' density theorem via the WLLN. However, this question is about an alternative definition of the sequence $$(B_n)_{n\ge0}$$.

Let $$D:C^1(I)\to C^0(I)$$ be the derivative operator, and for all $$n\ge1$$, let $$D_n:C^0(I)\to C^0(I)$$ be the approximate discrete derivative given by the incremental ratio $$D_nf(x):=\frac{f\big( \frac{n-1}{n} x+\frac{1}{n}\big)-f\big( \frac{n-1}{n} x\big)}{\frac{1}{n}},$$ (which is well-defined for $$f\in C^0(I)$$ and $$x\in I$$).

It is easy to check that definition \eqref{1} implies $$DB_n=B_{n-1}D_n\label{2}\tag{2}$$ together with:
$$B_0f(x)=B_nf(0)=f(0)\label{3}\tag{3}$$
Conversely these two imply formula \eqref{1}, as it follows immediately by induction, at least, if we already have it (quite a common situation of formulas proven by induction). Thus, since \eqref{2} and \eqref{3} characterize $$(B_n)_n$$, we may even take them as an inductive definition of $$(B_n)_n$$. Note that replacing $$D_n$$ with $$D$$ in \eqref{2} gives the analogous inductive definition for the Taylor polynomials in $$0$$. (Incidentally, formula \eqref{2} is relevant in the approximation theory, in that it implies that for $$f\in C^k(I)$$ one has $$B_nf\to f$$ in $$C^k$$: this by induction from the case $$k=0$$, since $$D_n$$ converges strongly to $$D$$. Also, it says that if some derivative $$f^{(k)}$$ is non-negative on $$I$$, so is $$(B_nf)^{(k)}$$.)

Question: How can we deduce naturally formula \eqref{1} (i.e., assuming we don't know it, and we do not have a crystal ball to guess it) from \eqref{2} and \eqref{3}?

$$\newcommand{\De}{\Delta}$$ Iterating your condition \eqref{2}, for $$k=0,\dots,n$$ we have $$\begin{equation*} D^kB_n=\frac{n!}{(n-k)!}\,B_{n-k}P_{n,k},\label{a}\tag{a} \end{equation*}$$ where $$\begin{equation*} P_{n,k}:=\De_{n-k+1}\cdots\De_n,\quad \De_j:=\tfrac1j\,D_j. \end{equation*}$$ By induction on $$k=0,\dots,n$$, $$\begin{equation*} (P_{n,k}f)(x)=\sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac{n-k}n\,x+\frac in\Big),\label{b}\tag{b} \end{equation*}$$ whence, using \eqref{a} and taking your condition \eqref{3} into account, we have $$\begin{equation*} \frac{(n-k)!}{n!}\,(D^kB_n f)(0)=(B_{n-k}P_{n,k}f)(0)=(P_{n,k}f)(0) =\sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac in\Big). \end{equation*}$$ Also, using again \eqref{a} and \eqref{b}, and again taking your condition \eqref{3} into account, we have $$\begin{equation*} \frac1{n!}\,(D^nB_n f)(x)=(B_0P_{n,n} f)(x)=(P_{n,n} f)(0) =\sum_{i=0}^n(-1)^{n-i}\binom ni f\Big(\frac in\Big), \end{equation*}$$ a constant. So, $$B_n f$$ is a polynomial of degree $$\le n$$, and hence \begin{align*} (B_n f)(x)&=\sum_{k=0}^n \frac{(D^kB_n f)(0)}{k!}\,x^k \\ &=\sum_{k=0}^n\binom nk x^k \sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac in\Big) \\ &=\sum_{i=0}^n f\Big(\frac in\Big)\sum_{k=i}^n (-1)^{k-i}\binom nk \binom ki x^k \\ &=\sum_{i=0}^n f\Big(\frac in\Big)\binom ni x^i(1-x)^{n-i}, \end{align*} as desired.
A comment on Josif Pinelis' formula $$(b)$$ for $$\Delta_{n-k+1} \dots\Delta_{n-1}\Delta_{n}$$, which is a main point of the computation. Let $$\{\tau_{a}\}_{a\in\mathbb{R}}$$ and $$\{\delta_{b}\}_{a\in\mathbb{R}_+}$$ denote respectively the linear group of translations on functions (that we may think defined on the whole real line w.l.o.g.), $$f(\cdot)\mapsto f(\cdot+a)$$, and the linear group of dilations, $$f(\cdot)\mapsto f(\cdot b)$$. So $$\tau_{a+b}=\tau_a\tau_b,$$ $$\delta_{ab}=\delta_a\delta_b,$$ $$\tau_{ab}=\delta_b^{-1}\tau_a\delta_b$$ Since $$\Delta_n:=\delta_{\frac{n-1}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)$$, moving all dilations on the left by the above relations imply nicely $$\Delta_{n-k+1} \dots\Delta_{n-1}\Delta_{n}=\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k},$$ whence $$\frac{1}{k!} D^kB_n=\frac{1}{k!}B_{n-k} D _{n-k+1} \dots D _{n-1} D _{n}=\Big({n\atop k}\Big)B_{n-k}\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k},$$ which we can expand to formula $$(b)$$.
edit. In fact we may skip the last expansion too, keeping all Josif's formulas on the level of operators. Since the $$D_k$$'s lower the degree of polynomials, $$(2)$$ and $$(3)$$ imply that $$B_n$$ takes values on polynomials of degree less than or equal to $$n$$, as said. So, for any $$x$$, denoting $$e_x$$ the evaluation form, $$e_xB_n=e_0\bigg[\sum_{k=0}^n \frac{x^k}{k!}D^kB_n\bigg]=e_0\bigg[\sum_{k=0}^n x^k \Big({n\atop k}\Big)B_{n-k}\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k}\bigg]=$$ $$=e_0\bigg[\sum_{k=0}^n \Big({n\atop k}\Big)x^k\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k}\bigg]=e_0\bigg( \mathbb{1} + x \big(\tau_{\frac{1}{n}}-\mathbb{1}\big) \bigg)^n =e_0\bigg( x \tau_{\frac{1}{n}} + (1-x)\mathbb{1} \bigg)^n$$ $$=e_0\bigg(\sum_{k=0}^n \Big({n\atop k}\Big)x^k(1-x)^{n-k}\tau_{\frac{k}{n}} \bigg)$$ which indeed takes $$f$$ to the original $$(B_nf)(x)$$ given by $$(1)$$.