  We construct and study a certain zeta function which interpolates
multi-poly-Bernoulli numbers at non-positive integers and whose values at
positive integers are linear combinations of multiple zeta values. This
function can be regarded as the one to be paired up with the $\xi$-function
defined by Arakawa and the first-named author. We show that both are closely
related to the multiple zeta functions. Further we define multi-indexed
poly-Bernoulli numbers, and generalize the duality formulas for poly-Bernoulli
numbers by introducing more general zeta functions.
