  We consider the chemotaxis-fluid system
\begin{align}\label{star}\tag{$\diamondsuit$} \left\{
\begin{array}{r@{\,}c@{\,}c@{\ }l@{\quad}l@{\quad}l@{\,}c}
n_{t}&+&u\cdot\!\nabla n&=\Delta n^m-\nabla\!\cdot(n\nabla c),\ &x\in\Omega,&
t>0,\\ c_{t}&+&u\cdot\!\nabla c&=\Delta c-c+n,\ &x\in\Omega,& t>0,\\
u_{t}&+&(u\cdot\nabla)u&=\Delta u+\nabla P+n\nabla\phi,\ &x\in\Omega,& t>0,\\
&&\nabla\cdot u&=0,\ &x\in\Omega,& t>0, \end{array}\right. \end{align} in a
bounded domain $\Omega\subset\mathbb{R}^3$ with smooth boundary and $m>1$.
Assuming $m>\frac{4}{3}$ and sufficiently regular nonnegative initial data, we
ensure the existence of global solutions to the no-flux-Dirichlet boundary
value problem for \eqref{star} under a suitable notion of very weak
solvability, which in different variations has been utilized in the literature
before. Comparing this with known results for the fluid-free setting of
\eqref{star} the condition appears to be optimal with respect to global
existence. In case of the stronger assumption $m>\frac{5}{3}$ we moreover
establish the existence of at least one global weak solution in the standard
sense.
  In our analysis we investigate a functional of the form $\int_{\Omega}\!
n^{m-1}+\int_{\Omega}\! c^2$ to obtain a spatio-temporal $L^2$ estimate on
$\nabla n^{m-1}$, which will be the starting point in deriving a series of
compactness properties for a suitably regularized version of \eqref{star}. As
the regularity information obtainable from these compactness results vary
depending on the size of $m$, we will find that taking $m>\frac{5}{3}$ will
yield sufficient regularity to pass to the limit in the integrals appearing in
the weak formulation, while for $m>\frac{4}{3}$ we have to rely on milder
regularity requirements making only very weak solutions attainable.
