About a chain rule for Wronskians
Computing the cells of the matrix row by row, starting at the top, for the first three rows we have with $y=g(x)$ for brevity, the matrix $M(f_1\circ g,...,f_n\circ g)$ involved in the Wronskian is $$ M(f_1, \ldots, f_n)(g(x)) = \begin{pmatrix} f_1(y) & f_2(y) & \cdots & f_n(y) \\ f_1'(y)g'(x) & f_2'(y)g'(x) & \cdots & f_n'(y)g'(x) \\ f_1^{(2)}(y)g'^2(x) & f_2^{(2)}(y)g'^2(x) & \cdots & f_n^{(2)}(y)g'^2(x) \\\vdots & \vdots & \vdots & \vdots \end{pmatrix} + \begin{pmatrix} f_1(y) & f_2(y) & \cdots & f_n(y) \\ f_1'(y)g'(x) & f_2'(y)g'(x) & \cdots & f_n'(y)g'(x) \\ f_1'(y)g''(x) & f_2'(y)g''(x) & \cdots & f_n'(y)g''(x) \\\vdots & \vdots & \vdots & \vdots \end{pmatrix}. $$ The second and third rows are linearly dependent (e.g., take the weighted sum of the second row weighted with $g''(x)$ and the third weighted with $-g'(x)$). So for the first three rows, by properties of the determinant we can remove the second term above without changing the value of the determinant. Only the first term survives, which corresponds to the formula you mention, with powers of $g'(x)$ but no higher derivatives of $g$.
For the next rows (4th, 5th etc), say the $k$-th row, only the term in $g'(x)^{k-1}$ survives as all others lead to rows that are linearly dependent with the previous rows (that we have "cleaned" one by one, as for the third row above).
The rest is a combinatorial problem. The second row has $g'(x)$, third row has $g'^2(x)$, fourth row has $g'^3(x)$ etc. As multiplying a row by $C$ multiplies the determinant by $C$, we find a total of $g'(x)^{1+2+...+n-1}$ that can be factored out of the determinant. Since $1+2+...+n-1=n(n-1)/2$, I very much agree with your nice formula!
What an interesting result! I did find this result as an exercise on pg. 217 of A short course in the theory of determinants.
In searching for a reference, I wanted to note similar product and quotient rules that arise from the same approach:
$$W[f_1\cdot g,\ldots,f_n\cdot g] = g^n\; W[f_1,\ldots,f_n]$$ $$W[f_1/g,\ldots,f_n/g] = g^{-n}\; W[f_1,\ldots,f_n]$$
In this setting, the combinatorics come from Pascal's triangle (binomials) instead of the Bell polynomials.
Perhaps this homogeneity property of the Wronskian will help track down the result.
The earliest reference I could find for this identity is a paper of Hurwitz from 1892 titled Über algebraische Gebilde mit eindeutigen Transformationen in sich, which can be found here. Here's a screenshot of the Wronskian identity appearing on page 407 of the paper:
Note: I'll try to translate this if I have some time, but my German is very poor. If anybody has some free time and can do better than Google Translate, they are very welcome to do so!
In formula (2), the expressions $\Delta_u$ and $\Delta_t$ are the Wronskian of holomorphic functions $u_1, \ldots, u_p$ expressed in two different complex coordinates $u$ and $t$ on a Riemann surface. The coordinates are related by a change-of-coordinates map $t = g(u)$ and the notation $\frac{dt}{du}$ means $g'(u)$.
Notice that this formula is not quite the same as the one you wrote: the exponent is $n(n+1)/2$ instead of $n(n-1)/2$. Notice also that the definition of the Wronskian seems odd: the first row has the derivatives of the functions instead of the functions themselves!
This is because Hurwitz is not taking the Wronskian of the functions $u_i$, but rather of the one-forms $du_i = \frac{du_i}{dz} dz$. By definition, the Wronskian of one-forms is obtained by plugging their coordinate functions inside the usual Wronskian, whence the formula. This also explains the different exponent: one-forms do not transform the same way as functions under changes of coordinates and so the Wronskian picks up an extra factor of the derivative when changing coordinates.
Here is some context for the appearance of the Wronskian in Hurwitz's work.
The Wronskian appears in the theory of Riemann surfaces because it is closely related to the Weierstrass points of a Riemann surface. These are finitely many points on the surface which are very special because they are the only points $p$ for which there exist non-constant meromorphic functions with a pole of order at most $g$ at $p$ and no other poles. Here $g$ is the genus of the surface.
The "chain rule" that you wrote means that the Wronskian of $n$ holomorphic functions defines a "$N$-differential", where $N=n(n-1)/2$. A $N$-differential is a section of the $N$-th tensor power of the holomorphic cotangent bundle of the surface. This is just a fancy way of saying that it's a collection of functions defined on coordinate charts that transform like $W(z) = g'(z)^N \tilde{W}(g(z))$ under changes of coordinates. The $N$-differential is then the global tensor obtained by patching together these local functions.
There are no non-constant holomorphic functions on a compact Riemann surface, so one is lead to consider the Wronskian of holomorphic $1$-forms instead of functions. This is defined locally by taking the Wronskian of the coordinate functions of the forms and again this defines a $N$-differential, where now $N=n(n+1)/2$ for the reasons explained above. For a surface of genus $g \geq 1$, the Wronskian of a basis of holomorphic $1$-forms vanishes only on a finite set of points, and these turn out to be exactly the Weierstrass points of the surface.
There are a lot of modern books on Riemann surfaces that cover this story. For example, it is discussed in section 4 of chapter 7 of Miranda's book Algebraic curves and Riemann surfaces. The Wronskian identity appears in the proof of Lemma 4.9 on page 237. This is also covered on pages 82-86 of the book Riemann surfaces by Farkas and Kra.