Generalizing the sum of consecutive cubes $\sum_{k=1}^n k^3 = \Big(\sum_{k=1}^n k\Big)^2$ to other odd powers

This is a partial answer, it just establishes the existence.

We have $$s_m(n) = \sum_{k=1}^n k^m = \frac{1}{m+1} \left(\operatorname{B}_{m+1}(n+1)-\operatorname{B}_{m+1}(1)\right)$$ where $\operatorname{B}_m(x)$ denotes the monic Bernoulli polynomial of degree $m$, which has the following useful properties: $$\begin{align} \int_x^{x+1}\operatorname{B}_m(t)\,\mathrm{d}t &= x^m \quad\text{(from which everything else follows)} \\\operatorname{B}_{m+1}'(x) &= (m+1)\operatorname{B}_m(x) \\\operatorname{B}_m\left(x+\frac{1}{2}\right) &\begin{cases} \text{is even in $x$} & \text{for even $m$} \\ \text{is odd in $x$} & \text{for odd $m$} \end{cases} \\\operatorname{B}_m(0) = \operatorname{B}_m(1) &= 0 \quad\text{for odd $m\geq3$} \end{align}$$

Therefore, $$\begin{align} s_m(n) &\text{has degree $m+1$ in $n$} \\s_m(0) &= 0 \\s_m'(0) &= \operatorname{B}_m(1) = 0\quad\text{for odd $m\geq3$} \\&\quad\text{(This makes $n=0$ a double zero of $s_m(n)$ for odd $m\geq3$)} \\s_m\left(x-\frac{1}{2}\right) &\begin{cases} \text{is even in $x$} & \text{for odd $m$} \\ \text{is odd in $x$} & \text{for even $m\geq2$} \end{cases} \end{align}$$

Consider the vector space $V_m$ of univariate polynomials $\in\mathbb{Q}[x]$ with degree not exceeding $2m+2$, that are even in $x$ and have a double zero at $x=\frac{1}{2}$. Thus $V_m$ has dimension $m$ and is clearly spanned by $$\left\{s_j^2\left(x-\frac{1}{2}\right)\mid j=1,\ldots,m\right\}$$ For $m>0$, we find that $s_{2m+1}(x-\frac{1}{2})$ has all the properties required for membership in $V_m$. Substituting $x-\frac{1}{2}=n$, we conclude that there exists a representation $$s_{2m+1}(n) = \sum_{j=1}^m a_{m,j}\,s_j^2(n) \quad\text{for $m>0$ with $a_{m,j}\in\mathbb{Q}$}$$ of $s_{2m+1}(n)$ as a linear combination of squares of sums.

Adding to @Hecke's answer
First step was to build the matrix $A$ of coefficients of the squares of the Faulhaber-polynomials. Next step, to build the matrix $B$ of the coefficients of the un-squared Faulhaber-polynomials, but only of the odd orders.

Then the matrix $C$ of coefficients for the composition, such that $C \cdot A = B$ occurs by inversion of $A$ such that $C = B \cdot A^{-1}$. This can be done with any truncation of size. Important: Note that the first two columns of $A$ and $B$ are completely zero, so we must reduce that matrices to the remaining columns.

