A very well described question. In mathematical terms we have two binomial populations:
$$x|npI\sim bin(n,p)$$
$$y|mqI\sim bin(m,q)$$
Where $x$ is the number of deaths in the control group with $n$ observations, $y$ is the number of deaths in the group of size $m$ which received the enhancer, and $I$ is the prior information - which includes the binomial assumption. The hypothesis is that $q>p$. The alternative which is most obvious is $q\leq p$ (i.e. we are not quesitoning the binomial model). I presume your "naive" approach is to consider $\hat{p}=\frac{x}{n}$ as fixed and then test if $q>\hat{p}$ using the cummulative binomial distribution.
To do the test and account for the uncertainty in $p$, the best approach is to integrate them out using a prior distribution which describes what is known about them (independently of what the data tells you about them). To do this we require the joint posterior for $p,q$.
$$p(pq|xynmI)\propto p(pq|I)p(xy|pqnmI)$$
Now I think it is reasonable to suppose the likehood factors, as you have assumed a binomial distribution, which is conditionally independent. So we have:
$$p(xy|nmpqI)\propto p^{x}(1-p)^{n-x}q^{y}(1-q)^{m-y}$$
Now because you have a large $n$ then the particular prior used won't matter unless it is quite strong. Assuming this is the case, then we can use a uniform prior $p(p|I)=1$. However we also require the conditonal prior $p(q|pI)$. This is likely to be non-uniform, but as before unless it is on the order of $m$ (i.e an $m$ degree polynomial) it won't be able to greatly affect the result. So i will use a flat prior also. This means the posterior is proportional to the likelihood which is is the product of two beta distributions. The probability we are after is given by:
$$Pr(q>p|xynmI)=\frac{
\int_0^1\int_p^1 q^{y}(1-q)^{m-y}dq p^{x}(1-p)^{n-x} dp
}{
\int_0^1\int_0^1 q^{y}(1-q)^{m-y} dq p^{x}(1-p)^{n-x} dp
}$$
$$=\frac{(n+1)!}{x!(n-x)!}
\int_0^1\left[\frac{(m+1)!}{y!(m-y)!}\int_0^{1-p} (1-q)^{y}q^{m-y}dq\right] p^{x}(1-p)^{n-x} dp
$$
The term inside the brackets can be converted to a sum via repeated use of integration by parts (see e.g. wikipedia), so we can write:
$$Pr(q>p|xynmI)=\frac{(n+1)!}{x!(n-x)!}
\int_0^1\left[\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!}p^{k}(1-p)^{m+1-k}\right] p^{x}(1-p)^{n-x} dp
$$
$$=\frac{(n+1)!}{x!(n-x)!}\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!}
\int_0^1p^{x+k}(1-p)^{n-x+m+1-k} dp
$$
$$=\frac{(n+1)!}{x!(n-x)!}\sum_{k=0}^{y}\frac{(m+1)!}{(m+1-k)!k!}\frac{(x+k)!(n+m+1-x-k)!}{(n+m+2)!}
$$
We can arrange this into 3 combinatorical functions (bit like hypergeometric pdf), and show how the uncertainty affects the answer as follows:
$$Pr(q>p|xynmI)=\frac{1}{{n+m+2 \choose n+1}}\sum_{k=0}^{y}{x+k \choose k}{n-x+m+1-k \choose m+1-k}$$
I think there is an intuitive way that one could reason this out as a good test from a frequentist point of view, but I can't think of it myself. For the "sanity check" note that it is a strictly increasing function of $y$ - any additional death in the enhanced group results in an increased probability for $q>p$. Additionally, because it is a Bayesian test, you can use it sequentially, and continue testing the enhanced cases until the probability is sufficiently high or low. Further, for large $a$ and fixed $b$ we have ${a+b\choose b}\approx \frac{a^b}{b!}$, and so as $n,x$ become large such that $\hat{p}=\frac{x}{n}\approx p$ we get:
$$Pr(q>p|xynmI)\approx\frac{(m+1)!}{n^{m+1}}\sum_{k=0}^{y}\frac{x^k}{k!}\frac{(n-x)^{m+1}}{(m+1-k)!}=\sum_{k=0}^{y}{m+1\choose k}\hat{p}^k(1-\hat{p})^{m+1-k}$$
Which is basically the "naive" approach just by plugging in the estimate (and another nice equivalence between bayes and frequentist methods). For your case where $n=300$ using open office on my pc had no problem calculating the choose functions with $m=300$ also (even though it can't calculate above $170!$). So there isn't any need to use the approximation, apart from gaining an intuitive understanding the result. They are both simple sums.
When $n>>m$ and $x>>y$, then the binomial plug-in method will work fine. When $m$ is also large, you are probably better off using the normal approximation to the binomial sum above, and then you get the standard test for difference in normal means from two independent populations - just compare the difference in proportions to the standard deviation of the difference.
Note that if you do wish to use a non-uniform prior for the test, a computationally easy way is to assign two independent beta distributions, so we have $p(pq|I)\propto p^{a_p-1}(1-p)^{b_p-1}q^{a_q-1}(1-q)^{b_q-1}$ (with all parameters integers) and the test becomes:
$$\frac{1}{{n+m+a_p+b_p+a_q+b_q-2 \choose n+a_p+b_p-1}}\sum_{k=0}^{y+a_q-1}{x+k+a_p-1 \choose k}{n-x+m-k+b_p+a_q+b_q-2 \choose m-k+a_q+b_q-1}$$
which is effectively adding pseudo observations to your data set.
UPDATE
In response to the comments, it is quite simple to write an R function to calculate the above probability as follows:
prob.q.gt.p <- function(x,y,n,m,a_p=1,b_p=1,a_q=1,b_q=1){
prob <- 0
sum.to <- y+a_q-1
for(k in 0:sum.to){
prob <- prob + choose(x+k+a_p-1,k) * choose(n+m+b_p+a_q+b_q-2-x-k,m+a_q+b_q-1-k)
}
prob <- prob / choose(n+m+a_p+b_p+a_q+b_q-2,n+a_p+b_p-1)
return(prob)
}
#example run - should return value of 0.0336
prob.q.gt.p(20,0,300,50)
Now I don't know what observed values you intend to use, except that $n=300$. Now suppose $x=20$ ($20$ deaths) for death rate of $6.7\text{%}$, and also that $m=50$ (expose $50$ new organisms to increased toxin). The table below shows the probability that $q>p$ for various values of $y$
$$\begin{array}{c|c}y & \frac{y}{m}\times 100\text{%}& Pr(q>p|y,x=20,n=300,m=50,I)\\ 0 & 0\text{%} & 0.0336 \\ 1 & 2\text{%} & 0.1424 \\ 2 & 4\text{%} & 0.3238\\ 3 & 6\text{%} & 0.5308 \\ 4 & 8\text{%} & 0.7126 \\ 5 & 10\text{%} & 0.8433 \\ 6 & 12\text{%} & 0.9231 \\ 7 & 14\text{%} & 0.9658 \\ 10 & 20\text{%} & 0.9982 \\ 15 & 30\text{%} & >0.9999
\end{array}$$
