In section 1-5 of the third edition ofFoundations of Electromagnetic Theoryby Reitz, Milford and Christy, the authors give a coordinate-system-independent definition of the divergence of a vector field:

$$\nabla\cdot\mathbf{F} = \lim_{V\rightarrow 0}\frac{1}{V}\int_S\mathbf{F\cdot n}da$$

with $\displaystyle S$ the surface enclosing the volume $\displaystyle V$.

In deriving the Cartesian form of the divergence they use $\displaystyle \Delta x\Delta y\Delta z$ as a volume element over which to make what appear to be first order Taylor expansions:

$$F_x(x_0 + \Delta x, y, z) = F_x(x_0, y, z) + \Delta x\frac{\partial F_x}{\partial x}\Bigg\arrowvert_{x_0, y, z}$$

$$F_y(x, y_0 + \Delta y, z) = F_y(x, y_0, z) + \Delta y\frac{\partial F_y}{\partial y}\Bigg\arrowvert_{x, y_0, z}$$

$$F_z(x, y, z_0 + \Delta z) = F_z(x, y, z_0) + \Delta z\frac{\partial F_z}{\partial z}\Bigg\arrowvert_{x, y, z_0}$$

I say "appear to be" because the zero order term and all derivatives in a Taylor expansion are evaluated at the initial point $\displaystyle (x_0, y_0, z_0)$. In the next step they plug the above and $\displaystyle \mathbf{F}(x_0, y_0, z_0)$ into the divergence equation, with the six sides of $\displaystyle \Delta x\Delta y\Delta z$ as the surface:

$$\nabla\cdot\mathbf{F} = \lim_{V\rightarrow 0}\frac{1}{\Delta x\Delta y\Delta z}\Bigg\{\int F_x(x_0, y, z)dydz + \Delta x\Delta y\Delta z\frac{\partial F_x}{\partial x}\\ + \int F_y(x, y_0, z)dxdz + \Delta x\Delta y\Delta z\frac{\partial F_y}{\partial y}\\ + \int F_z(x, y, z_0)dxdy + \Delta x\Delta y\Delta z\frac{\partial F_z}{\partial z}\\ - \int F_x(x_0, y, z)dydz - \int F_y(x, y_0, z)dxdz - \int F_z(x, y, z_0)dxdy\Bigg\}$$

Is there anything to leaving terms in the form $\displaystyle F_x(x_0, y, z)$ instead of $\displaystyle F_x(x_0, y_0, z_0)$ and $\displaystyle \int F_x(x_0, y, z)dydz$ instead of $\displaystyle F_x(x_0, y_0, z_0)\Delta y\Delta z$, especially since the partials are treated like constants which get multiplied by the areas of the sides outright?