When trying to optimize a function, one always should know what is the bottle-neck of this function - without you will spend a lot of time running in the wrong direction.

Let's use your python-function as baseline (actually I use result=np.zeros(shape,dtype=a.dtype) otherwise your method returns floats which is probably a bug):

>>> import numpy as np
>>> a=np.random.randint(1,1000,(300,300), dtype=np.int)
>>> b=np.random.randint(1,1000,(300,300), dtype=np.int)
>>> %timeit subtractPython(a,b)
274 ms ± 3.61 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The first question we should ask ourselves is: Is this task memory or CPU bound? Obviously, this is a memory-bound task - a subtraction is nothing compared to needed memory-read- and write-accesses.

This means, all above we have to optimize the memory layout in order to reduce cache-misses. As a rule of thumb, our memory accesses should access one consecutive memory address after another.

Is this the case? No, the array result is in C-order, i.e. row-major-order and thus the access

results[:, :, index] = subtracted

isn't consecutive. On the other hand,

results[index, :, :] = subtracted

would be a consecutive access. Let's change the way information is stored in result:

def subtract1(a, b):
    xAxisCount = a.shape[0]
    yAxisCount = a.shape[1]

    shape = (xAxisCount,  xAxisCount, yAxisCount) #<=== Change order
    results = np.zeros(shape, dtype=a.dtype)
    for index in range(len(b)):
        subtracted = (a - b[index])
        results[index, :, :] = subtracted   #<===== consecutive access
    return results

The timings are now:

>>> %timeit subtract1(a,b)
>>> 35.8 ms ± 285 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

There are also 2 more small improvements: we don't have to initialize result with zeros and we can save some python overhead, but this gives us just about 5%:

def subtract2(a, b):
    xAxisCount = a.shape[0]
    yAxisCount = a.shape[1]

    shape = (xAxisCount,  xAxisCount, yAxisCount) 
    results = np.empty(shape, dtype=a.dtype)        #<=== no need for zeros
    for index in range(len(b)):
        results[index, :, :] = (a-b[index])   #<===== less python overhead
    return results

>>> %timeit subtract2(a,b)
34.5 ms ± 203 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Now this is about factor 8 faster than the original version.

You could use Cython to try to speed-up this even further - but the task is probably still memory-bound, so don't expect to get it significantly faster - after all cython cannot make the memory work faster. However, without proper profiling it is hard to tell, how much improvement is possible - would not be surprised, if someone would come up with a faster version.