Manually fetching 8 neighboring texels from 3D texture greatly decreases performance
What's probably happening is that you're doing indirect texture fetches. Are you doing something like this:

vec4 a = texture3D(sampler, coord);
vec4 b = texture3D(sampler, coord + vec4(0, 0, 1));
vec4 c = texture3D(sampler, coord + vec4(0, 1, 0));
vec4 d = texture4D(sampler, coord + vec4(1, 0, 0));
//...etc.


If so, that can slow you down. What you can do instead is to pass the texture coordinates for each sampler call into your vertex shader from your CPU code. You'd set it up as if you were doing multi-texturing, but all the texture units would point to the same texture, just sampled at different coordinates. 

Taking the 2D case, instead of passing in a single texture coordinate like this:

float2 texCoords[] = { 
    { 0, 0 },
    { width, 0 },
    { width, height },
    { 0, height }
};


You'd pass in 4 of them:

float2 texCoords[] = {
    // texture 0 coords
    { 0, 0 },
    { width, 0 },
    { width, height },
    { 0, height },

    // texture 1 coords
    { 1, 0 },
    { width + 1, 0 },
    { width + 1, height },
    { 1, height },

    // Texture 2 coords
    { 0, 1 },
    { width, 1 },
    { width, height + 1 },
    { 0, height + 1 },

    // Texture 3 coords
    { 1, 1 },
    { width + 1, 1 },
    { width + 1, height + 1 },
    { 1, height + 1 }
};


In your fragment shader, you'd do something like this:

vec4 pixel0 = texture2D(sampler0, coord0);
vec4 pixel1 = texture2D(sampler1, coord1);
vec4 pixel2 = texture2D(sampler2, coord2);
vec4 pixel3 = texture2D(sampler3, coord3);


Now you have the 4 pixels needed for interpolation, but there was no fragment shader math done to get the texture coords, so there are no indirect look-ups which are slower than direct look-ups.