We study the equation $M_Ψ(z+1)=\frac{-z}{Ψ(-z)}M_Ψ(z), M_Ψ(1)=1$ defined on a subset of the imaginary line and where $Ψ$ is a negative definite functions. Using the Wiener-Hopf method we solve this equation in a two terms product which consists of functions that extend the classical gamma function. These functions are in a bijection with Bernstein functions and for this reason we call them Bernstein-gamma functions. Via a couple of computable parameters we characterize of these functions as meromorphic functions on a complex strip. We also establish explicit and universal Stirling type asymptotic in terms of the constituting Bernstein function. The decay of $|M_Ψ(z)|$ along imaginary lines is computed. Important quantities for theoretical and applied studies are rendered accessible. As an application we investigate the exponential functionals of Levy Processes whose Mellin transform satisfies the recurrent equation above. Although these variables have been intensively studied, our new perspective, based on a combination of probabilistic and complex analytical techniques, enables us to derive comprehensive and substantial properties and strengthen several results on the law of these random variables. These include smoothness, regularity and analytical properties, large and small asymptotic behaviour, including asymptotic expansions, bounds, and Mellin-Barnes representations for the density and its successive derivatives. We also study the weak convergence of exponential functionals on a finite time horizon when the latter expands to infinity. As a result of new factorizations of the law of the exponential functional we deliver important intertwining relation between members of the class of positive self-similar semigroups. The derivation of our results relies on a mixture of complex-analytical and probabilistic techniques.