(Revamped. Majority of entry at my mini-arxiv now.)

The natural, unique extension of the Bernoulli numbers is the Bernoulli Appell polynomial sequence, which can be operationally defined by action on polynomials and analytic functions, such as the exponential and logarithm, when convergent, or order by order for a formal power series, through the umbral relation $f(B.(x+1))-f(B.(x))={f}'(x)=D_x \; f(x)$, where D is the derivative. Action on $exp(xt)$ gives the e.g.f., of the Bernoulli polynomials without any resort to numerical values of the Bernoulli numbers. Action on the Mellin transform of $exp(-xt)$ defines the Bernoulli numbers in terms of the Riemann zeta values. The e.g.f. of their umbral compositional inverses is the reciprocal of the one for the Bernoullis, which gives the "reciprocal polynomials", based on reciprocal natural numbers, very naturally associated with not only the exponential divided by its argument, but also the logarithm.

From grafting the Bernoulli and reciprocal polynomials together stem a pair of Lie operator derivatives for powers of the state number, or Euler op, and associated normal ordered ops. The Lie ops are related to the compositional inverse pair of functions that are the log and exp functions for the multiplicative formal group law associated to the Todd class. The matrix reps are conjugates of the infinitesimal generator of the Pascal triangular matrix by the mutually orthogonal Stirling number matrices and encode the combinatorics of simplices. The multiplicative, compositional, and umbral compositional inversions are inextricably bound together and reveal myriad associations to combinatorics, Lie theory, and topology.

The interplay of the Bernoulli and reciprocal polynomials reveal this. It also provides easy proofs of many, if not most, identities involving the Bernoullis and a way of looking at the Riemann zeta function that can not be readily achieved from the perspective of the e.g.f. operators of the Euler-Maclaurin expansion. For example, regarding the Mellin transform as a means of interpolation, the natural extension of the Bernoulli polynomials is the Hurwtz zeta function.

$$B_{-s}(x)=s \sum_{n \ge 0}\frac{1}{(n+x)^{s+1}},$$ which with $x=1$ becomes $s\cdot \zeta(s+1)$, and for the reciprocal integers,

$$\bar{B}_{-s}(x)=\frac{(x+1)^{1-s}-x^{1-s}}{1-s}.$$ The two are related through umbral composition and inversion, so that the pole singularities are reflected in each other and, in fact, the Gauss-Newton series and umbral composition leads to

$$\zeta(s)=\sum_{n \ge 0}(-1)^{n+1}\;\frac{(-s)!}{n!(2-s-n)!} \frac{2^{2-s}-2^n}{2^n}C_n,$$

for $Re(s)<1$ (but gives good results with just eight terms over the range of reals $ -6 \le s \le 2$--it's capturing the dependence of zeta on the singularity, the falling factorials of $s$, and zeta's first three simple zeroes up to that approx.--ten terms captures the dependence on the next zero) where $C.=(1,1,5/6,1/2,1/10,-1/6,-5/42,1/6,...)$ are determined by $C_n=(1-G.)^n$ and $G.=(1,0,-1/6,0,1/10,0,-5/42,0,...)$ come from the umbral composition of the Bernoulli polynomials with the Bernoulli numbers $ B_n(1)=(-1)^nB_n(0)$.

For relations to Hirzebruch genera / Todd class, see this MOQ.