I'm sure there may be a more compact answer, but that would require extensive time. The transform is realized by first realizing the identity between associated Legendre polynomials and Hypergeometric functions. Then converting the Hypergeometric function to its summation formula provides the requisite simplification to enable an inverse Mellin transform. See Mathematica output, attached.

This might give you some better ideas.

Very clever.

Is there any way to preserve the Legendre polynomial in the final expression (instead of the hypergeometric function)?

Unfortunately, I don't think there's any easy means to do so... I was intrigued by your problem, but from experience, special functions rarely enable easy representation. If you think about it, what you call Legendre polynomials, they are themselves a class of polynomials, actually, an infinite summation; but, we assign a nice name to it and-viola-it seems like a simple "function."

Personally, the minute I see hypergeometric functions, Meijer G-functions, &c, I consider the solution "out-of-hand;" because, these representations are not compact--generally--do not converge very fast, do not allow easy manipulation, obscure the underlying physics and behavior, etc.

Do remember that most of these special functions arise from transforming a function from some domain punctured with poles and branches, hence, the unavoidable infinite sums, power series, etc. Thus, the solution is truly a series solution, inept, a polynomial solution. Although, if the solution were a Green's function, that too could be made into a series solution, but its representation is very compact, complete. The problem I have in general with solutions containing special functions is the radius of convergence is quite small and further complications become unwieldy, in fact, they can easily lead one astray.

Alright, so, my first stab at the problem was a bit lazy...so be it!

Here are some other routes of attack:

1] use the definition of a generating function for Legendre polynomials and perform inverse Mellin transform on that.

2] use fact that a Mellin transform is related to a Laplace transform, where one might find the inverse Laplace transform of Legendre polynomials somewhere or by hand.

3] use fact that Mellin transforms are related to Fourier transforms...hopefully find inverse, etc, etc.

4] use convolution theorem for Mellin transforms, therefore, break up problem into known pieces, transform each piece, then hope you can find the convolutions of all the pieces.

Now, another approach is to find something like I did and manipulate the formulae to force out certain definitions and forms so desired, for example, if you multiply and divide by k factorial, I think you could think about coaxing the representations back into Legendre polynomials, only with k! as a multiplicant.

Sometimes, you can get your cake but you can't eat it to, so, think about asymptotics, expansions, etc; these solutions often provide the salient information sought and enable you to further investigate the physics.

I'm gonna give you two important books. The second will be of particular interest to you, because it lists Mellin transforms.

1] https://drive.google.com/file/d/1jt7zDKd3ao1fSgVtBtwYoCWQLTVWw4Y7/view?usp=sharing

2] https://drive.google.com/file/d/1yowf9EcoG5pJOLdWdhRU9LL5gLzCR8_e/view?usp=sharing

I will tell you that these tables of integrals rarely have the "one" I'm looking for at the time, so, I usually get stuck, in the end, working the problem out by hand.

Your problem is an intriguing one and it seems the material you are working on is equally interesting. I don't want to just give the answer--not that I am implying I have the answer--rather, I find particular joy in giving the tools, suggestions and allowing you to claim victory. You should own the answer. As it stands, the variety of integrals are simply infinite, hence, you find over time building a set of functional transforms particular to your field of research and focus. These tables could one day fill a book like the two linked above.

Remember the convolution theorem, this can be a real life safer in these types of problems.

Here's a bit more:

Article Mellin transforms with only critical zeros: Legendre functions

...gets thick, eh?

So, here is the inverse Mellin transform for the generating function that represents Legendre's polynomials, see attached.

Here's another book I think you might like:

https://drive.google.com/file/d/1mduTWY_mtCgYCiPin-9dERhKmVUH4XeN/view?usp=sharing

From Fritz Oberhettinger's book, page 79, entry 8.45, just might be the ticket!

Wolfram can be a great resource.

https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1Regularized/26/01/01/0002/

I think this would relate the hypergeometric function PFQ regularized to Legendre ploynomial.

This is my guess. Another option is to look for Bivariate Meijer-G function, which I believe will eliminate the infinite sum. Also, the MATLAB code for bivariate Meijer G-function is available on the Internet (you can just google it).