diff --git a/src/constants.F90 b/src/constants.F90
index d75ee4e54a..4c1e3013c2 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -222,7 +222,7 @@ module constants
        N_3HEC  = 799, N_A0    = 800, N_AC   = 849, N_2N0   = 875, N_2NC   = 891
 
   ! Depletion reactions
-  integer, parameter :: DEPLETION_RX(6) = [N_2N, N_3N, N_4N, N_GAMMA, N_P, N_A]
+  integer, parameter :: DEPLETION_RX(6) = [N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N]
 
   ! ACE table types
   integer, parameter :: &
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 17880f6e1b..7a6b2444c7 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -887,7 +887,7 @@ contains
         micro_xs % reaction(:) = ZERO
 
         ! Only non-zero reaction is (n,gamma)
-        micro_xs % reaction(4) = sig_a - sig_f
+        micro_xs % reaction(1) = sig_a - sig_f
       end if
 
       ! Ensure these values are set
@@ -979,12 +979,24 @@ contains
 
       ! Depletion-related reactions
       if (need_depletion_rx) then
-        do j = 1, 6
-          ! Initialize reaction xs to zero
-          micro_xs % reaction(j) = ZERO
+        ! Initialize all reaction cross sections to zero
+        micro_xs % reaction(:) = ZERO
 
-          ! If reaction is present and energy is greater than threshold, set
-          ! the reaction xs appropriately
+        ! Physics says that (n,gamma) is not a threshold reaction, so we don't
+        ! need to specifically check its threshold index
+        i_rxn = this % reaction_index(DEPLETION_RX(1))
+        if (i_rxn > 0) then
+          associate (xs => this % reactions(i_rxn) % xs(i_temp))
+          micro_xs % reaction(1) = (ONE - f) * &
+               xs % value(i_grid - xs % threshold + 1) + &
+               f * xs % value(i_grid - xs % threshold + 2)
+          end associate
+        end if
+
+        ! Loop over remaining depletion reactions
+        do j = 2, 6
+          ! If reaction is present and energy is greater than threshold, set the
+          ! reaction xs appropriately
           i_rxn = this % reaction_index(DEPLETION_RX(j))
           if (i_rxn > 0) then
             associate (xs => this % reactions(i_rxn) % xs(i_temp))
@@ -992,12 +1004,17 @@ contains
                 micro_xs % reaction(j) = (ONE - f) * &
                      xs % value(i_grid - xs % threshold + 1) + &
                      f * xs % value(i_grid - xs % threshold + 2)
+              elseif (j >= 4) then
+                ! One can show that the the threshold for (n,(x+1)n) is always
+                ! higher than the threshold for (n,xn). Thus, if we are below
+                ! the threshold for, e.g., (n,2n), there is no reason to check
+                ! the threshold for (n,3n) and (n,4n).
+                exit
               end if
             end associate
           end if
         end do
       end if
-
     end if
 
     ! Initialize sab treatment to false
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index ef8b5915c1..b245bf5249 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -1074,17 +1074,17 @@ contains
         else
           ! Determine index in NuclideMicroXS % reaction array
           select case (score_bin)
-          case (N_2N)
-            m = 1
-          case (N_3N)
-            m = 2
-          case (N_4N)
-            m = 3
           case (N_GAMMA)
-            m = 4
+            m = 1
           case (N_P)
-            m = 5
+            m = 2
           case (N_A)
+            m = 3
+          case (N_2N)
+            m = 4
+          case (N_3N)
+            m = 5
+          case (N_4N)
             m = 6
           end select