To answer your question, a google search of part of the formula shows no other results but this one, suggesting it is original
https://www.google.com/search?q=%22376698240%22&rlz=1C1ONGR_enUS980US980&oq=%22376698240https://www.google.com/search?q=%22376698240%22
It' possible to make arbitrary series by just coming up with bigger and bigger polynomials p(n) and r(n). Such polynomials arise by starting with a kernel, which is the original integral that yields the desired constant, and then using polynomial division on the integral to make new binomial identities. It's possible to make an arbitrary number of these, but the most impressive series are those that have a fast convergence and small polynomials . For Zeta(3), you would use three integrals for the underlying 4F3 hypergeometric function. This is done with Catalan's Constant on page 13 here, which requires two integrals. https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3919892
The more interesting question is how to derive these. Just saying "WZ algorithm" is not that helpful. There needs to be a mathematically precise way to finding the fastest series ofconstrained by a certain hypergeometric form. For example, the fastest non-reciprocal pi formula (if we limit to reciprocal of pi, Ramanujan-type formulas would obviously win) that can be expressed as the composition of two 3f2 series in which all the parameters are rational.