The inverse of the coefficients-matrix $C$ (which itself can be seen in the inner braces of Hecke's answer) might also be interesting; at least it has a slightly simpler structure: $$ C^{-1}= \tiny \begin{bmatrix} 2 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1/3 & 2/3 & . & . & . & . & . & . & . & . & . & . & . \\ . & . & 1/2 & 1/2 & . & . & . & . & . & . & . & . & . & . \\ . & . & -1/15 & 2/3 & 2/5 & . & . & . & . & . & . & . & . & . \\ . & . & . & -1/6 & 5/6 & 1/3 & . & . & . & . & . & . & . & . \\ . & . & . & 1/21 & -1/3 & 1 & 2/7 & . & . & . & . & . & . & . \\ . & . & . & . & 1/6 & -7/12 & 7/6 & 1/4 & . & . & . & . & . & . \\ . & . & . & . & -1/15 & 4/9 & -14/15 & 4/3 & 2/9 & . & . & . & . & . \\ . & . & . & . & . & -3/10 & 1 & -7/5 & 3/2 & 1/5 & . & . & . & . \\ . & . & . & . & . & 5/33 & -1 & 2 & -2 & 5/3 & 2/11 & . & . & . \\ . & . & . & . & . & . & 5/6 & -11/4 & 11/3 & -11/4 & 11/6 & 1/6 & . & . \\ . & . & . & . & . & . & -691/1365 & 10/3 & -33/5 & 44/7 & -11/3 & 2 & 2/13 & . \\ . & . & . & . & . & . & . & -691/210 & 65/6 & -143/10 & 143/14 & -143/30 & 13/6 & 1/7 \end{bmatrix}$$ ... and the numerator of $691$ in rows $13,14,15,\ldots$ surely rings some bell concerning Bernoulli-numbers
Note: this matrix gives the compositions of the squared sums-of-powers by the higher odd orders of the original sums-of-powers, so helps for the inverse of your problem: $$ \begin{array}{rrrrr} S(0,x)^2 &=& 2 S(1,x) & - 1 S(0,x) &\\ \hline \\ S(1,x)^2 &=& 1 S(3,x) \\ S(2,x)^2 &=& 1/3 \; S(3,x) &+2/3 \; S(5,x) \\ S(3,x)^2 &=&& 1/2 \; S(5,x) &+1/2 \; S(7,x) \\ S(4,x)^2 &=&& -1/15 \; S(5,x)&+2/3 \; S(7,x)&+2/5 \; S(9,x) \\ \vdots \end{array} $$ The first row does not fit into the pattern. The coefficients are compositions of binomials and Bernoulli-numbers, which might best be written as below according to the analysis of the matrix $E$

[update] The following matrix displays coefficients after the Bernoulli-numbers and the column-index are removed. $$ D = \Tiny \begin{bmatrix} 2 & . & . & . & . & . & . & . & . & . & . & . \\ . & 2 & . & . & . & . & . & . & . & . & . & . \\ . & 4 & 2 & . & . & . & . & . & . & . & . & . \\ . & . & 9 & 2 & . & . & . & . & . & . & . & . \\ . & . & 6 & 16 & 2 & . & . & . & . & . & . & . \\ . & . & . & 20 & 25 & 2 & . & . & . & . & . & . \\ . & . & . & 8 & 50 & 36 & 2 & . & . & . & . & . \\ . & . & . & . & 35 & 105 & 49 & 2 & . & . & . & . \\ . & . & . & . & 10 & 112 & 196 & 64 & 2 & . & . & . \\ . & . & . & . & . & 54 & 294 & 336 & 81 & 2 & . & . \\ . & . & . & . & . & 12 & 210 & 672 & 540 & 100 & 2 & . \\ . & . & . & . & . & . & 77 & 660 & 1386 & 825 & 121 & 2 \end{bmatrix} $$ We have $$ C^{-1}_{r,c} = D_{r,c}\cdot {B_{2(r-c)} \over c } $$ where $B_k$ are the Bernoulli-numbers and the column- and rowindices $c,r$ begin at $1$.

Perhaps an even better display is this: $$ E= \Tiny \begin{bmatrix} 1 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\ . & 3 & 1 & . & . & . & . & . & . & . & . & . & . & . \\ . & . & 6 & 1 & . & . & . & . & . & . & . & . & . & . \\ . & . & 5 & 10 & 1 & . & . & . & . & . & . & . & . & . \\ . & . & . & 15 & 15 & 1 & . & . & . & . & . & . & . & . \\ . & . & . & 7 & 35 & 21 & 1 & . & . & . & . & . & . & . \\ . & . & . & . & 28 & 70 & 28 & 1 & . & . & . & . & . & . \\ . & . & . & . & 9 & 84 & 126 & 36 & 1 & . & . & . & . & . \\ . & . & . & . & . & 45 & 210 & 210 & 45 & 1 & . & . & . & . \\ . & . & . & . & . & 11 & 165 & 462 & 330 & 55 & 1 & . & . & . \\ . & . & . & . & . & . & 66 & 495 & 924 & 495 & 66 & 1 & . & . \\ . & . & . & . & . & . & 13 & 286 & 1287 & 1716 & 715 & 78 & 1 & . \\ . & . & . & . & . & . & . & 91 & 1001 & 3003 & 3003 & 1001 & 91 & 1 \end{bmatrix}$$ where $ C^{-1}_{r,c} = E_{r,c} \cdot 2 B_{2(r-c)} / r $ and it seems $E_{r,c} = \binom{r}{2(r-c)} $ where $ c \le r \lt 2c$ so it seems $$ C^{-1}_{r,c} = B_{2(r-c)} \cdot \binom{r}{2(r-c)}\cdot \frac 2{2c-r+1} \qquad \text{ where } c \le r \lt 2c \tag 1$$ [added]
If I use this description of the composition of $C^{-1}$ then I can symbolically invert this again and get (the top-left range here) $$ C= \Tiny \begin{bmatrix} 1/2 & . & . & . & . & . \\ 0 & 1 & . & . & . & . \\ 0 & -3 \; B_2 & 3/2 & . & . & . \\ 0 & 18 \; B_2^2 & -9 \; B_2 & 2 & . & . \\ 0 & -180 \; B_2^3+15 \; B_4 \; B_2 & 90 \; B_2^2-15/2 \; B_4 & -20 \; B_2 & 5/2 & . \\ 0 & 2700 \; B_2^4-495 \; B_4 \; B_2^2 & -1350 \; B_2^3+495/2 \; B_4 \; B_2 & 300 \; B_2^2-30 \; B_4 & -75/2 \; B_2 & 3 \end{bmatrix} \tag 2$$ and where the first row/column should be deleted since they are out of the pattern. For example, using the third and the fourth row in this matrix we get $$ \begin{array} {rrlll} S(5,x) &=& -3 \; B_2 \; S(1,x)^2 &+ 3/2 \; S(2,x)^2 \\ S(7,x) &=& 18 \; B_2^2 \; S(1,x)^2 &- 9 \; B_2 \; S(2,x)^2 &+ 2 \; S(3,x)^2 \end{array} \tag 3$$

Not an answer, but only some results I found.

Denote $s_m=\sum_{k=1}^nk^m.$ Suppose that $$s_{2m+1}=\sum_{i=1}^m a_is_i^2,a_i\in\mathbb Q.$$

Use the method of undetermined coefficients, we can get a list of $\{m,\{a_i\}\}$: $$\begin{array}{l} \{1,\{1\}\} \\ \left\{2,\left\{-\frac{1}{2},\frac{3}{2}\right\}\right\} \\ \left\{3,\left\{\frac{1}{2},-\frac{3}{2},2\right\}\right\} \\ \left\{4,\left\{-\frac{11}{12},\frac{11}{4},-\frac{10}{3},\frac{5}{2}\right\}\right\} \\ \left\{5,\left\{\frac{61}{24},-\frac{61}{8},\frac{28}{3},-\frac{25}{4},3\right\}\right\} \\ \left\{6,\left\{-\frac{1447}{144},\frac{1447}{48},-\frac{332}{9},\frac{595}{24},-\frac{21}{2},\frac{7}{2}\right\}\right\} \\ \left\{7,\left\{\frac{5771}{108},-\frac{5771}{36},\frac{5296}{27},-\frac{2375}{18},56,-\frac{49}{3},4\right\}\right\} \\ \left\{8,\left\{-\frac{53017}{144},\frac{53017}{48},-\frac{60817}{45},\frac{21817}{24},-\frac{3861}{10},\frac{1127}{10},-24,\frac{9}{2}\right\}\right\} \\ \left\{9,\left\{\frac{2755645}{864},-\frac{2755645}{288},\frac{632213}{54},-\frac{1133965}{144},\frac{13379}{4},-\frac{11725}{12},208,-\frac{135}{4},5\right\}\right\}\end{array}$$

Through some observation, we can find that $a_i a_{i+1}<0,a_m=\dfrac{m+1}2,a_2=-3a_1,a_{m-1}=\dfrac{1}{24} m^2 (m+1).$